Males and Females Have Distinct Molecular Events in the Articular Cartilage during Knee Osteoarthritis

Osteoarthritis (OA) is a major public health challenge that imposes a remarkable burden on the affected individuals and the healthcare system. Based on the clinical observation, males and females have different prevalence rates and severity levels of OA. Thus, sex-based differences may play essential roles in OA’s prognosis and treatment outcomes. To date, the comprehensive understanding of the relationship between sex and OA is still largely lacking. In the current study, we analyzed a published transcriptome dataset of knee articular cartilage (GSE114007) from 18 healthy (five females, 13 males) and 20 OA (11 females, nine males) donors to provide a slight insight into this important but complex issue. First, comparing female healthy cartilage samples with those of males revealed 36 differential expression genes (DEGs), indicating the fundamental sex-related differences at the molecular level. Meanwhile, 923 DEGs were distinguished between OA and healthy female cartilage, which can be enriched to 15 Reactome pathways. On the other hand, when comparing OA and healthy male cartilage, there are only 419 DEGs were identified, and only six pathways were enriched against the Reactome database. The different signaling response to OA in the male and female cartilage was further enforced by recognizing 50 genes with significantly different OA-responsive expression fold changes in males and females. Particularly, 14 Reactome pathways, such as “Extracellular matrix organization”, “Collagen biosynthesis and modifying enzymes”, “Dissolution of fibrin clot”, and “Platelet Aggregation (Plug formation)”, can be noted from these 50 sex-dependent OA-responsive genes. Overall, the current study explores the Sex as a Biological Variable (SABV) at the transcriptomic level in the knee articular cartilage in both healthy status and OA event, which could help predict the differential OA prognosis and treatment outcome of males and female patients.


