Integrative analysis of GTEx data reveals associations between multi-tissue transcriptomes and human spermatogenic dysfunction

2.1 Data collection from GTEx

In this study, all analyses were conducted using publicly accessible data from GTEx. This dataset includes RNA sequencing data from various tissues, as well as metadata for each donor subject (Table 1). The GTEx Analysis V10 dataset was used in this study. Gene expression data and metadata were converted into AnnData format (version 0.11.3) using Jupyter Notebook [10].

2.2 Testis tissue clustering

To perform cluster analysis of the GTEx testis samples, we used bulk RNA-seq data of the testis. Testis clustering was conducted using the Scanpy (version 1.10.4), a tool can be used for clustering analysis of bulk omics data [11]. The gene expression values for each testis, reported as TPM (transcripts per million), were normalized using the sc.pp.normalize_total() function and transformed to log(TPM + 1) using the sc.pp.log1p() method. The ranked genes for each testis cluster were identified using the Wilcoxon method. Additional methods for testis clustering were executed with the default parameters in Scanpy.

2.3 Differential gene expression screening and gene set enrichment analysis

The differentially expressed genes were screened using Diffxpy (v0.7.4). Gene set enrichment analysis in this study was conducted using ShinyGO (version 0.82). The genes exhibiting the most significant differential expression, particularly those involved in critical biological pathways such as cellular metabolism, neuroendocrine regulation, and inflammatory responses, were selected for analysis and presentation.

2.4 Cell type enrichment analysis

To perform cell type enrichment analysis, we first obtained gene expression patterns of different cells in the testis from a testicular scRNA-seq dataset (GSM8037606) [12] using Scanpy. Afterwards, we generated a reference for the cell type enrichment analysis tool xCell [13] using this scRNA-seq data. Then, we performed cell type enrichment analysis on the GTEx testis RNA-seq data to estimate the potential proportions of different testicular cell types within the testis tissue.

2.5 Statistics analysis

The de.test.t_test() method in Diffxpy, which performs Welch’s t-test (a parametric test to compare the means of two independent groups, without requiring equal variances) for genes expressed in two groups, was utilized for screening differentially expressed genes in testicular data. When the testicular data were analyzed alongside other tissues, as the donor subject number reduced, the de.test.two_sample() method in Diffxpy with the parameter test = “rank” was used, which use Wilcoxon rank sum test (a non-parametric test comparing two independent groups by ranking data to detect distributional differences) to identify the differentially expressed genes between the two groups.

3.1 GTEx data collection and testis clustering

The GTEx V10 dataset contains bulk RNA-seq data from various tissues of human subjects. We first extracted the testicular data from GTEx and classify the testicular tissues into two groups, those with good spermatogenic capacity and those with impaired capacity, by comparing the expression of genes related to spermatogenesis and somatic cell proportions within the testis. Then, based on the classification of testicular function (good vs. poor), we divided the subjects’ extragonadal tissues into two corresponding groups and compare the differences between them. This approach aimed to identify which extragonadal tissues show gene expression patterns that may be associated with testicular function (Fig. 1).

Fig. 1.The main strategy of this study. The quality of spermatogenic function was first evaluated based on bulk RNA-seq data from the testis. Based on this assessment, the human subjects were divided into two groups (Testis.good and Testis.poor). RNA differences in non-testicular tissues were then compared between these two groups to identify which extra-testicular tissues show gene expression patterns associated with testicular function. GTEx, Genotype-Tissue Expression database.

As a result, 414 testes from GTEx were collected and clustered into six distinct groups, labeled T1–T6. Within these testicular clusters, T5 and T6 exhibited high expression levels of genes such as the transcription regulator CRTC1 and the metabolic regulator ACOT9. In contrast, clusters T1, T2, T3, and T4 showed elevated expression of genes like PCYT2 and S1PR2. Notably, the SLC9C2 gene was highly expressed in T3, while the TDRD15 gene was highly expressed in T4 (Fig. 2A). We then analyzed the gene markers of somatic cells (VIM marked all somatic cells, CD14 marked macrophage cells, VWF marked endothelial cells, DLK1 marked Leydig cells, and ACTA2 marked myoid cells) and differentiating spermatids (maker genes include HOOK1, TNP2, PRM3, SPATA18, SPATA7 and SPEM3) [14] in these testis clusters, and found that somatic cell marker genes were highly expressed in T5 and T6, whereas maturing spermatid marker genes were highly expressed in T1 and T2 (Fig. 2B). High expression of maturating spermatid genes in T1 and T2 indicates that the cell ratio of spermatids in testes is high. In contrast, in T5 and T6 testis clusters, high expression of somatic cell markers indicates a high ratio of somatic cells in these testes.

