Female Anterior Cruciate Ligaments Exhibit a Muted Mechanobiological Response to Mechanical Loading

TL;DR

Female ACLs respond less to mechanical stress compared to male ACLs, which may explain why women are more prone to ACL tears.

Contribution

This study identifies muted mechanotransduction pathways in female ACLs and explores the role of estrogen in tissue remodeling.

Findings

Female ACLs showed fewer differentially expressed genes in response to cyclic loading compared to males.

Mechanotransduction pathways were enriched only in male ACLs after loading.

Estrogen treatment did not significantly alter tissue remodeling marker expression in ACL explants.

Abstract

Female athletes are significantly more likely to tear their anterior cruciate ligament (ACL) compared to their male counterparts. While there are several potential reasons for this, previous data from our lab demonstrated that female ACL explants have an impaired remodeling response to loading, which may prevent the repair of fatigue damage and lead to increased ACL rupture. The objective of this study was to identify the biological mechanisms driving the impaired remodeling of female ACLs to cyclic loading, including the role of estrogen. ACLs were harvested from male and female New Zealand white rabbits and cyclically loaded in a tensile bioreactor, followed by bulk RNA‐sequencing. Additional ACL explants treated with or without estradiol were analyzed using RT‐qPCR to determine the regulatory effect of estrogen on markers for tissue remodeling and inflammatory cytokines with cyclic loading. We found that female ACLs exhibited significantly fewer differentially expressed genes in response to loading compared to male ACLs. Additionally, multiple mechanotransduction pathways were enriched with loading only in the male ACLs. While a few estrogen‐related pathways were enriched in both male and female ACLs with loading, the expression of tissue remodeling markers was not different between estrogen treatment and vehicle control. Together, our findings highlight specific mechanotransduction pathways that may be responsible for the muted biological response of female ACLs to load, which provides a potential explanation for the increased rate of ACL tears in women.

Figures

Differentially expressed genes. (A) Response of male and female ACLs in response to cyclic loading and comparing the response between sexes (sex effect load) (n = 3 males, n = 2 females). Male and female loads are relative to their sex‐matched static load. Sex effect load is comparing the female ACL response to load to the male ACL response to load. (B) Comparison of baseline gene expression differences of female ACLs relative to male ACLs (n = 3 static, n = 2 fresh). All DEGs had an absolute log2FC > 1 and an adjusted p‐value < 0.05.

Basal remodeling gene expression response to estrogen treatment. RT‐qPCR analysis of statically loaded estrogen‐treated ACLs relative to their sex‐matched ethanol vehicle control (represented by the red box and whisker plots). (A) Statistically loaded male ACL response. (B) Statistically loaded female ACL response. (C) Comparison of statically loaded male and female ACL response to estrogen treatment (n = 6 for all male samples, n = 7 for female vehicle, n = 5–6 for female estrogen‐treated). Data represented as a box and whisker plot with the whiskers representing the min and max data.

Effect of estrogen treatment on ACL remodeling in response to load. RT‐qPCR analysis of cyclically loaded male and female ACLs relative to their sex and treatment‐matched 0.1 MPa static load (represented by the red box and whisker plots). (A) Gene expression changes in vehicle control male ACLs. (B) Gene expression changes in estrogen‐treated male ACLs. (C) Gene expression changes in vehicle control female ACLs. (D) Gene expression changes in female estrogen‐treated ACLs (n = 6 for all male samples, n = 7 for female vehicle, n = 6 for female estrogen‐treated). Data represented as a box and whisker plot with the whiskers representing the min and max data.

Basal estrogen receptor target genes response to estrogen treatment. RT‐qPCR analysis of statically loaded estrogen‐treated ACLs relative to their sex‐matched ethanol vehicle control (represented by the red box and whisker plots). (A) Statistically loaded male ACL response. (B) Statistically loaded female ACL response. (C) Comparison of statically loaded male and female ACLs response to estrogen treatment (n = 5 for male vehicle, n = 6 for male estrogen treated, n = 6–7 for female vehicle, n = 6 for female estrogen treated). Data represented as a box and whisker plot with the whiskers representing the min and max data.

Effect of estrogen treatment on estrogen receptor target genes in response to load. RT‐qPCR analysis of cyclically loaded male and female ACLs relative to their sex and treatment‐matched 0.1 MPa static load. (A) Gene expression changes in vehicle control male ACLs. (B) Gene expression changes in estrogen‐treated male ACLs. (C) Gene expression changes in vehicle control female ACLs. (D) Gene expression changes in female estrogen‐treated ACLs (n = 6 for all male samples, n = 6–7 for female vehicle, n = 6 for female estrogen‐treated). Data represented as box and whisker plot with the whiskers representing the min and max data.

Tables

