Transcriptomic differences between fibrotic and non-fibrotic testicular tissue reveal possible key players in Klinefelter syndrome-related testicular fibrosis

Klinefelter syndrome (KS; 47,XXY) affects 1–2 in 1000 males. Most men with KS suffer from an early germ cell loss and testicular fibrosis from puberty onwards. Mechanisms responsible for these processes remain unknown. Previous genomics studies on testis tissue from men with KS focused on germ cell loss, while a transcriptomic analysis focused on testicular fibrosis has not yet been performed. This study aimed to identify factors involved in the fibrotic remodelling of KS testes by analysing the transcriptome of fibrotic and non-fibrotic testicular tissue. RNA sequencing was performed to compare the genes expressed in testicular samples with (KS and testis atrophy) and without (Sertoli cell-only syndrome and fertile controls) fibrosis (n = 5, each). Additionally, differentially expressed genes (DEGs) between KS and testis atrophy samples were studied to reveal KS-specific fibrotic genes. DEGs were considered significant when p < 0.01 and log2FC > 2. Next, downstream analyses (GO and KEGG) were performed. Lastly, RNA in situ hybridization was performed to validate the results. The first analysis (fibrotic vs non-fibrotic) resulted in 734 significant DEGs (167 up- and 567 down-regulated). Genes involved in the extracellular structure organization (e.g. VCAM1) were found up-regulated. KEGG analysis showed an up-regulation of genes involved in the TGF-β pathway. The KS vs testis atrophy analysis resulted in 539 significant DEGs (59 up- and 480 down-regulated). Chronic inflammatory response genes were found up-regulated. The overlap of X-linked DEGs from the two analyses revealed three genes: matrix-remodelling associated 5 (MXRA5), doublecortin (DCX) and variable charge X-Linked 3B (VCX3B). RNA in situ hybridization showed an overexpression of VCAM1, MXRA5 and DCX within the fibrotic group compared with the non-fibrotic group. To summarize, this study revealed DEGs between fibrotic and non-fibrotic testis tissue, including VCAM1. In addition, X-linked fibrotic genes were revealed, e.g. MXRA5, DCX and VCX3B. Their potential role in KS-related testicular fibrosis needs further study.