Introduction
As the most common form of arthritis, osteoarthritis (OA) is a series of pathology that causes persistent pain, swelling, and reduced motion in the affected joints. For years, OA was identified as an age-related pathology; thus, it has been called "wear and tear" arthritis. During the past few years, OA is increasingly recognized as a highly heterogeneous group of diseases characterized by variable clinical phenotypes, which may contribute to the inconsistency of clinical prognosis and treatment response [1].
With the growing recognition of Sex as a Biological Variable (SABV) in the pathophysiology of a diversity of diseases [2], the impact the sex on OA has also attracted more and more attention. To date, it is well known that OA has a higher prevalence in women than men, as 62% of OA patients are women [3]. Indeed, women have a consistently higher OA prevalence rate than men in all age groups between the 30s to 95 plus [4]. Worldwide motif 4 (ADAMTS4), and ADAMTS5, are not significantly different between the healthy male and female cartilage (Supplemental Table S1).
On the other hand, we identify 10 DEGs with a p-value less than 0.05 that are highly expressed in healthy female cartilage than their male counterparts, and 26 DEGs with a p-value less than 0.05 whose expression level is lower in females ( Figure 1). For all these 36 DEGs, only TSIX transcript, XIST antisense RNA (TSIX) has an adjusted p-value less than 0.05 ( Figure 1C and Supplemental Table S1, highlighted in red). Among the latter 26 genes whose expression levels are lower in females, 15 genes are Y-chromosome linked ( Figure 1C), demonstrating the reliability of the current study. Thus, the different expression levels of non-Y-chromosome-linked genes between males and females may present the SABV at a molecular level (Figure 1).
Pathway enrichment was used to uncover the potential functional interaction among these 36 DEGs, while only 14 genes could be recognized by the Reactome knowledgebase. DEGs that were not recognized by the current Reactome database are summarized in Supplemental Table S2. The Reactome recognized genes were clustered into "chromatin organization", "hemostasis", "disease", "metabolism", "transport of small molecules", "metabolism of proteins", and "extracellular matrix organization." Among them, nine identified pathways have a p-value less than 0.05, but none of them qualified as a significant enrichment that should have an FDR smaller than 0.05 (Table 1 and Supplemental Table S3). We then analyzed the cartilage gene expression changes during OA of males and females separately. First, in the female cartilage, there were 923 DEGs in total, among which 382 were downregulated and 541 were upregulated during OA (Figure 2 and Supplemental  Table S4). Among these genes, 30 significantly downregulated DEGs and 45 upregulated ones were identified with an adjusted p-value less than 0.05 (Supplemental Table S4, highlighted in red). Ranking based on the p-values, the top 15 significantly downregulated genes were summarized in Figure 2C, while the top 15 significantly upregulated in Figure 2D. The list of genes that are significantly differentially expressed in healthy male and female cartilage. DEGs with a statistically significant higher level in females have a negative log 2 FC value, while those highly expressed in males have a positive log 2 FC value. The gene with an adjusted p-value less than 0.05 is highlighted in red. The Y-chromosome linked genes are highlighted in blue font. In the Reactome knowledgebase, 424 of the 923 DEGs could not be matched (Supplemental Table S5); thus, the pathways were enriched based on the other 499 DEGs. Overall, there were 68 pathways with a p < 0.05, among which 15 pathways with an FDR less than 0.05 (Table 2 and Supplemental Table S6). Nine of the 15 pathways are related to ECM organization (Table 2). For the other six pathways, "FOXO-mediated transcription of cell cycle genes", "FOXO-mediated transcription", and "RUNX3 regulated immune response and cell migration" belong to the event "gene expression (transcription)", "Response of EIF2AK1 (HRI) to heme deficiency" belongs to the event "cellular responses to external stimuli", "Interleukin-4 and Interleukin-13 signaling" belongs to the event "immune system", and "Gap junction assembly" belongs to the event "vesicle-mediated transport" ( Table 2). In particular, 55 of 499 identified DEGs were enriched in "extracellular matrix organization", which is the most significant event in the female cartilage in response to OA.  Second, we analyzed the male cartilage in the same way. Male samples have much less OA-responsive DEGs compared with female samples. There were 419 DEGs in total, 186 upregulated and 233 downregulated, among which 18 downregulated and four upregulated DEGs have an adjusted p-value less than 0.05 ( Figure 3 and Supplemental Table S7, highlighted in red). In addition, the top 15 significant upregulated and downregulated genes based on p-value in male cartilage during OA were not as same as those in female cartilage. The top 15 significantly downregulated genes in male cartilage were listed in Figure 3C, while the top 15 significantly upregulated genes in Figure 3D.
In the Reactome knowledgebase, 202 of the 419 DEGs could not be matched (Supplemental Table S8). Thus, the pathways enrichment based on the other 217 DEGs dispersed the molecular events including "immune system", "signal transduction", "neuronal system", "hemostasis", "gene expression (transcription)", "metabolism", "DNA replication", "transport of small molecules", "disease", "metabolism of proteins", "cell cycle", "autophagy", "vesicle-mediated transport", "cellular responses to external stimuli", and "extracellular matrix organization". There are 79 pathways that have a p-value less than 0.05, among which six have an FDR less than 0.05 (Table 3 and Supplemental Table S9). Here, "Response of EIF2AK1 (HRI) to heme deficiency" belongs to the event "cellular responses to external stimuli", "ATF4 activates genes in response to endoplasmic reticulum stress" and "PERK regulates gene expression" belong to the event "metabolism of proteins", "NGF-stimulated transcription" and "Nuclear Events (kinase and transcription factor activation)" belong to the event "signal transduction", and "MECP2 regulates neuronal receptors and channels" belongs to the event "gene expression (transcription)". None of these six pathways are categorized in the event of "extracellular matrix organization".