Table 1. List of top 20 DEGs in response to load in male ACLs (and the sex effect of load (female ACLs response to load relative to male ACLs response to load).

Male load			Sex effect load
Gene ID	Log2FC	p adj	Gene ID	Log2FC	p adj
ALK	−6.791	0.04	BARHL1	6.258	0.02
COX6	−6.299	0.01	SOSTDC1	−6.100	0.02
ENSOCUG00000013331	5.807	0.01	COL21A1	6.073	0.04
DSG2	−5.457	0.04	IL1RL1	−5.616	0.01
MASP1	−5.340	0.03	NALCN	−5.561	0.01
COL21A1	−5.169	0.02	CDO1	−5.527	0.02
NEK2	−5.097	0.02	GBP4	5.455	0.02
ENSOCUG00000023797	−5.060	0.02	CYTL1	−5.116	0.01
NPS	4.549	0.04	OLFM3	4.759	0.03
MANEAL	4.402	0.04	ANGPTL7	−4.747	0.02
IL1RL1	4.353	0.01	GALNT15	−4.625	0.02
IL2RA	−4.289	0.02	CADM4	4.375	0.05
RSPO1	−4.214	0.03	ENSOCUG00000029275	4.334	0.04
CAPN3	−4.209	0.03	ENSOCUG00000029534	−4.116	0.02
SPON1	−4.151	0.02	CCDC3	−3.965	0.01
KIF23	−4.130	0.03	MYOC	−3.838	0.02
PENK	−4.075	0.05	ADAMTSL3	−3.733	0.05
RASL12	−4.013	0.04	PRELP	−3.687	0.01
APLN	−4.004	0.04	NOVA1	−3.663	0.02
NALCN	3.958	0.01	MAMDC2	−3.595	0.01

Table 2. Differential expression of ECM remodeling genes in all comparisons.

	Male load		Female load		Sex effect load		Sex effect static		Sex effect fresh
Gene name	Log2FC	p adj	Log2FC	p adj	Log2FC	p adj	Log2FC	p adj	Log2FC	p adj
COL1A1	0.13	0.96	−0.86	0.50	−0.99	0.58	0.98	0.51	0.56	0.89
COL1A2	0.26	0.85	−0.91	0.36	−1.17	0.40	0.95	0.40	1.22	0.60
LOX	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
COL3A1	−0.04	0.98	−1.62	0.17	−1.58	0.28	0.65	0.62	1.19	0.60
TGFB1	0.06	0.92	−0.06	0.91	−0.12	0.85	0.10	0.86	−0.08	0.95
ACTA2	−2.40	0.11	−0.22	0.88	2.18	0.24	−0.72	0.59	−1.22	0.61
TIMP1	0.55	0.47	−0.54	0.40	−1.09	0.25	0.12	0.92	−0.62	0.59
TIMP3	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
MMP1	−2.27	0.11	0.96	0.49	3.23	0.12	−1.12	0.49	0.38	0.90
MMP2	−1.62	0.03	−0.27	0.66	1.36	0.14	−0.82	0.23	0.29	0.85
MMP10	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
MMP13	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
IL1B	−2.13	0.06	−0.53	0.56	1.60	0.25	−1.22	0.23	−0.65	0.94
PTGS2	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA

Table 3. Differential expression of estrogen receptor target genes in all comparisons.

	Male load		Female load		Sex effect load		Sex effect static		Sex effect fresh
Gene name	Log2FC	p adj	Log2FC	p adj	Log2FC	p adj	Log2FC	p adj	Log2FC	p adj
ESR1	3.542	0.02	−0.662	0.61	−4.204	0.05	2.386	0.16	0.339	0.81
IGFBP4	−1.989	0.10	−0.750	0.49	1.238	0.42	−0.583	0.62	−1.360	0.52
TPD52L1	2.396	0.10	−1.047	0.61	−3.444	0.17	0.158	0.96	1.560	0.49
GPER1	−2.192	0.12	0.360	0.75	2.552	0.16	−0.620	0.59	−0.213	0.96
STC2	−1.870	0.35	1.350	0.37	3.219	0.20	−0.188	0.95	0.348	0.95
RBBP8	−0.212	0.70	0.677	0.16	0.889	0.17	−0.031	0.97	−0.097	0.97
SIAH2	0.944	0.18	−0.025	0.98	−0.968	0.29	0.587	0.41	1.010	0.40
ADCY9	0.249	0.70	−0.201	0.76	−0.450	0.56	−0.020	0.98	−0.661	0.45
SERPINE1	−1.830	0.08	0.259	0.82	2.089	0.15	−0.706	0.55	0.810	0.61
EPHA4	−1.150	0.29	−0.348	0.73	0.802	0.55	−0.059	0.97	−0.906	0.60

Table 4. Mechanotransduction pathways enriched in response to load.

Male load		Sex effect load
Pathway	Z‐score	Pathway	Z‐score
Integrin cell surface interactions	−2	FAK signaling	−2
G‐protein coupled receptor signaling	−2	Gap junction signaling	−1
Class A/1 (Rhodopsin‐like receptors)	−1	WNT/β‐catenin signaling	0
RHO GTPase cycle	−1	RHO GTPase cycle	0
Extra‐nuclear estrogen signaling	−1	Integrin signaling	NA
cAMP‐mediated signaling	−1	Actin nucleation by ARP‐WASP complex	NA
FAK signaling	−1	Integrin cell surface interactions	NA
WNT/SHH axonal guidance signaling pathway	−1	Signaling by MET	NA
RAF/MAP kinase cascade	0	PI3K/AKT signaling	NA
RHO GTPases activate CIT	NA
Signaling by ALK	NA
Signaling by NOTCH3	NA
Signaling by PDGF	NA
Signaling by MET	NA
PI3K/AKT signaling	NA

Table 5. Mechanotransduction pathways enriched in both male ACL response to load and sex effect load.

Pathway	Male load Z‐score	Sex effect load Z‐score
FAK signaling	−1	−2
Integrin cell surface interactions	−2	NA
RHO GTPase cycle	−1	0
Signaling by MET	NA	NA
PI3K/AKT signaling	NA	NA

Table 6. Extracellular matrix pathways enriched in response to load.

Male load		Sex effect load
Pathway	Z‐score	Pathway	Z‐score
Collagen degradation	−2	Elastic fiber formation	NA
Inhibition of matrix metalloproteases	2	Collagen biosynthesis and modifying enzymes	NA
Degradation of the extracellular matrix	−2	Collagen chain trimerization	NA
Extracellular matrix organization	−2
Collagen biosynthesis and modifying enzymes	−2
Collagen chain trimerization	−2
Assembly of collagen fibrils and other multimeric structures	NA
Activation of matrix metalloproteinases	NA

Table 7. Extracellular matrix pathways enriched in both male ACLs response to load and sex effect load.

Pathway	Male load Z‐score	Sex effect load Z‐score
Collagen biosynthesis and modifying enzymes	−2	NA
Collagen chain trimerization	−2	NA

Table 8. List of top 20 predicted upstream regulators of response to load in female and male ACLs (relative to sex‐matched static control) and the sex effect of load (female ACL response to load relative to male ACL response to load).

Male load		Female load		Sex effect load
Upstream regulator	Activation z‐score	Upstream regulator	Activation z‐score	Upstream regulator	Activation z‐score
TNF	−4	DMPK	NA	ARID1A	−2
TGFB1	−3			cycloheximide	−2
Lipopolysaccharide	−3			AR	−2
VEGF	−3			dexamethasone	−2
CSF2	−3			anisomycin	−1
IL1B	−3			TGFB1	−1
AGT	−2			TNF	1
ESR1	−2			PLAU	1
IGF1	−2			SORL1	1
Demethasone	2			FAM3C	0
PEAR1	3			CDCA8	0
GUCY1B1	−1			PDGF‐BB (complex)	0
SOLR1	−1			LAMA4	0
Prostaglandin E2	−1			NR1I3	NA
ERBB2	−1			GNB2	NA
PTEN	1			UBE3D	NA
AR	0			BRMS1	NA
Immunoglobulin (complex)	0			sesaminol	NA
Beta‐estradiol	0			BCL9/CREBBP/EP300/CTNNB1/LEF/TCF (complex)	NA
Dihydrotestosterone	0			candesartan cilexetil	NA

Full text

Introduction

1

While anterior cruciate ligament (ACL) tears are common in general during athletic activity [1], female athletes are 2–8 times more likely to tear their ACLs compared to males [2, 3, 4]. Potential reasons for this include differences in anatomy, sex hormones, body size, and neuromuscular control [5, 6, 7]. However, women may also be more prone to ACL injury because of sex‐specific ACL remodeling in response to fatigue loading. Recent studies suggest that ACL injuries occur from fatigue failure and accumulation of tissue damage in the ACL [8, 9, 10, 11, 12, 13]. Previous data from our lab demonstrated that female rabbit ACL explants have an impaired remodeling response to load [14]. Specifically, we found that female ACLs fail to increase anabolic gene expression with loading and have increased catabolic gene expression compared to male ACL explants, which may prevent the repair of ACL fatigue damage and may explain the increased rate of ACL tears in women. However, the cause of this sex difference in ACL mechanobiology is unknown.

Differences in sex hormones (e.g., estrogen) are one obvious potential reason for the sex‐dependent ACL response to loading. Estrogen receptors (ERs) are present in both male and female ACL fibroblasts [15]. Additionally, estrogen influences multiple genes responsible for soft tissue remodeling [16, 17, 18, 19, 20, 21, 22, 23]; however, the magnitude and directionality of this effect on the ACL are inconsistent [16, 18, 19, 24, 25, 26, 27, 28, 29, 30, 31]. Despite these discrepancies, clinical data demonstrate that there are substantial changes in knee laxity during the menstrual cycle (average increase of 3 mm) that are correlated with estradiol concentrations [32]. Given that a 1.3 mm asymmetry in knee laxity is associated with a fourfold increase in ACL injury risk [33], these hormonal changes in knee laxity may help explain the increased ACL injury risk in women. Furthermore, changes in knee laxity can be attenuated through oral contraceptives that suppress spikes in estrogen levels [34, 35]. Together, these data suggest that estrogen has a direct effect on ACL tissue mechanics and may contribute to ACL rupture.

Previous data also suggest that estrogen regulates the response of cells to mechanical loading. A prior study on ACL fibroblasts found that mechanical loading completely reverses the effect of estrogen on gene expression [18]. However, the mechanisms underlying this interplay between estrogen and mechanical loading are unknown. In other tissues, ER signaling has been shown to interact with multiple mechanotransduction pathways (i.e., ERK, YAP/TAZ) [36, 37, 38, 39, 40, 41, 42, 43, 44]. However, the directionality of this effect is inconsistent. Specifically, studies on cancer cells show that ER signaling both increases [42, 43, 44] and decreases [39, 40, 41] YAP activation. While previous studies demonstrate that ER signaling regulates cellular mechanotransduction, it is not possible to predict the effect of estrogen on ACL fibroblast mechanotransduction from prior data alone.

This study builds upon our prior findings in which we demonstrated that female ACLs have an impaired remodeling response to mechanical loading compared to male ACLs [14]. While those findings exhibited sex‐based differences in ACL mechanobiology, the underlying signaling mechanisms remain unclear. Therefore, the objective of this study was to identify potential biological mechanisms driving the impaired remodeling response of female ACLs. We conducted bulk RNA‐sequencing on cyclically loaded male and female rabbit ACL explants to determine the innate sex differences in ACL mechanobiology in the native tissue environment. Additionally, we investigated the effect of estrogen treatment on the expression of remodeling genes in male and female loaded ACL explants. Our overall hypothesis was that female ACLs would exhibit an impaired mechanobiological response to load (increased catabolism, decreased anabolism), and this response would be influenced by estrogen. Specifically, we hypothesized that female ACLs would have altered mechanotransduction (suppression of ERK and YAP) as well as increased ER activity in response to load compared to males. We also hypothesized that estrogen treatment would increase catabolic and inflammatory gene expression in both male and female ACL explants in response to load.

Materials and Methods

2

ACL Harvest

2.1

A total of 33 male and female white New Zealand rabbits (2.8–3.2 kg, 14–16 weeks old, Charles River) were euthanized (Supporting Information S1: Table 1), and ACLs were isolated under sterile conditions from a study approved by the Institutional Animal Care and Use Committee at Penn State University. This age range coincides with the onset of sex differences in serum estrogen levels in rabbits [45]. ACLs were harvested following a previously published protocol [14]. Briefly, all surrounding tissues were removed, and the femur and tibia were cut into bone blocks to grip during mechanical loading. Cross‐sectional area was determined by assuming an elliptical cross‐section, and major and minor diameters were determined using calipers.

Mechanical Loading

2.2

ACLs were placed in a custom tensile bioreactor [46] with culture media (low‐glucose DMEM, 2% penicillin‐streptomycin, 5% FBS, 25 mm HEPES, 4 mm GlutaMAX, 1 mm 2‐phospho‐L‐ascorbic acid trisodium) and kept at 37°C and 5% CO_2_. A subset of the ACLs was untreated to determine the innate mechanobiology of male and female ACL explants. To investigate the effect of estrogen, ACLs were treated with 10 nM 17β‐estradiol (E2) (Sigma) or ethanol (0.1% v/v) as a vehicle control. The E2 concentration was based on a previous study investigating the role of estrogen and cyclic loading on isolated ACL fibroblasts [18]. Since underloading is detrimental to tendon/ligament homeostasis [47, 48], a 0.1 MPa static load was applied for 18 and 20 h to acclimate the untreated and treated (vehicle/estrogen) ACLs to culture conditions, respectively. The 18‐h timepoint was chosen to replicate previous work that found differences in the expression of remodeling genes with cyclic loading [14] for the RNA‐sequencing, and the 20‐h timepoint was chosen to be more consistent with prior studies of estrogen treatment on ACL mechanobiology [18]. The bioreactor then cyclically loaded the samples to 4 MPa at 0.5 Hz for 8 h. This stress level was chosen since it represents the stress that exhibited the largest gene expression differences between male and female ACLs [14]. Control samples were maintained under a 0.1 MPa static load for the same duration. Previous data from our lab demonstrated that these loading parameters maintain cell viability in the explants [14, 46].

RNA Extraction

2.3

Immediately following loading, samples were removed, and RNA was extracted following a previously published protocol [14]. Briefly, ACLs were excised from the bone blocks and were rinsed with ice‐cold RNase‐free water and flash frozen in liquid nitrogen. The tissue was pulverized (6775 Freezer/Mill, SPEX) and total RNA was extracted with RNeasy minicolumns (RNeasy Fibrous Tissue Kit, Qiagen). RNA concentration was quantified using a Qubit 4 Fluorometer (Thermo Fisher), and RNA integrity numbers (RINs) were determined using a TapeStation (Agilent).

RNA‐Sequencing

2.4

Samples that passed quality control (RIN > 6, total RNA > 250 ng) were sent to GENEWIZ (South Plainfield, NJ) for bulk RNA‐Sequencing (Supporting Information S1: Table 1). cDNA libraries were prepared using the NEBNext Ultra II RNA Library Prep for Illumina using the manufacturer's instructions (NEB, Ipswich). The sequencing libraries were validated on the Agilent TapeStation (Agilent Technologies) and quantified using a Qubit 3.0 Fluorometer (Invitrogen) as well as qPCR (KAPA Biosystems). The sequencing libraries were clustered on a lane of a NovaSeq. 6000 S4 flow cell. The samples were sequenced with an average of 32 million reads per sample using a 2 × 150 bp paired‐end configuration. Image analysis and base calling were conducted by the Control software. Raw sequence data (.bcl files) generated by the sequencer were converted into fastq files and de‐multiplexed using Illumina's bcl2fastq 2.17 software. One mismatch was allowed for index sequence identification. Sequenced reads were trimmed to remove possible adapter sequences and nucleotides with poor quality using Trimmomatic v.0.36. The trimmed reads were mapped to the OryCun2.0 reference genome available on ENSEMBL using the STAR aligner v.2.5.2b. Unique gene hit counts were calculated by using FeatureCount from the Subread package v.1.5.2.

Differential Gene Expression Analysis

2.5

Analysis and visualization of the bulk RNA‐sequencing was performed utilizing the NIH Integrated Analysis Platform (NIDAP) bulk RNA‐sequencing pipeline that utilizes R programs developed by a team of NCI bioinformaticians (Plantair Technologies, Denver, Colorado). Raw hit counts were uploaded, and quality control measures were performed using principal component analysis and supervised hierarchical clustering to determine intergroup relationships and identify outliers. A female‐loaded sample (FL3) was identified as an outlier and removed (Supporting Information S2: Figure 1). After outlier removal and normalization, differentially expressed gene (DEG) analysis was conducted utilizing the Limma Voom R package. Five comparisons were made to compare the response of each sex to load (female load, male load) and the relationships between sex in the fresh (sex effect fresh), static (sex effect static), and loaded (sex effect load) samples. The static control and male samples were used as the baseline for the comparisons. A moderated t‐statistic was used to determine p‐values, and a Benjamini–Hochberg false discovery rate was used to obtain the adjusted p‐values. Genes with an adjusted p‐value < 0.05 and |log2 fold change| > 1 were identified as DEGs. A list of DEGs was extracted from the model for each of the five comparisons.

Pathway Analysis and Upstream Regulator Analysis

2.6

After differential gene expression analysis, all DEGs for each comparison were uploaded into Ingenuity Pathway Analysis (IPA, Qiagen, Version Q4 2024). Significantly enriched pathways were determined with a −log(p‐value) > 1.3. Additionally, upstream regulator analysis was performed to determine potential regulators that could explain the observed changes in gene expression. To identify significantly activated or inhibited pathways or upstream regulators, a z‐score was calculated in IPA for each pathway or upstream regulator. A positive z‐score suggests inhibition, while a negative z‐score suggests activation. The magnitude of the z‐score reflects the confidence of the predicted activation or inhibition state, with |z‐scores | > 2 considered a strong prediction and |z‐scores| > 1 considered a moderate prediction.

Reverse Transcription Quantitative Polymerase Chain Reaction (RT‐qPCR)

2.7

The RNA extracted from each sample was diluted to 1 ng/ul, and cDNA was synthesized (High‐Capacity cDNA Reverse Transcription Kit with RNase Inhibitor, Thermo Fisher). qPCR was performed using TaqMan probes (Supporting Information S1: Table 2) and a StepOne Plus Real‐Time PCR system to measure the expression of anabolic (COL1A1, COL1A2, LOX, COL3A1, TGFβ1, ACTA, TIMP1, TIMP3), catabolic (MMP1, MMP2, MMP10, MMP13), and inflammatory (IL‐1β, PTGS2) genes. Additionally, ER target genes (ESR1, GPER1, PGR) were measured along with GAPDH as the reference gene. PCR efficiency and cycle number quantification (Cq) were obtained for each individual reaction using PCR‐Miner (version 4.0) [49]. Samples that exhibited undetectable fluorescence for a given gene were excluded from the analysis for that gene. An outlier test (ROUT) was conducted to determine if any of the reaction efficiencies for each gene were an outlier. After removal of outliers, a single efficiency value was calculated for each probe by averaging all the individual sample efficiencies. Gene expression was quantified using the delta‐delta Cq method, correcting for primer efficiencies [50]. The fold change for each sample was calculated as:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

where the probe efficiencies are listed in Supporting Information S1: Table 2.

Statistical Analysis of RT‐qPCR Data

2.8

The fold change for each gene was log‐normalized, and Mann–Whitney tests were conducted on the normalized data to compare the gene expression between treatments and sex. To determine the effect of mechanical loading on treatment, the estrogen‐treated and ethanol vehicle control samples were compared to their treatment‐matched statically loaded control samples (for both males and females). To determine baseline treatment‐specific gene expression differences, the estrogen‐treated static control samples were compared to the vehicle control static control samples (for both males and females). To determine baseline sex‐specific gene expression differences, the female vehicle control samples were compared to the male vehicle control samples. Additionally, Mann–Whitney tests were used on the transformed data to determine the differential expression of each gene compared to the respective control condition. Given the numerous statistical comparisons, p‐values were corrected using a false discovery rate analysis within each gene category (anabolic, catabolic, inflammatory, and ER target genes) with a Q‐value of 5%. Statistical significance after corrections was set as p < 0.05, with statistical trends set at p < 0.10. All statistical analysis was performed using GraphPad Prism (version 10.3.1).

Results

3

Differential Gene Expression Analysis

3.1

In response to cyclic loading, male ACLs had 259 DEGs (122 upregulated, 137 downregulated) while female ACLs only had 10 DEGs (4 upregulated, 6 downregulated). When comparing the response to load between female and male ACLs (sex effect load), there are 116 DEGs (25 upregulated, 91 downregulated in female ACLs compared to males) (Figure 1A, Supporting Information S1: Tables 3–5). When looking at the top 20 DEGs in the male samples' response to load, 75% of the DEGs were downregulated. Similarly, for the sex effect load comparison, 70% of the DEGs were further downregulated in female ACLs in response to load compared to males (Table 1). These differences are despite the fact that there are minimal differences between female and male ACLs before loading. Specifically, there were only 40 DEGs (19 upregulated, 21 downregulated) in female static ACLs compared to male samples, and only 2 DEGs (1 upregulated, 1 downregulated) in female freshly harvested ACLs compared to males (Figure 1B, Supporting Information S1: Tables 6 and 7).

Differentially expressed genes. (A) Response of male and female ACLs in response to cyclic loading and comparing the response between sexes (sex effect load) (n = 3 males, n = 2 females). Male and female loads are relative to their sex‐matched static load. Sex effect load is comparing the female ACL response to load to the male ACL response to load. (B) Comparison of baseline gene expression differences of female ACLs relative to male ACLs (n = 3 static, n = 2 fresh). All DEGs had an absolute log2FC > 1 and an adjusted p‐value < 0.05.

Genes indicative of tissue remodeling (anabolic, catabolic, and inflammatory markers) [14] were searched in the list of DEGs for all five comparisons. In response to load, male ACLs downregulated catabolic gene MMP2 with inflammatory marker IL1β trending toward downregulation, while female ACLs exhibited no differentially expressed remodeling genes (Table 2). When comparing across sex, there were no differentially expressed remodeling genes for any of the comparisons. (Table 2).

Differentially Expressed ER Target Genes

3.2

Previously reported ER target genes [51, 52] were searched in the list of DEGs for all five comparisons. In response to load, male ACLs upregulated ESR1 and multiple genes trended towards downregulation (IGFBP4, TPD52L1, SERPINE1) while female ACLs exhibited no differentially expressed ER target genes (Table 3). When comparing the DEGs across sex, there were surprisingly no differentially expressed ER target genes for the freshly harvested or static samples. In response to load, ESR1 trended towards downregulation in female‐loaded ACLs compared to male‐loaded ACLs (Table 3).

Pathway Analysis of ACLs Response to Cyclic Load

3.3

In response to cyclic loading, IPA found that male ACLs had 139 signaling pathways enriched, while female ACLs only had 31 enriched pathways (Supporting Information S1: Tables 8 and 9). When directly comparing across sex (sex effect load), there were 82 pathways enriched (Supporting Information S1: Table 10).

We then investigated which estrogen‐related and mechanotransduction pathways were enriched with loading. For estrogen, extra‐nuclear estrogen signaling was enriched in male ACLs in response to load and was trending towards inhibition (z‐score: −1), while ESR‐mediated signaling was enriched in female ACLs in response to load (Supporting Information S1: Tables 8 and 9).

For mechanotransduction, male ACLs had 15 pathways enriched in response to load (Table 4) with 8 pathways predicted to be downregulated and zero pathways predicted to be upregulated (Table 4). The female ACLs exhibited no enriched pathways. Additionally, there were 9 pathways enriched when directly comparing across sex (sex effect load). Prominent pathways of interest included integrin signaling, FAK signaling, MAPK signaling, PI3K/AKT signaling, Rho GTPase cycle, and G‐protein coupled receptor signaling. When comparing enriched pathways between male and female ACLs in response to load, five pathways were enriched in both male load and sex effect of load, which were generally predicted to be inhibited (Table 5). Interestingly, FAK signaling was predicted to be inhibited in the sex effect load comparison, indicating an increased inhibition of FAK in female ACLs compared to male ACLs in response to load.

Finally, we also investigated pathways associated with extracellular matrix (ECM) organization. Male ACLs had 8 ECM pathways enriched in response to load (Table 6), while female ACLs had no enriched pathways. Three pathways were enriched when directly comparing across sex (sex effect load). These included collagen synthesis and degradation, activation and inhibition of matrix metalloproteinases, and elastic fiber formation. Interestingly, most of these pathways were predicted to be inhibited in the male ACLs' response to load. When looking at what pathways are different between male ACLs and female ACLs in response to load, two pathways related to collagen synthesis were enriched in both the male load and sex effect of load comparisons (Table 7).

Upstream Regulator Analysis of ACLs Response to Cyclic Load

3.4

IPA upstream regulator analysis identified several key regulators of male ACL response to load (Table 8), including multiple regulators related to sex hormones (beta‐estradiol, ERBB2, dihydrotestosterone, AR, ESR1), growth factors (TGFβ1, VEGF, IGF1) and inflammatory cytokines (TNF, LPS, IL1β, prostaglandin E2), with the majority of these regulators predicted to be inhibited. Only DMPK was predicted to be a regulator of the female ACLs' response to load (Table 8). In the sex effect load, AR and TGFβ1 were predicted to be inhibited, while TNF was predicted to be activated.

Effects of Estrogen Treatment on the Remodeling Response of ACLs to Cyclic Load

3.5

When comparing the expression of remodeling genes between male and female statically loaded ACLs cultured without estrogen (vehicle control), there were multiple anabolic (COL1A2, LOX, TGFβ1, TIMP3), catabolic (MMP1, MMP2, MMP13), and inflammatory (IL1β, PTGS2) markers significantly upregulated in the female ACLs compared to males and MMP10 was trending toward upregulation (Supporting Information S2: Figure 2A). However, there was no differential expression of any of the anabolic, catabolic, or inflammatory markers with estrogen treatment in the statically loaded male or female ACLs (Figure 2A,B). This was also true when comparing the response of estrogen treatment between male and female statically loaded ACLs (Figure 2C).

Basal remodeling gene expression response to estrogen treatment. RT‐qPCR analysis of statically loaded estrogen‐treated ACLs relative to their sex‐matched ethanol vehicle control (represented by the red box and whisker plots). (A) Statistically loaded male ACL response. (B) Statistically loaded female ACL response. (C) Comparison of statically loaded male and female ACL response to estrogen treatment (n = 6 for all male samples, n = 7 for female vehicle, n = 5–6 for female estrogen‐treated). Data represented as a box and whisker plot with the whiskers representing the min and max data.

No differential gene expression was observed with cyclical loading in the male ACLs compared to static controls with vehicle or estrogen treatment (Figure 3A,B). Direct comparisons between the vehicle control and estrogen‐treated male samples found that there was no effect of estrogen on gene expression (Supporting Information S2: Figure 3A). Consistently, the female cyclically loaded ACLs exhibited no differential gene expression compared to static controls with either vehicle or estrogen treatment (Figure 3C,D). Direct comparisons between vehicle and estrogen‐treated cyclically loaded female samples found no effect of estrogen on gene expression (Supporting Information S2: Figure 3B). Finally, direct comparison between the male and female cyclically loaded ACLs found that there were no sex differences in the vehicle or estrogen‐treated groups (Supporting Information S2: Figure 3C,D).

Effect of estrogen treatment on ACL remodeling in response to load. RT‐qPCR analysis of cyclically loaded male and female ACLs relative to their sex and treatment‐matched 0.1 MPa static load (represented by the red box and whisker plots). (A) Gene expression changes in vehicle control male ACLs. (B) Gene expression changes in estrogen‐treated male ACLs. (C) Gene expression changes in vehicle control female ACLs. (D) Gene expression changes in female estrogen‐treated ACLs (n = 6 for all male samples, n = 7 for female vehicle, n = 6 for female estrogen‐treated). Data represented as a box and whisker plot with the whiskers representing the min and max data.

Effects of Estrogen Treatment on Estrogen Target Genes of ACLs to Cyclic Load

3.6

In contrast to our RNA‐sequencing data (sex effect static), we found GPER1 to be significantly upregulated in the statically loaded female ACLs compared to males without estrogen treatment (Supporting Information S2: Figure 2B). However, there was no differential expression of any of the ER genes in the statically loaded male or female ACLs in response to estrogen treatment (Figure 4A,B). Additionally, there were no differences in gene expression when comparing the response of estrogen treatment between male and female statically loaded ACLs (Figure 4C).

Basal estrogen receptor target genes response to estrogen treatment. RT‐qPCR analysis of statically loaded estrogen‐treated ACLs relative to their sex‐matched ethanol vehicle control (represented by the red box and whisker plots). (A) Statistically loaded male ACL response. (B) Statistically loaded female ACL response. (C) Comparison of statically loaded male and female ACLs response to estrogen treatment (n = 5 for male vehicle, n = 6 for male estrogen treated, n = 6–7 for female vehicle, n = 6 for female estrogen treated). Data represented as a box and whisker plot with the whiskers representing the min and max data.

Consistent with our findings with the markers of tissue remodeling, there was no differential expression of the estrogen target genes with cyclic loading in the male (Figure 5A,B) or female samples (Figure 5C,D) with or without estrogen treatment. Similarly, there was no effect of estrogen (Supporting Information S2: Figure 4A,B) or sex (Supporting Information S2: Figure 4C,D) on the expression of estrogen target genes in response to loading.

Effect of estrogen treatment on estrogen receptor target genes in response to load. RT‐qPCR analysis of cyclically loaded male and female ACLs relative to their sex and treatment‐matched 0.1 MPa static load. (A) Gene expression changes in vehicle control male ACLs. (B) Gene expression changes in estrogen‐treated male ACLs. (C) Gene expression changes in vehicle control female ACLs. (D) Gene expression changes in female estrogen‐treated ACLs (n = 6 for all male samples, n = 6–7 for female vehicle, n = 6 for female estrogen‐treated). Data represented as box and whisker plot with the whiskers representing the min and max data.

Discussion

4

The objective of this study was to investigate potential mechanisms driving the impaired remodeling response of female ACLs in response to load. Consistent with our hypothesis, our RNA‐sequencing data demonstrated that female ACLs had a muted response to mechanical loading compared to male ACLs. Specifically, female ACLs had significantly fewer DEGs in response to load compared to male ACLs (Figure 1A). This is despite the fact that there were minimal differences in gene expression between freshly harvested female and male ACLs (Figure 1B). This suggests that the differences in the number of DEGs in response to load between male and female ACLs are due to differences in mechanobiology rather than basal biological differences. Furthermore, multiple mechanotransduction pathways were significantly enriched in male ACLs in response to load, while none were enriched in female ACLs (Table 4). Surprisingly, half of these mechanotransduction pathways were predicted to be downregulated in the male samples. However, a number of prominent mechanotransduction pathways, including MAPK (i.e., ERK) and PI3K/AKT, were enriched without likely being downregulated, suggesting that these pathways may be important for the mechanobiological response of male ACLs. Additionally, when comparing the response to load between male and female ACLs, FAK (focal adhesion kinase) signaling was predicted to be further inhibited with load in female ACLs compared to male ACLs (Tables 4 and 5). Therefore, FAK signaling may also help explain sex differences in ACL mechanobiology, given that FAK plays a prominent role in cell‐matrix mechanotransduction and regulates the mechanoresponse of tendon cells [53]. Together, our RNA‐sequencing data support the hypothesis that female ACLs have a muted mechanobiological response to load, potentially inhibiting their ability to repair tissue damage that predisposes the ACL to rupture.

Our secondary hypothesis was that the impaired mechanobiological response in female ACLs is due to estrogen signaling. Surprisingly, there were no basal differences in ER target genes between male and female ACLs (Table 3). This could be due to the fact that we did not measure estrogen concentrations to confirm that estrogen levels were higher in female rabbits at the time of euthanasia. While older rabbits closer to sexual maturity (16–24 weeks according to the supplier) may have produced a greater difference in gene expression, the age of the rabbits used in this study coincides with the onset of sex differences in serum estrogen levels [45]. Additionally, since rabbits are induced ovulators, female rabbits exhibit relatively stable serum estrogen concentrations over time [54]. These data suggest that the estrogen concentration should be consistently higher in the female ACL samples, suggesting that the chosen estrogen targets are not representative of estrogen signaling in rabbits. Downstream targets of estrogen activation are tissue and species‐specific [55, 56], and there are no data available on ER target genes in the ACL or rabbit, providing difficulties in validating ER activation. Additionally, we did not observe differential expression of the chosen ER targets in the female ACLs in response to load. However, the ESR‐mediated signaling pathway was enriched in the female ACL response to load (Supporting Information S1: Table 9). Additionally, the extra‐nuclear estrogen signaling was enriched and trending toward inhibition in the male ACL response to load (Supporting Information S1: Table 8). Furthermore, upstream regulator analysis revealed inhibition of ESR1 in the male ACL response to load (Table 8). Therefore, while rabbit‐specific markers of estrogen signaling may need to be identified, our RNA‐sequencing data suggest that estrogen‐related signaling may help explain sex differences in ACL mechanobiology.

To further investigate this, we cyclically loaded male and female ACLs in the presence of exogenous estradiol and measured gene expression of markers of ECM remodeling. Unfortunately, we were not able to replicate our previous work that utilized the identical loading protocol and identified a sex dependent response in these genes with loading [14]. Given that there was no effect of cyclic loading on our chosen remodeling genes, we were unable to evaluate the effect of estrogen (Figure 3, Supporting Information S2: Figure 3A,B). Power calculations were conducted to determine if this study is underpowered; however, the number of samples required to reach significance was infeasible (Supporting Information S1: Table 12). It is possible that the lack of effect with cyclic loading is due to a poor choice of gene markers for tissue remodeling. Specifically, while our RNA‐sequencing data found 259 DEGs with loading in male ACLs, only two of those genes were chosen as markers of tissue remodeling for our PCR experiments (MMP2, IL1β) (Table 2). However, when looking at the pathway analysis related to extracellular organization, there were only three pathways enriched when comparing the response to load between female and male ACLs (Table 6). Furthermore, only a few genes enriched in these pathways are relevant to remodeling of the bulk ECM (MMP2, MMP14, COL5A3), whereas the others are more relevant to the local pericellular matrix (ADAM10, ADAM12, COL4A1, COL4A2, LAMB1) (Supporting Information S1: Table 11). Together, this suggests that the mechanobiological difference between male and female ACLs may not involve bulk ECM remodeling. Still, the RNA‐sequencing data suggest that ACLs do have a sex‐dependent mechanobiological response, and future work is needed to identify the role of estrogen in driving this differential response.

It is also possible that the absence of an effect of estrogen on ACL explant gene expression was due to insufficient estradiol treatment. However, we used the same concentration (10^−8^ M) that increased expression of COL1 and COL3 in a previous study of isolated ACL fibroblasts [18]. Additionally, our chosen concentration is three orders of magnitude greater than serum estrogen concentrations found in rabbits (3.5 × 10^−11^ M) [25, 26]. It is possible that the ERs are already bound with endogenous estradiol in the female rabbits and that the male rabbits may have had a greater estradiol uptake. However, since there was no change in gene expression with treatment in the male explants, this suggests that there was no effect of estrogen on these genes, even if there was greater estradiol uptake in the male explants. Another possibility is that there was sufficient estradiol in the culture media to saturate estrogen signaling. While we did not use charcoal‐stripped FBS in our culture media, there is only approximately 2 pg/mL (< 10^−11^ M) of estradiol in 10% FBS, which was not sufficient to saturate the response of isolated ACL fibroblasts to estradiol [18]. A final possibility is that the effect of estrogen is no longer present at the end of our culture period (28 h), given that the effect of estrogen treatment on collagen gene expression decreased after 24 h [18].

Two final limitations to this study were that we did not evaluate cell viability in our explants and that our RNA‐sequencing sample size was small for the loaded female samples due to the removal of an outlier (n = 2) (Supporting Information S2: Figure 1). While prior experiments demonstrated that our protocols maintain cell viability in rabbit ACL explants [14], it is possible that the cell death may have affected our current results. Still, our RINs for the RNA‐sequencing data were high (7.5 ± 0.3) and our housekeeping genes were comparable between groups (Supporting Information S2: Figure 5), suggesting that there was comparable RNA quality and quantity between samples. Second, while it is possible that having a smaller sample size could reduce statistical power and lead to fewer DEGs, this is likely not the reason for the low number of DEGs in the female ACLs' response to load. Specifically, it is clear from the PCA plot of the samples included in the sequencing analysis (Supporting Information S2: Figure 1B) that the male‐loaded samples separated from the male static samples to a greater amount than the female samples, which is consistent with the fewer number of DEGs in the female ACLs in response to load.

In conclusion, our findings suggest that female ACLs have a muted mechanobiological response to load compared to male ACLs. This may be due to decreased activation of mechanotransduction pathways, including FAK, MAPK (i.e., ERK), and PI3K/AKT, as well as increased estrogen signaling in female ACLs in response to load. These data support the novel hypothesis that females exhibit an impaired mechanobiological response to load, making female athletes more susceptible to ACL tears. Future work could determine if these changes in gene expression lead to changes in tissue structure and mechanics by comparing the structural and mechanical changes that occur between male and female ACLs in response to load. This will enable a direct comparison of the male and female ACLs' ability to repair tissue damage. Future work is needed to understand better the role of estrogen signaling in ACL mechanobiology. Together, these data can provide insight into the ability of female ACLs to repair tissue damage and possible approaches for reducing the risk of ACL ruptures in females by modulating estrogen levels (e.g., contraceptive use).

Author Contributions

Lauren Paschall carried out experiments and developed the protocols needed. Maxwell Konnaris provided input on the bioinformatic analysis. Lauren Paschall wrote the manuscript with support from all other authors. Erdem Tabdanov, Aman Dhawan, and Spencer E. Szczesny contributed to the interpretation of the results. Lauren Paschall and Spencer E. Szczesny conceived the original idea. All authors approve of the final submitted manuscript.

Supporting information

Supplementary Figure 1: Principal component analysis. Supplementary Figure 2: Comparison of vehicle control female and male ACL gene expression at baseline. Supplementary Figure 3: Comparisons of the male and female ACLs remodeling response to load. Supplementary Figure 4: Comparisons of the male and female ACLs estrogen receptor target genes in response to load. Supplementary Figure 5: Comparison of housekeeping gene expression between groups.

Supplementary Table 1: Summary of sample sizes for each experimental condition. Supplementary Table 2: List of Taqman probes and their efficiencies. Supplementary Table 3: List of differentially expressed genes in male ACL response to load. Supplementary Table 4: List of differentially expressed genes in female ACL response to load. Supplementary Table 5: List of differentially expressed genes in sex effect load (female ACL response to load compared to male ACL response to load). Supplementary Table 6: List of differentially expressed genes in sex effect static (statically loaded female ACLs relative to male statically loaded ACLs). Supplementary Table 7: List of differentially expressed genes in sex effect fresh (freshly harvested female ACLs relative to male freshly harvested ACLs). Supplementary Table 8: List of significantly enriched pathways in male ACL response to load. Pathways with a z‐score > 0 are predicted to be activated and z‐score < 0 are predicted to be inhibited. A prediction was unable to be made for pathways with a NA z‐Score. Supplementary Table 9: List of significantly enriched pathways in female ACL response to load. A prediction was unable to be made for pathways with a NA z‐Score. Supplementary Table 10: List of significantly enriched pathways in sex effect load (female ACL response to load compared to male ACL response to load). Pathways with a z‐score > 0 are predicted to be activated and z‐score < 0 are predicted to be inhibited. A prediction was unable to be made for pathways with a NA z‐Score. Supplementary Table 11: List of DEGs present in pathways associated with ECM organization in male ACL response to load and sex effect load. Supplementary Table 12: Sample sizes required based on post‐hoc power calculations conducted on the PCR remodeling genes assuming a power > 0.8 and type I error of 0.05.