For testicular transcriptome analysis, clusters T1 and T2 were grouped as “Good” (superior spermatogenesis) and clusters T5 and T6 as “Poor” (diminished spermatogenesis) (Fig. 2C). We then screened for differentially expressed genes between the Good and Poor testes using Diffxpy. As a result, in the Poor group, 1846 genes were highly expressed and 2190 genes were lowly expressed (Supplementary Table 1). Genes upregulated in Poor testes were primarily linked to chemical stimulus response and cell communication; downregulated genes were linked to spermatogenesis and sperm motility (Fig. 2C). To further analyze the differences between the two groups of testes, we performed cell type enrichment analysis on the RNA-seq data of these testes. The results showed that spermatogonia were more enriched in the Poor group, while late primary spermatocytes/round spermatids and elongating spermatids were more enriched in the Good group. Early primary spermatocytes showed no difference between the two groups (Fig. 2D). For some somatic cells, such as Leydig cells, Sertoli cells, testicular peritubular cells, smooth muscle cells, epithelial cells, and macrophages, these were all more enriched in the Poor group of testes (Fig. 2D). The results of the cell-type enrichment analysis suggest that in the Good group of testes, there are more germ cells at the maturating stage, and the level of inflammation may be lower.

Notably, among these differentially expressed genes, we identified genes associated with cell communication and chemical response, such as GJA1, HIF1A, and INSR, which were highly expressed in Poor testes. Conversely, genes associated with sperm maturation, such as SPAG6, SPEM2, and ZPBP, were expressed at lower levels in Poor testes (Fig. 2E). These differentially expressed genes may not be attributed to the upregulation or downregulation in specific cell types but rather to changes in the cell counts of particular cell populations.

Fig. 2.Classification of spermatogenic capacity based on testicular transcriptomes. (A) Clustering of testes according to their gene expression profiles. (B) The expression patterns of somatic cell marker genes and sperm marker genes in the testis clusters. (C) Gene set enrichment analysis of differentially expressed genes in testes exhibiting reduced spermatogenic ability. The testes were categorized into “good” and “poor” groups based on their spermatogenic capacity. (D) Cell type enrichment analysis of the two groups of testes. (E) Genes that are either upregulated or downregulated in testes with diminished spermatogenic ability. GJA1, gap junction alpha-1; HIF1A, hypoxia-inducible factor 1-alpha; INSR, insulin receptor; SPAG6, sperm-associated antigen 6; SPEM2, SPEM family member 2; ZPBP, zona pellucida binding protein.

3.2 Hypothalamus-pituitary-testis transcriptional relationships

To analyze the changes in gene expression in the hypothalamus and pituitary gland between the Good Testis (T.good) group and the Poor Testis (T.poor) group, we extracted RNA-seq data from human subjects whose testis, hypothalamus, and pituitary gland were all sequenced. A total of 83 subjects were selected, including 45 in the T.good group and 38 in the T.poor group (Supplementary Table 2). Using the differentially expressed (DE) gene screening criteria of p-value < 0.01 and mean log2(TPM + 1) >0.01, we identified 20,567 DE genes in the testes, 324 DE genes in the pituitary, and 74 DE genes in the hypothalamus (Fig. 3A and Supplementary Table 3). Pituitary DE genes were enriched in copper response and glycolysis pathways (Fig. 3B), while hypothalamic DE genes were enriched in cadmium response and thyroid hormone transport (Fig. 3C). Notably, the pituitary showed ~4.4 times more DE genes than the hypothalamus, suggesting a stronger transcriptional link to testicular function.