Male and Female Cartilage Have Significant Different Alteration Genes during OA
To confirm the differences between male and female cartilage in response to OA as observed above, we also compared the gene expression fold change in both sexes and identified 63 DEGs with a p < 0.05 (Supplemental Table S10). By referencing the single-sex OA-healthy cartilage comparison results, genes that do not have OA-responsive alteration(s) in either gender were excluded to eliminate the false positive result and lead to the identification of 50 DEGs (Table 4). Note that none of these genes were detected with an adjusted p-value less than 0.05, while 23 of these 50 DEGs could not be recognized by Reactome (Supplemental Table S11). Based on the 27 Reactome-recognized genes, 60 pathways were enriched (p < 0.05; Supplemental Table S12). Among them, 14 pathways have an FDR less than 0.05, which could be clustered in the events of "Extracellular matrix organization" (including "Extracellular matrix organization", "Collagen biosynthesis and modifying enzymes", "Collagen chain trimerization", "Collagen formation", "Assembly of collagen fibrils and other multimeric structures", "Collagen degradation", "ECM proteoglycans", "Integrin cell surface interactions", and "Anchoring fibril formation"), "Hemostasis" (including "Dissolution of Fibrin Clot", "GP1b-IX-V activation signaling", "Platelet Aggregation (Plug Formation)", and "Platelet Adhesion to exposed collagen"), and "Disease" (including "Diseases of glycosylation") ( Table 5). These results further validate male and female cartilage differences at the molecular event level in response to OA.