More information about the patients included in this study can be found in Table 1. Immediately after the biopsy was retrieved, samples were snap frozen in liquid nitrogen after which they were stored at − 80 °C for next generation sequencing. In addition, biopsy samples were fixed in alcohol-formalin acetic acid mixture (AFA; Q02022; international medical products, Oudergem, Belgium) for at least one hour and further dehydrated and fixed with TISSUE TEK-VIP (5990; Sakura, Berchem, Belgium) overnight after which the tissue was embedded in paraffin. For validation by RNA in situ hybridization, sections of 5 mm thick were cut at three different depths using a microtome (SM2010R; Leica, Brussels, Belgium).
Histological examination. The fibrotic score of each testicular tissue sample included in this study was determined as previously described 8 . In brief, paraffin-embedded tissue was used to perform a haematoxylinperiodic acid Schiff 's reagent (H/PAS) staining. The histology of the tissue was scored according to the degree of fibrosis. A mean score of five evaluated image fields was given to each tissue sample: (0) normal histology; (1) normal tubules but interstitial fibrosis; (2) degenerated tubules and interstitial fibrosis; (3) hyalinized tubules and severe fibrosis; (4) severe fibrosis, complete loss of tubules. If the mean score was not a round number, it www.nature.com/scientificreports/ was rounded up to the nearest round number. Only samples with a fibrotic score ≤ 2 were included in the nonfibrotic groups while for the fibrotic groups, only samples with a score of 3 or 4 were included for KS patients and samples with a score of 2-4 for the TA patients were included. Supplementary Fig. 1 shows the H/PAS staining of every testicular tissue piece included in this study.
To examinate the histology of the included testicular tissue more thoroughly, the percentage of tubules with ongoing spermatogenesis, tubules with an SCO image and hyalinized tubules was calculated. An immunohistochemical staining against melanoma associated antigen 4 (provided by Giulio Spagnoli, University of Basel, Switzerland) was performed to evaluate the presence of spermatogonial stem cells and spermatogonia as described before 8 . The number of tubules, as well as the area of interstitium, was evaluated in three tissue pieces at three different depts from every patient included, using the Fiji software 25 . The results of this analysis can be found in Table 1.
RNAseq: Wetlab protocol. Total RNA was extracted from the testicular tissue using the RNeasy Micro kit (74004; Qiagen, Antwerp, Belgium) following the manufacturer's guidelines. RNA concentration was determined with a Nanodrop spectrophotometer (Isogen Life Science, De Meern, the Netherlands). The subsequent steps of the RNAseq procedure were carried out at the BrightCore facility (UZ Brussel). Extracted RNA from KS, SCO, TA and FC biopsies was analysed on the Agilent 2100 Bioanalyzer system (G2939BA; Agilent technologies), leading to an RNA integrity number (RIN) score. Samples with RIN scores above 4.7 were considered pure enough to continue with the experiments. An overview of the RIN scores for the testicular tissue included in this study is shown in Aupplementary Table I. Subsequently, an input of 150 ng RNA per sample was processed using the KAPA RNA HyperPrep with RiboErase (08098131702; Roche, Vilvoorde, Belgium) library preparation kit with the following steps: oligo hybridization and RNA depletion, 2.2 × bead-based RNA depletion clean-up, DNase digestion, 2.2 × bead-based DNase digestion clean-up, RNA elution, fragmentation and priming, first strand synthesis, second strand synthesis and A-tailing, adaptor ligation, first 0.63 × bead-based post-ligation clean-up, second 0.7 × bead-based post-ligation clean-up, library amplification and library amplification cleanup. These steps were carried out according to the manufacturer's protocol. Library validation was performed through the measurement of cDNA fragment concentrations, using the Victor Nivo Multimode Microplate Reader (HH35000500; PerkinElmer, Zaventem, Belgium), and the length of the different fragments was measured by the Fragment-analyzer (M5311AA; Agilent technologies). These steps were carried out according to the manufacturer's protocol in order to obtain a library insert size of 300-400 bp.
Libraries were pooled according to the two different patient groups (fibrotic vs non-fibrotic tissue) as well as pooled within the fibrotic group (KS vs TA). The pooled libraries were sequenced with the Illumina NovaSeq 6000 (Illumina, Eindhoven, the Netherlands). The sequencing stats (coverage and percentage of aligned reads) are shown in Supplementary Table I. In addition, the sequencing are publicly available in the GEO database (GSE200680).
RNAseq: differential expression and downstream analysis. The reads were mapped to the human genome (hg19) using the tool STAR v2.5.0c 26 . The mapped reads were then translated into a quantitative measure of gene expression, counting the number of reads that map to each exon/gene using HTSeq v0.11.0 27 .
Following quantification of expression levels, the differential expression between the different groups was determined. Focus was put on the protein coding genes (19,265 genes, downloaded from Ensemble biomart in May 2019). In a first step, a PCA analysis was run to check for outliers ( Supplementary Fig. 2). Next, the differential expression analysis was run after which the genes that were significantly up-or down-regulated with an absolute log2 fold change > 2 after shrinkage (type apeglm) and an adjusted p-value < 0.01 were extracted. This was done with the DESeq2 package v1.24.0 in R v3.6.1 28 . The correction of multiple testing was done with the Benjamini-Hochberg FDR implementation in DESeq2.
Then, two gene ontology (GO) analyses were run using (i) the set of up-regulated genes and (ii) the set of down-regulated genes using the enrichGO function of the Bioconductor package clusterProfiler. On these two gene sets, a kyoto encyclopedia of genes and genomes (KEGG) analysis was run using the enrichKEGG function of the same package 29 . RNA in situ hybridization. RNA in situ hybridization was used as validation method to confirm the overexpression of VCAM1, DCX and MXRA5 in the fibrotic tissue. Two samples of one patient per diagnosis category (KS3, TA4, FC5, SCO1) were included for the validation. This was performed by the 'extended RNA analysis' facility at VUB (https:// lmmo. resea rch. vub. be/ en/ era-facil ity). The RNAscope 2.5 HD Duplex Assay (322430; Advanced Cell Diagnostics, Newark, CA, USA) was used according to the manufacturer's instructions with adapted pretreatment conditions (15 min target retrieval and 15 min protease plus). Probes used were B. subtilis duplex negative control (320751; Advanced Cell Diagnostics), H. sapiens duplex positive control (321641; Advanced Cell Diagnostics) and H. sapiens MXRA5 (419691; Advanced Cell Diagnostics), VCAM1 (440371-C2; Advanced Cell Diagnostics) and DCX (489551; Advanced Cell Diagnostics). Imaging was performed with the Aperio slide scanner (GT 450; Leica). Quantification of the RNA in situ hybridization was performed with the HALO software (Indicalabs) 30 . Two testicular tissue pieces from one patient of each diagnosis group (KS3, SCO1, FC5, TA1) were included for the quantification. A mean percentage of cells expressing the genes of interest (VCAM1, MXRA5 and DCX) was calculated.

Results
Transcriptomic analysis: fibrotic vs non-fibrotic tissue. No significant difference (p = 0.568) was found between the fibrotic and non-fibrotic group concerning age at biopsy (Supplementary Fig. 3). Transcriptomic analysis resulted in a total of 734 significant DEGs, of which 167 were up-regulated and 567 were down-regulated in the fibrotic group (KS + TA) in comparison to the non-fibrotic group (SCO + FC) (Fig. 1A,B; Supplementary Information 1-Significant DEGs). Results showed that all FC samples clustered together, as did the KS samples, while the gene expression in SCO and TA samples was more diverse. Nevertheless, the TA samples with a fibrotic score of 3 clustered together with the KS samples (Fig. 1C). GO analysis showed diverse biological functions for the up-regulated DEGs. The top ten of the most up-regulated functions is shown in Fig. 2A. The most interesting function found through this GO analysis included the extracellular structure organization, since testicular fibrosis can be attributed to an excessive deposition of ECM components. The genes involved the extracellular structure organization, were, among others, decorin, lumican, vascular cell adhesion molecule 1 (VCAM1) and several collagen-related genes. The network of genes involved in this biological function is displayed in Fig. 2B. GO analysis for the down-regulated genes resulted in biological functions which are mainly related to spermatogenesis and fertilization, such as spermatid development and differentiation ( Supplementary Fig. 4). KEGG analysis showed the up-regulation of genes involved in different pathways. Figure 2C shows the 8 most up-regulated pathways and their intersection size, meaning the number of up-regulated genes which are involved in these pathways, separately or overlapping. An up-regulation of genes involved in the transforming growth factor beta-1 (TGF-β1) pathway was found. The up-regulation of TGF-β1 in KS patients compared with SCO and TA patients and fertile controls has been validated by RT-qPCR analysis and was found significant ( Supplementary Information 2 and Supplementary Fig. 5). Moreover, a down-regulation of genes involved in the extracellular matrix (ECM)-receptor interaction pathway was seen. A total of 15 down-regulated www.nature.com/scientificreports/ and 11 up-regulated DEGs were located on the X-chromosome (Fig. 2D). An overview of the number of DEGs per chromosome is shown in Supplementary Fig. 6A.
Transcriptomic analysis-KS vs TA. Next, a transcriptomic analysis to find DEGs between the two fibrotic groups, KS and TA samples, was performed. The heatmap showed a clustering of KS samples and a clustering of TA samples (Fig. 3A). A total of 482 significant DEGs were found, of which 57 were up-regulated and 425 were down-regulated in the KS samples compared to the TA samples (Fig. 3B). Enrichment analysis revealed the most significant biological functions for both the up-and down-regulated genes. Among the up-regulated biological functions, the chronic inflammatory response (PTGES/GJA1/CYP19A1) was found (Fig. 3C). Processes involved in cilium movement, meiotic cell cycle and the cellular process involved in reproduction in multicellular organisms were found as most down-regulated biological functions ( Supplementary Fig. 7). KEGG analysis for the up-regulated genes showed the steroid hormone biosynthesis as most up-regulated pathway. For the down-regulated genes, several pathways were found, including ABC transporters, tight junction and glycosaminoglycan degradation (Supplementary Fig. 8). www.nature.com/scientificreports/ Twenty-one down-regulated genes were located on the X-chromosome, while three up-regulated DEGs could be found on the X-chromosome (Fig. 3D). An overview of the number of DEGs per chromosome for the second analysis is shown in Supplementary Fig. 6B. Concerning the X-linked genes, three genes were found in both the first analysis, comparing fibrotic and non-fibrotic tissue, and the second analysis, comparing fibrotic tissue with and without a normal karyotype: up-regulated matrix-remodelling associated 5 (MXRA5), up-regulated doublecortin (DCX) and down-regulated variable charge X-Linked 3B (VCX3B).
Validation. RNA in situ hybridization was performed to visualize three selected genes from our two analyses, namely VCAM1, MXRA5 and DCX. VCAM1 was selected due to its role in the inflammation pathway, while MXRA5 and DCX were chosen as the two up-regulated X linked genes which overlapped in the two analyses of this study. The duplex array combined VCAM1 with MXRA5 and with DCX. In both combinations, an overexpression of VCAM1 was clearly observed in KS samples compared with FC samples. A total of 84% of cells within the KS samples expressed VCAM1 while this expression was only seen in 17% of all cells in the FC samples. VCAM1 RNA molecules were mainly present in the interstitial tissue in FC and SCO samples, while in KS samples, VCAM1 expression seems to correlate with Leydig cell hyperplasia. VCAM1 expression was detected in 36% and 16% of all cells in the TA and SCO samples, respectively. MXRA5 expression was barely seen in FC (7% of cells) and SCO (2% of cells) samples, while visibly more expression was detected in KS (41%) and TA (14%) samples. Quite some DCX RNA molecules could be distinguished in KS (28%) and TA (16%) samples, in contrast to SCO (4%) and FC (7%) samples. In FC samples, some positive patches were seen in the centre of the www.nature.com/scientificreports/ seminiferous tubules, however, this was also visible in the negative control of all FC samples included, confirming these patches did not correspond with real DCX expression ( Fig. 4 and Supplementary Fig. 9).

Discussion
The exact association between the supplementary X chromosome and the pathophysiology seen in men with KS has not yet been elucidated. In this study, a differential gene expression analysis was performed on KS and control testicular tissue to reveal a number of genes which may be responsible for the KS-related testicular fibrosis. First, a gene expression analysis of fibrotic versus non-fibrotic testis tissue was performed. In aging men, testicular fibrosis occurs naturally 31 , however, a statistical analysis showed no significant differences between the two patient groups in their respective age periods. The most down-regulated biological functions were associated with fertilization and spermatogenesis. For the up-regulated biological functions, both processes involving hormonal changes and processes of extracellular structure organization were found. Testicular fibrosis is characterized by excessive ECM synthesis and it has been shown that the tubular wall is remodelled in men with impaired spermatogenesis 32 . A recent study of our group revealed significant changes in the ECM of KS compared to SCO and FC, i.e. a loss of alpha smooth muscle actin 2 and an altered expression of collagen I and IV 33 . Several genes are involved within the ECM organization, including some collagen genes, VCAM1 and small leucine-rich proteins such as decorin and lumican. Decorin has already been shown to be highly expressed in testis tissue from men with KS, especially in hyalinized tubules 34 . Loss of lumican has been shown to ameliorate cardiac fibrosis 35 . VCAM1 is expressed by the Leydig cells and can also be found in the basal parts of the Sertoli cells and is known to play a role in the testicular immunoregulation 36,37 . Recently, VCAM1 has been shown to be crucial for the reconstruction of testicular interstitial tissue in vitro 38 . In addition, VCAM1 is a key element of the inflammation pathway, contributing to several immunological disorders as well as cancer. Furthermore, www.nature.com/scientificreports/ VCAM1 has already been associated with pulmonary fibrosis 39,40 . Our results also showed that DEGs involved in the TGF-β1 pathway were significantly up-regulated. It is known that testicular fibrosis is mediated by TGF-β1 and that its expression is increased in men with impaired spermatogenesis 41,42 . Moreover, TGF-β1 production is regulated by decorin and vice versa 43 . In addition, it has been shown that VCAM1 is a TGF-β1 inducible gene 39 .
In accordance, our research also revealed a significantly higher expression of TGF-β1 in KS biopsies compared with fertile controls. Recently, a study from Mahyari et al. revealed different up-regulated genes related to fibrosis, specifically expressed by Leydig cells, through a single-cell analysis of KS biopsies 24 . Three of these profibrotic genes (IGFBP5/WFIKKN2/ INHBA) were also found up-regulated in our data. Secondly, a gene expression analysis between the two fibrotic groups, KS and TA, was performed. This analysis, in combination with the first analysis, revealed an overlap of three genes: MXRA5, DCX and VCX3B. MXRA5 is a proteoglycan which has a function in cell adhesion and ECM remodelling. MXRA5 has been associated with several diseases and conditions, including preeclampsia and colorectal cancer 44,45 . A few years ago, MXRA5 was identified as an anti-inflammatory and anti-fibrotic molecule, regulated by TGF-β1 46 . However, no exact function of MXRA5 within the testis has been elucidated. DCX has been elaborately studied within the brain, as it is a microtubule-associated protein required for cerebral cortex development. In addition, it has been correlated with several disorders of the neuronal structures, including lissencephaly, epilepsy and developmental dyslexia 47 . A novel human DCX domain containing gene which is expressed mainly in the adult testis was discovered several years ago, however, the role of this protein within the testis is unknown 48 . VCX3B is a member of the VCX/Y gene family, of which all the members are exclusively expressed in male germ cells. The exact function of these genes is currently unknown, although the deletion of one of the family members, VCX-A, has already been associated with mental retardation, a condition also found in some patients with KS 49,50 . The single cell data set from the KS tissue included in Laurentino et al. showed expression of MXRA5 in PTMCs, expression of DCX in Leydig cells and PTMCs and expression of VCX3B in germ cells 22 . Moreover, the KEGG analysis revealed genes involved in steroidogenesis pathogenesis to be up-regulated in the KS samples, which is in accordance with research showing an increased intratesticular testosterone level in men with KS 51 . Although no changes in the blood vessel density were found in KS testicular biopsies compared to controls, as reported in our previous research 34 , the above finding may contribute to the hypothesis that a disturbed vascularization is present within the KS testis 52,53 . In men with KS, as in females, one of the two X-chromosomes is silenced to compensate for the dosage of X-linked genes. However, it is known that about 15% of these genes escape X-inactivation, resulting in an overexpression of several X-linked genes. This overexpression could be related to the clinical features seen in men with KS, including the testicular fibrosis [54][55][56] . In this study, the only known escapee X-linked DEG found was MXRA5 57 , substantiating its role in KS-related testicular fibrosis. In addition, down-regulation of TEX11 was found in our second analysis. Mutations of this X-linked gene were previously shown to result in meiotic arrest and azoospermia in infertile men 58 . Furthermore, down-regulation of C14ORF39 and SYCE1 were found when comparing KS and TA samples. Variations of these genes have been identified as the reason for infertility in NOA patients 59 .
Winge et al. recently compared all published studies concerning differentially expressed analyses and revealed an overlap of only eight transcripts which were up-regulated in at least two studies on testis tissue, while no down-regulated genes overlapped between the studies 60 . Of these eight genes, only one was found differentially expressed in our analysis of fibrotic versus non-fibrotic tissue, e.g. DLK1. This marker for immature Leydig cells has already been shown to be up-regulated in adult men with testicular pathologies, including males with KS 61 . When also including the transcriptomics studies performed on blood samples from men with KS, a total of 139 genes were found to be up-or down-regulated in at least two studies. Of these genes, five could be found in our first analysis (fibrotic vs non-fibrotic): DLK1, RAB34, PLSCR2, GAPDHS and GPA33. Nevertheless, GPA33 was found to be up-regulated in our results while it was found down-regulated in the studies by D' Aurora et al. and Belling et al. 15,20 . When comparing the overlapping genes of the transcriptomic studies with our second analysis (KS vs TA), only one overlapping gene was found: CNGA2. Nevertheless, this was down-regulated while it was found up-regulated in the study by Winge et al. 21 . The six genes that were found in our analyses and the previously published studies, were also found through the single cell analysis of a KS testis sample by Laurentino et al. 22,60 . Since the other studies compared KS samples with fertile samples, an additional analysis was performed, comparing the differentially expressed genes between the KS and FC samples. In this analysis, 26 of the 139 overlapping genes published in Winge et al. were found 60 . A comparison of the results between the different published papers and our results has been included as Supplementary Information 3-Comparison studies.
In this study, snap-frozen testicular biopsies were used to perform the gene expression analysis, in contrast to most other studies using blood samples or extracted RNA from paraffin-embedded tissue. Nevertheless, this study has a few limitations. Only a small number of samples could be included, as KS samples for research are rather scarce. Since ideal control samples with a fibrotic score of zero were not available, we included samples with fibrotic scores one and two in the non-fibrotic group. A last limitation includes the use of bulk RNA sequencing, which is sensitive to differences in cell-type proportions. Meaning that, each time a differentially expressed gene is reported, it should be taken into account that the differential expression can occur due to cellularity differences between the samples.

Conclusions
The gene-expression in fibrotic and non-fibrotic testis tissue was compared, revealing genes such as VCAM1, which may play a role in testicular fibrosis. In addition, this analysis revealed a pertinent role for genes involved in the TGF-β1 pathway. Secondly, a differential gene expression analysis between the fibrotic groups, KS and TA samples, was performed. An overlap of X-linked genes from our two analyses, revealed the up-regulation of MXRA5 and DCX as well as the down-regulation of VCX3B. The altered expression of these genes may lead www.nature.com/scientificreports/ to mechanisms leading to the KS-related testicular fibrosis. Nevertheless, the function of these genes needs further study.

Data availability
The datasets used and analysed during the current study have been deposited in the GEO public database (GSE200680).