Within the HPT axis, hypothalamic GnRH and pituitary FSH/LH play critical endocrine roles in regulating testicular spermatogenesis. Interestingly, our analysis revealed no significant expression change in the core neuroendocrine genes GNRH1 (hypothalamus) or FSHB/LHB/CGA (pituitary) between the T.good and T.poor groups (Fig. 3D), suggesting that spermatogenic impairment may involve alternative regulatory mechanisms. The homeobox protein gene HMX3 is essential for the development of the human inner ear and the hypothalamic-pituitary axis [15]. The expression of glucokinase in the pituitary may function as a glucose sensor to regulate pituitary secretion [16], while the glucokinase regulatory protein encoded by GCKR can inhibit the activity of glucokinase [17]. In the T.poor group, the HMX3 gene was found to be downregulated in the pituitary, whereas the GCKR gene was upregulated (Fig. 3E). Neuromedin-S peptides are believed to be involved in the regulation of the hypothalamic-pituitary-adrenal axis [18]. The vasoconstrictor endothelin-1, expressed in the hypothalamus, can stimulate the secretion of GnRH [19]. In subjects from the T.poor group, the neuromedin-S coding gene NMS was downregulated, while the endothelin-1 coding gene EDN1 was upregulated in the hypothalamus (Fig. 3E). In the brain, the angiopoietin-like 4 protein, encoded by ANGPTL4, primarily functions in metabolic crosstalk between glia cells and neurons [20]. The voltage-dependent calcium channel gamma coding genes CACNG1–8 are associated with schizophrenia [21]. In the T.poor group, the ANGPTL4 gene was found to be upregulated in both the hypothalamus (p < 0.01) and pituitary (p < 0.05), whereas CACNG6 was downregulated (Fig. 3F). All these data suggest that the metabolism and secretion of the hypothalamus and pituitary may be associated with testicular function in T.poor subjects.

Fig. 3.Testis-associated genes in the pituitary and hypothalamus. (A) Differentially expressed (DE) genes in the pituitary and hypothalamus of subjects with varying testicular spermatogenesis abilities. (B) Pathways enriched by DE genes in the pituitary. (C) Pathways enriched by DE genes in the hypothalamus. (D) No significant changes in the coding genes for GnRH, LH, and FSH in the hypothalamus and pituitary of subjects with poorer testicular function. (E) Genes that are differentially expressed in either the pituitary or hypothalamus. (F) Genes that are differentially expressed in both the pituitary and hypothalamus. GNRH1, progonadoliberin-1; FSHB, follitropin subunit beta; LHB, lutropin subunit beta; CGA, glycoprotein hormones alpha chain; ANGPTL4, angiopoietin-related protein 4; CACNG6, voltage-dependent calcium channel gamma-6 subunit; HMX3, homeobox protein HMX3; GCKR, glucokinase regulatory protein gene; NMS, neuromedin-S; EDN1, endothelin-1.

3.3 Transcriptional associations between the adrenal and thyroid glands and testicular function

Beyond the HPT axis, we investigated two other major endocrine regulators, the adrenal and thyroid glands, for their potential associations with testicular function. From GTEx, we identified 85 subjects with RNA-seq data for testes, adrenal, and thyroid glands: 36 with high spermatogenic capacity (T.good) and 59 with impaired spermatogenesis (T.poor) (Supplementary Table 2). Using stringent criteria (q < 0.05; mean log2(TPM + 1) >0.01), we observed tissue-specific patterns: testes showed extensive alterations (20,972 DE genes), adrenal glands moderate changes (260 DE genes), and thyroid glands minimal variation (21 DE genes) (Fig. 4A and Supplementary Table 4). The magnitude of transcriptional changes revealed a striking difference in endocrine-testis associations, with adrenal glands exhibiting 12-fold more DE genes than thyroid glands. This indicates that the adrenal gland has a closer link with testicular status than that of the thyroid gland.