Discussion
It is broadly accepted that exploring the OA-responsive biomarkers shared by both genders will pave the path for developing the therapeutics that benefit both male and female OA patients [32]. On the other hand, the distinguished clinical appearance between male and female patients warrants the mechanistic investigation at the molecule level. In the current study, the global gene expression profiles of knee joint articular cartilage from male and female donors of a well-accepted dataset [33][34][35][36][37][38][39][40][41][42], GSE114007, were comprehensively compared to gain insight into the understanding of the SABV not only in the healthy status, but also in the response of OA stimulation.
Firstly, the 36 identified male-vs.-female DEGs in healthy cartilage confirmed the hypothesis that the SABV is not limited to the thickness and articular surface areas [21,43] but extended to the static transcriptomic level. In particular, besides the 15 Y-chromosomelinked genes, several genes among the 36 male-vs.-female DEGs in healthy cartilage have been correlated with OA development and progression. For example, as an intensively investigated long non-coding RNA (lncRNA), XIST is highly expressed in OA cartilage tissue and IL-1β-treated chondrocytes [44] and has anti-apoptosis and chondroprotective effects [45]. On the other hand, another lncRNA, MIR4435-2HG, is downregulated in OA [46] and may have inhibition effects on the progression of OA [47]. Regarding the ECM components, a small leucine-rich proteoglycan (SLRP), epiphycan, plays an important role in maintaining joint integrity, and epiphycan-deficient mice spontaneously develop OA with age [48]; Col1A2 is one of the typical markers for fibrocartilage [49] and MXRA5 is highly expressed in the synovial fluid of OA patients [50]. Some other DEGs identified in our current studies have also been associated with OA in previous investigations. For instance, PDLIM1 is downregulated in IL-1β-treated chondrocytes [51], THY1 is highly expressed in OA cartilage and could be induced by IL-1β [52], and EIF1AY has been identified as one of the 9 OA diagnostic biomarkers [53]. In addition, AQP1 promotes caspase-3 activation and thereby contributes to chondrocyte apoptosis [54], and thus the activation of AQP1 induced by OA process can be used to control the tissue degeneration [55].
In addition, IGFBP4 has been identified as the late response gene of parathyroid hormone-related protein (PTHrP) in chondrocytes [56]. It functions as an IGF-1 inhibitor and participates in the inflammatory response [57]. Meanwhile, IGFALS encodes a serum protein that binds IGFs to increase their half-life and vascular distribution [58]. As the male healthy articular cartilage has a lower expression level of IGFALS and higher expression level of IGFBP4 than female cartilage, we infer that IGF-1 signaling is less activated in male cartilage than their female counterpart.
Note that among these 36 DEGs, only TSIX has an adjusted p-value less than 0.05, indicating the significance of TSIX for gender-dependent biological differences in the articular cartilage. However, the detailed function of TSIX in cartilage remains blank. In addition, the limited available sample could lead to only one DEG identified with an adjusted p-value less than 0.05 identified, while more DEGs with a p-value less than 0.05 (36 DEGs) were recognized. Thus, further studies are undoubtedly encouraged to fully understand the SABV in healthy knee articular cartilage at the molecular level, which warrants a worldwide collaboration for more database collection in a diverse of populations.
Interestingly, when we profile OA-responsive transcriptional changes in male and female cartilage separately, the amount of OA-responsive DEGs with an adjusted p-value less than 0.05 in female cartilage is triple that in male cartilage, indicating more intense OA-induced molecular changes in female cartilage than that in male counterparts. This transcriptomic difference could be correlated with the clinical observation that women experience more severe OA symptoms than men [59,60]. Considering the different total amounts of OA-responsive DEGs, it is no surprise to find that the top 15 OA-responsive upregulated and downregulated DEGs are not identical in male and female cartilage. In fact, male and female cartilage do share some top OA-responsive DEGs with an adjusted p-value less than 0.05, such as CISH, ADM, HLPDA, DDIT3, DDIT4, CFI, ST6GALNAC5, SPOCK1, and TNFSF15. Regarding the Reactome-enriched pathways, "response of EIF2AK1 (HRI) to heme deficiency" is the common significant pathway with adjusted p-value less than 0.05 in response to OA stimulation shared by male and female cartilage. These shared genes and pathways could be considered as potential targets for OA diagnosis and treatment, which can benefit both genders.
The OA-responsive molecular events in female cartilage are tightly clustered in the "extracellular matrix organization", which could explain the reason that female patients have more severe OA-related cartilage defects than males [60,61]. Meanwhile, "FOXOmediated transcription", "RUNX3-regulated immune response and cell migration", and "Interleukin-4 and Interleukin-13 signaling", the pathways with FDR less than 0.05, might be additional key pathways to regulate OA in females. In fact, recent studies demonstrate that FOXO transcription factors modulate autophagy and proteoglycan 4 in cartilage, and conditional knockout FOXOs could induce OA-like changes in the mice [62,63]. On the other hand, ECM degradation does not present as the leading OA-responsive event in the male cartilage. Instead, "ATF4 activates genes in response to endoplasmic reticulum stress", "NGF-stimulated transcription", "MECP2 regulates neuronal receptors and channels", "PERK regulates gene expression", and "Nuclear Events (kinase and transcription factor activation)" were enriched from the OA-responsive DEGs in male cartilage with FDR less than 0.05, indicating a distinct molecular response to OA between male and female cartilage. The activation of the PERK-ATF4-CHOP axis is especially known to mediate impaired cartilage function [64]; however, the effects of these male-specific OA-responsive pathways in arthritis are still unknown.
SABV of cartilage in response to OA was further evaluated by comparing the OAresponse DEGs from both genders directly, by which 50 genes with significantly different expression fold changes were identified, but none of the genes has an adjusted p-value less than 0.05. As expected, "Extracellular matrix organization" is the major sex-relative differential event harboring 9 of the 14 enriched pathways. There are also differences in "Hemostasis" and "Diseases of glycosylation" events. Note that several genes clustered in the event "Hemostasis" (including pathways "Dissolution of Fibrin Clot", "GP1b-IX-V activation signaling", "Platelet Aggregation (Plug Formation)", and "Platelet Adhesion to exposed collagen") have also been investigated in OA-related area. For example, SER-PINE1 has been identified as one of the OA-specific genes in human joint fibroblast-like synoviocytes [65]. While SERPINE2, a contributor for both "Hemostasis" and "Diseases of glycosylation" events, upregulated by IL-1α stimulation in human chondrocytes, and recombinant SERPINE2 may protect chondrocytes by inhibiting MMP-13 expression [66]. Besides, high platelet counts within the normal range are significantly associated with knee and hip OA in women aged above 50 [67].
Considering aging may be an important indicator of OA, it is not a surprise that the donors of the OA groups are older than the healthy group when the dataset was built [33]. Interestingly, specifically grouping the samples in the same dataset GSE114007 by donor age, Chen et al. concluded that age is not a dependent variable for differentially expressed gene identification [41]. Here, as demonstrated in Table 6, no difference regarding donor age between males and females was found in healthy cartilages nor OA samples. Thus, the age contribution on OA-responsive differentially expressed genes, if any, has already been considered in parallel for both genders. Note that comparing healthy cartilage of different age stages for each gender would be an interesting and important topic for gaining more insight on the molecular events in senescence, particularly in a gender-dependent manner. Besides, an inter-cohort validation should be conducted in the future to verify the genes and pathways discovered in the current study. Last but not least, it is the first time that multiple genes and pathways mentioned above are associated with chondrogenic differentiation, maintenance, and pathology. The underlying mechanistic and functional details are largely unknown. No doubt, a huge amount of effort should be devoted on a global base to transferring the discovery here to the real world.