In the T.poor group, we found the immune response-related genes, such as the attractin coding gene ATRN and the C3/C4 complement encoding genes C3 and C4B, were upregulated in the adrenal gland. Additionally, only C4B was upregulated in the thyroid gland (Fig. 4B). This suggests that increased inflammatory activity in both adrenal and thyroid tissues may be linked to impaired spermatogenic function in the testes. Among the DE genes, we identified several genes exhibiting consistent expression patterns between the adrenal and thyroid glands. For instance, the RAMP2 gene, which encodes a protein that regulates the activity of receptor for adrenomedullin, calcitonin, and amylin [22, 23, 24], along with the membrane protein coding gene KANTR, were downregulated in both the adrenal and thyroid glands, but were upregulated in the testis of the T.poor group (Fig. 4C). Another gene, NPFF, which encodes a morphine-modulating peptide, was found to be downregulated in both the adrenal and thyroid glands of the T.poor group, but showed no significant change in the testes (Fig. 4D). These data indicate that the inflammatory state and endocrine regulation in both the adrenal and thyroid glands may be associated with the spermatogenic capacity of the testes.

Fig. 4.Differentially expressed genes in the adrenal and thyroid glands associated with testicular function. (A) Differentially expressed (DE) genes in the adrenal and thyroid glands of subjects with varying testicular spermatogenic abilities. (B) Immune response genes were upregulated in the testis, adrenal gland, and thyroid of subjects in the T.poor group. (C) Genes that were downregulated in both the adrenal gland and thyroid but upregulated in the testis. (D) The gene NPFF was downregulated in both the adrenal gland and thyroid, with no significant change observed in the testis. ATRN, attractin; C3, complement C3; C4B, complement C4-B; RAMP2, receptor activity-modifying protein 2; KANTR, KANTR integral membrane protein gene; NPFF, pro-FMRFamide-related neuropeptide FF.

3.4 Multi-tissue transcriptomic associations with testicular function

To investigate systemic associations with spermatogenic capacity, we analyzed bulk RNA-seq data from various somatic tissues in addition to testes, including the lung, liver, stomach, skin (sun-exposed), visceral adipose tissue, left ventricle of the heart, whole blood, tibial artery, and tibial nerve. To control for potential confounding effects of mechanical ventilation on pulmonary transcriptomes, we excluded subjects with a history of ventilator use. This stringent filtering resulted in 11 eligible subjects (6 T.good subjects and 5 T.poor subjects; Supplementary Table 2). Differential expression analysis (p-value < 0.01; mean log2(TPM + 1) >0.01) revealed 6595 DE genes in the testes, 189 DE genes in the lung, 90 DE genes in the liver and fewer than 90 DE genes in each of the other tissues (Fig. 5A and Supplementary Table 5). These data suggest that the lung and liver may be more closely related to testicular function compared to the other tissues or organs.

In the lungs of the T.poor group, the PNMT gene, which encodes a methyltransferase that converts norepinephrine into epinephrine, and acts on phenylethanolamine and octopamine, was downregulated only in the lungs. Conversely, the CCR1 gene, which regulates immune responses, was upregulated in both the lungs and testis (Fig. 5B). In the liver of the T.poor group, the SLC2A8 gene, which encodes a membrane protein responsible for the transport of glucose and fructose, was downregulated in both the liver and testis. Additionally, the GABRE gene, which encodes the gamma-aminobutyric acid receptor subunit epsilon, was upregulated in the liver and testis (Fig. 5C). These data suggest that immune responses, metabolism, and endocrine regulation in the lung or liver may be associated with the spermatogenic capacity of the testis.

Fig. 5.Genes in lung and other non-testicular tissues associated with testicular function. (A) Differentially expressed (DE) genes in various tissues/organs based on differing testicular spermatogenesis abilities. (B) Genes that are differentially expressed in the lung. (C) Genes that are differentially expressed in the liver. PNMT, phenylethanolamine N-methyltransferase; CCR1, C-C chemokine receptor type 1; SLC2A8, solute carrier family 2, facilitated glucose transporter member 8; GABRE, gamma-aminobutyric acid receptor subunit epsilon.

3.5 Transcriptomic associations between testicular function and multiple brain regions

Beyond the hypothalamic-pituitary axis, our study reveals significant correlations between testicular function and various brain regions involved in emotional and cognitive processing. We analyzed bulk RNA-seq data from testes and 11 distinct brain regions in 14 human subjects (6 T.good, 8 T.poor; Supplementary Table 2), identifying tissue-specific patterns of differential gene expression (p-value < 0.01; mean log2(TPM + 1) >0.01): 5965 DE genes in the testes, 163 DE genes in the amygdala, 139 DE genes in the substantia nigra, 133 DE genes in the anterior cingulate cortex, 124 DE genes in the frontal cortex, and fewer than 100 DE genes in each of the other brain regions (Supplementary Table 6 and Fig. 6A).

In these brain regions, the amygdala is involved in emotional processes such as fear. The substantia nigra is associated with voluntary movement and reward processing in humans, and the loss of dopaminergic neurons in the substantia nigra may lead to Parkinson’s disease [25]. The anterior cingulate cortex and frontal cortex also play a role in regulating emotional functions [26]. In the amygdala, we found that the IRX3 gene was upregulated in the T.poor group, while the ALDH1A1 gene was downregulated (Fig. 6B). The 5-hydroxytryptamine receptor gene HTR1E was downregulated in the substantia nigra of the T.poor group, and the neuronal acetylcholine receptor gene CHRNA2 was downregulated in the anterior cingulate cortex and amygdala of the T.poor group (Fig. 6C). In both the anterior cingulate cortex and frontal cortex, we found that the gamma-aminobutyric acid receptor gene GABRG1 and the peroxisomal acyl-coenzyme A oxidase gene ACOX2 were upregulated in the T.poor group (Fig. 6D). These data suggest that testicular function may be associated with brain regions involved in emotional processing.

Fig. 6.Brain regions associated with testicular spermatogenesis. (A) Differentially expressed (DE) genes in brain regions based on varying abilities of testicular spermatogenesis. (B) Genes that are differentially expressed in the amygdala. (C) Genes that are differentially expressed in the substantia nigra. (D) Genes that are differentially expressed in the anterior cingulate cortex and frontal cortex. IRX3, iroquois-class homeodomain protein IRX-3; ALDH1A1, aldehyde dehydrogenase 1A1; HTR1E, 5-hydroxytryptamine receptor 1E; CHRNA2, neuronal acetylcholine receptor subunit alpha-2; GABRG1, gamma-aminobutyric acid receptor subunit gamma-1; ACOX2, peroxisomal acyl-coenzyme A oxidase 2.

In the brain, the activation, inhibition, and communication between neurons rely on neurotransmitters that function at synapses. We analyzed additional neurotransmitter receptors that may be associated with testicular function (Supplementary Table 7). As a result, we found that the dopamine receptor gene DRD4 was downregulated in the hippocampus of the T.poor group (Fig. 7A), while the glutamate receptor gene GRIN2D was downregulated in the amygdala of the T.poor group (Fig. 7B). In the substantia nigra, in addition to HTR1E, the acetylcholine receptor gene CHRM4, the glycine receptor gene GLRB, the gamma-aminobutyric acid receptor genes GABRA2 and GABRB3, as well as the glutamate receptor genes GRIK2 and GRM1 were all downregulated in samples of T.poor group (Fig. 7C). These data indicate a significant correlation between testicular function and neuronal communication in specific brain regions, particularly in the substantia nigra.

Fig. 7.Differentially expressed brain neurotransmitter receptors associated with testis function. (A) The dopamine receptor gene DRD4 was downregulated in the hippocampus of subjects in the T.poor group. (B) The glutamate receptor gene GRIN2D was downregulated in the amygdala of T.poor subjects. (C) Several neurotransmitter receptor genes were downregulated in the substantia nigra of T.poor subjects. DRD4, D(4) dopamine receptor; GRIN2D, glutamate receptor ionotropic; CHRM4, muscarinic acetylcholine receptor M4; GLRB, glycine receptor subunit beta; GABRA2, gamma-aminobutyric acid receptor subunit alpha-2; GRIK2, glutamate receptor ionotropic, kainate 2; GABRB3, gamma-aminobutyric acid receptor subunit beta-3; GRM1, metabotropic glutamate receptor 1.