Materials and Methods
By using the keywords "osteoarthritis" and "cartilage" in the NCBI GEO DataSets website [68] with the selection of "Homo sapiens" under the column of "Top Organisms" and "Expression profiling by high throughput sequencing" under the column of "study type", 31 series were identified. After reviewing all these datasets to check if they provided the sex information of the donors, one series (GSE114007) containing transcriptome data of human knee cartilage samples was included in the current study [33]. In this dataset, there were samples from 5 healthy female donors (age 27-57, mean 42), 13 healthy male donors (age 18-61, mean 34.5), 11 OA female donors (age 52-82, mean 66.3), and 9 OA male donors (age 51-71, mean 64.9) ( Table 6 and Supplemental Table S13). According to the original study of this dataset, there is no significant difference between healthy and OA samples in other factors, such as the health condition of the donors, tissue sampling location, and body mass index [33,41]. SRA data of all the samples were downloaded from NCBI SRA website [69]. Following comparisons were conducted: (1) male healthy (HM) cartilage with female healthy (HF) cartilage to explore the baseline molecular differences in the articular cartilage between genders, (2) male OA (OM) cartilage with HM cartilage to detect the molecular changes in response to OA in males, (3) female OA (OF) cartilage with HF cartilage to detect the molecular changes in response to OA in females, and (4) OA-responsive DEGs in males (OM-HM) with that in females (OF-HF) to find the genes altered significantly different between genders during OA (OM-HM vs. OF-HF). Data analyses were performed on the Galaxy platform (UseGalaxy.org; [70]) with an established, broadly validated protocol [71][72][73]. Briefly, the FASTQC RNA-seq reads were aligned to the human genome (GRCh38) using HISAT2 aligner (Galaxy Version 2.1.0+galaxy 5) with default parameters [74]. Raw counts of sequencing read for the feature of genes were extracted by featureCounts (Galaxy Version 1.6.4+galaxy1) [75]. Then, the limma package (Galaxy version 3.38.3 + galaxy3) was used to identify DEGs with its voom method [76,77]. Expressed genes were selected as their counts per million (CPM), value not less than 1 in at least two samples across the entire experiment, while lowly expressed genes were removed for the flowing analyses. The parameters were set as 1 for minimum log 2 Fold change and 0.05 for p-value adjustment threshold. As our current investigation is an explorative study, Benjamini-Hochberg correction was employed in the limma-voom analysis for p value adjustment [78], which is highly recommandated by the limma user guide [79]. To provide FDR control, the limma Test significance relative to a fold-change threshold (TREAT) function was applied to select genes that are more likely to be biologically significant [80], accompanied by the Robust Setting to protect against outlier genes [81]. A trimmed mean of M values (TMM) method was used for normalization among RNA samples. Quasi-likelihood F-tests (ANOVA-like analysis) were achieved to identify DEGs [82]. Genes with fold change (FC) more than 2 and p value less than 0.05 were assigned as DEGs. Heatmap diagrams were conducted in R (version 3.6.3) [83] with packages pheatmap (version 1.0.12), while volcano plots were generated by GraphPad Prism (version 8.2.1; GraphPad Software, Inc., San Diego, CA, USA). Pathway enrichment of identified DEGs was performed against the Reactome knowledgebase [84]. The enriched pathways with a false discovery rate (FDR) less than 0.05 were considered significantly meaningful.

Conclusions
In summary, our current study confirmed SABV in the knee cartilage at the transcriptomic level in both healthy and OA statuses. This study, at least partially, explains the clinical observed sex-relative differences of OA outcomes. Due to the lack of knowledge about some of the identified DEGs, further worldwide collaboration is necessary to comprehensively uncover the sex-relative differences of knee articular cartilage health and disease.

Conflicts of Interest:
The authors declare no conflict of interest.