Pathway Analysis of Testicular Dysfunction Microarray Data Implicates Alteration of ITGB1 Expression as a Key Determinant of Spermatogenic Impairment

1. Introduction

Infertility is becoming increasingly recognised as an important medical concern, as infertile phenotypes and their underlying causes are substantially diverse, markedly in men, and remain poorly understood despite available treatments [1]. Despite thorough clinical evaluation, a substantial proportion of infertile men have no identifiable cause and are therefore classified as having primary idiopathic impaired spermatogenesis (IS) [2]. The lack of a known aetiology limits diagnosis and the development of targeted therapies, highlighting the need to better understand the molecular mechanisms underlying impaired spermatogenesis. At the cellular level, impaired spermatogenesis can result from germ-cell death, cell-cycle arrest or both, and can potentially occur at any stage of life [3]. Testis maldescensus, failures of hormonal support, and a wide range of exogenous factors may contribute to this impairment of spermatogenesis [4], likely resulting a significant reduction in numbers of sperm in the ejaculate, leading to oligozoospermia or azoospermia [5]. The latter is a form of infertility that presents particular challenges, since the only available interventions are assisted conception using sperm extracted from the testis or epididymis– if they can be found – or donor sperm [6]. DNA microarray has been used extensively to detect RNA expression patterns, and the findings are usually deposited in publicly accessible databases. Among these, the Gene Expression Omnibus (GEO) database is a free microarray data resource that was created by the National Centre for Biotechnology Information (NCBI) in 1999 [7]. Over the last decades, over 20,000 microarray datasets have been published into the GEO database [8]. For instance, in 2007, Dubé et al., discovered using microarrays that 414 genes were expressed at least 2-fold higher or lower in IS patients compared with fertile men, in whom the genes CRISP1, SPINLW1, FAM12B and DEFB129 were up-regulated. These genes were involved in gene transcription, immune function and fertility but it was determined that genes involved in water and ion transport, including CFTR, AQP5 and KCNK4, were the most likely to be associated with IS because their expression was repressed in infertile patients [9]. In another study, the gene ART3 was identified as a susceptibility gene for IS using genome-wide gene expression analysis with microarray on IS and CT samples of testicular tissue [10].

All the studies mentioned above, however, suffer from the same problem, namely that any pathology that impairs development of a particular cell type or types within an organ, thereby also perturbs the relative abundance of the different cell types in that organ. Consequently, it is very difficult to interpret accurately expression data obtained from such studies and is particularly problematic in the testis, where spermatogenic cell-type specific effects are common. One solution to this could be to use pathway analysis on expression data from testis samples. To the best of our knowledge, the present report describes the first attempt to apply that approach to human testicular data.

A signalling pathway is a cascade of gene activity that changes a biological function in response to signal transduction from outside to inside the cell [11]. Since 1972 [12], such pathways have been one of the major projects in biological research. There had been over 10,000 such signalling pathways described up to 2017 in the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database and are important in the regulation of every type of cellular process [13]. It is logical therefore to use pathway analysis to investigate pathological alterations in gene expression. Pathological processes that result in the ablation of a particular cell type will distort the apparent relative expression levels of genes expressed in the organ as whole, given that different cell types have different profiles of gene expression. Therefore, in such cases, microarrays analyzed at the level of single gene expression are of limited value since it is impossible to determine whether changes are the result of actual alterations in expression or the result of the depletion of a particular cell type. Accordingly, the present work used statistical and computational procedures to analyze microarray data from a mixed cell system to determine alterations in pathways.

2. Materials and Methods

2.1. Microarray Data

Microarray datasets GSE9210 and GSE25518 were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). GSE9210 contains data of 58 testicular biopsy samples, including 47 non-obstructive azoospermia (NOA) samples and 11 obstructive azoospermia (OA) samples. OA was defined as the presence of motile spermatozoa from microsurgical epididymal sperm aspiration (MESA) and testicular sperm extraction (TESE) in the absence of any sperm in the ejaculate. NOA samples have no epididymal testicular or ejaculate spermatozoa [10]. A total of 2,611 transcripts showed statistically significant, differential expression between the two groups [10]. GSE25518 contained 19 IS samples from cryptorchid testes and 4 CT samples from descended, contralateral testes. The paper reported 483 genes absent or under-expressed in the cryptorchid group but their association with the impaired spermatogenesis needed to be confirmed with future work [14]. Comparable methodologies were used in both studies, full details of which are given in the appropriate references [10,14].

2.2. Four-Phase NBI Framework Methodology

The methodology developed for pathway analysis (54) of signalling gene expression in microarray data consists of four phases, summarized as: (1) integrated analysis of differential microarray expression (55, 56); (2) signalling network construction; (3) network analysis for the identification of key biomarkers; (4) corroboration using biological samples. Phase 1 selected those genes in human signalling pathways whose apparent expression levels significantly differ between IS and CT samples. Network degree and centrality analysis was then used to identify and functionally characterize the central nodes [15]. Candidate biomarkers were then selected by comprehensive consideration of their biological functions in pathways and in orthologous functional gene groups (DAVID Bioinformatics resources 6.8) (57-65). Biological simulation was conducted to determine the significance of individual proteins in the regulation of the network and to select and then verify the final target gene. Full details of the methodology are in the supplementary Materials.

3. Results

3.1. Sertoli Cell- and Germ Cell-Specific Genes

The GSE9210 [10] and GSE25518 [14] contain microarray data of 22575 genes and 54675 genes respectively. Although these two microarray datasets were both obtained from biopsied testes, it is not clear from the publications whether the samples contained germ cells. Therefore, the expression of genes specific for Sertoli cells and for germ cells in the two microarray datasets was examined, which revealed the presence of transcripts for 7 Sertoli cell-specific genes and 6 germ cell-specific genes in GSE9210, while GSE25518 showed expression of 10 Sertoli cell-specific and 7 germ cell-specific, genes as identified from PubMed searches (shown in Table S1). The heatmaps of these genes form each database is shown in Figure 1.

Individual genes show differences in expression among individual samples both within and between the two datasets, possibly indicating some variation in cellular composition. However, in general, the germ cell-specific genes (DAZL-HPRT1) in GSE9210 (Figure 1, left-hand panel) show markedly higher expression levels than Sertoli cell-specific, genes in CT samples (x-axis, listed in red) and lower levels in the IS samples (x-axis, listed in black). This suggests that while both the CT and IS samples contain both Sertoli cells and germ cells, the proportion of germ cells in IS is likely to have been lower than in CT. The latter pattern (higher Sertoli-cell gene expression than germ-cell gene expression in IS samples) can also be discerned in the right-hand panel of Figure 1, which shows the equivalent data for GSE25518.

3.2. Network Construction

In total, 1889 genes in the GSE9210 data were found to be significantly differently expressed (p-value <= 0.05 and fold change greater or equal to 2 in the IS compared with the CT samples (Table 1). Of those, 744 were up-regulated and 1145 down-regulated. In the same way, a further 331 genes (298 up regulated and 33 down regulated) were found to be significantly differentially expressed in the GSE25518 data (Table 2). Thus, 2220 apparently significantly differentially expressed genes were found by combining the two gene sets and therefore used in the subsequent pathway analysis. The genes found to be significantly differently expressed in the IS group compared with the CT group are represented in 229 KEGG pathways, in which 51 of the 229 pathways were found to be statistically significant (p <= 0.05) using Fisher’s exact test. Nine of these pathways belong to the cellular processes of the biological signalling pathway class (shown in Table 2), on which we chose to focus as being particularly likely to yield targets for intervention. A total of 161 genes were involved in the 9 pathways and showed significant difference in expression between the IS and CT groups.

A network that combined all of the 9 signalling pathways was built consisting of: 697 genes including 161 genes significantly differentially expressed between CT and IS; 536 downstream genes; and 13 different gene-to-gene interaction types including activation, binding, compound, phosphorylation, dephosphorylation, dissociation, expression, indirect, inhibition, missing, repression, state change and ubiquitination (KEGG database). This network was then modularity analyzed using the software package Gephi (https://gephi.org/) to cluster all genes in the network based on their interactions and to assess their functional characteristics. In a given module, all nodes share a specific function in a pathway process. Based on our network analysis, 16 models were identified, with modularity score of 0.776, which, according to this software package, indicates a high probability of clustering into distinct functional model. The number of genes in the 16 models is given in Table 3.

3.3. Network Analysis

Based on the above network, several network parameters for each node including eccentricity, closeness centrality, betweenness centrality, degree, in-degree, out-degree, clustering coefficient, number of triangles and modularity class were calculated according to the formula given in the methods. Afterwards, 78 genes were extracted from this network based on the above parameter calculation (eccentricity, closeness centrality, betweenness centrality, degree, in-degree, out-degree, clustering coefficient, number of triangles and modularity class) and their result condition (in-degree >= 1 & betweenness centrality >= 1000) (Table S4). These 78 genes showed high centrality in the network so a change in expression of any of them could potentially affect the whole network.

3.4. Target Gene and Male infertility Relative Analysis

A set of functional analyses were performed to validate whether the 78 target genes were positively related to azoospermia. DAVID Bioinformatics resources, the GeneCards database (http://www.genecards.org/) and the Human Mutation Database (http://www.hgmd.cf.ac.uk/ac/index.php) were used to reveal different types of information. DAVID analysis indicated that 11 of those 78 target genes were classified in 9 Reproduction Related clusters and all the other 67 genes were linked with the 11 in the 9 clusters (shown in Table 4).

Based on the 78 genes, a sub-network was then created (Figure 2). The red nodes indicate the associated genes that have previously been reported in the literature to have relevance to azoospermia or male infertility, while the genes in the green box are those identified here for the first time to be associated with azoospermia.

Figure 2. The suspected signalling network for IS based on 78 central genes. The signalling network for IS based on the 78 genes that appear to have different expression levels between IS and CT testes. Gene names in red are genes that have previously been found to be related to azoospermia or male infertility; the genes in the green box are novel candidate genes for involvement in azoospermia. The interaction types between genes were determined from biological pathway databases. Genes in same box are paralogs. Arrows = expression/activation; T-bars = inhibition; straight lines = binding/association; dashed arrow = indirect expression/activation; dashed T-bar = indirect inhibition; empty boxes represent other genes that were not identified by our methods; empty circles indicate other molecules.

Figure 2. The suspected signalling network for IS based on 78 central genes. The signalling network for IS based on the 78 genes that appear to have different expression levels between IS and CT testes. Gene names in red are genes that have previously been found to be related to azoospermia or male infertility; the genes in the green box are novel candidate genes for involvement in azoospermia. The interaction types between genes were determined from biological pathway databases. Genes in same box are paralogs. Arrows = expression/activation; T-bars = inhibition; straight lines = binding/association; dashed arrow = indirect expression/activation; dashed T-bar = indirect inhibition; empty boxes represent other genes that were not identified by our methods; empty circles indicate other molecules.

The 78 genes in the sub-network were analyzed using DAVID, GAD and PubMed to identify genes that could be associated with azoospermia if their expression were perturbed (susceptibility genes). These include the genes that had been reported previously in the literature shown in red in Figure 2. The DAVID and GAD analyses identified 18 genes in the network as susceptible: ITGB1, ITGB3, ITGB4, ITGB5, ITGB6, ITGB7, ITGB8, TP53, CDC20, CDK2, CCNA1, CCNA2, EGFR, FLT1, KDR, MET, PDGFRA and IGF1R. PubMed searches showed that 11 of them had been reported to be related to azoospermia or male infertility [16,17,18,19,20,21,22,23,24,25,26]. The remaining 8 were the ITGB genes, which according to the DAVID GO analysis, are clustered into the subgroup “reproduction”. The ITGB paralog genes have not previously been reported in the literature as related to azoospermia, so they are potentially novel, candidate biomarkers of IS (Figure 2a).

Figure 2a. The suspected signalling network for IS based on 33 proteins. The signalling network shown in Figure 2 showing the proteins, protein groups or isoforms produced by the genes at the nodes in Figure 2. The symbols and annotations are as defined in the key to Figure 2.

Figure 2a. The suspected signalling network for IS based on 33 proteins. The signalling network shown in Figure 2 showing the proteins, protein groups or isoforms produced by the genes at the nodes in Figure 2. The symbols and annotations are as defined in the key to Figure 2.

3.5. Candidate Biomarker Validation

We used leave-one-out support vector machine (SVM), which is a supervised learning algorithm for classification data analysis [27], to check the power of the 7 candidate genes in the classification of IS and CT samples based on the original expression data. The correct rate of classification for GSE9210 data was 97% and for GSE25518 was 83%. The heatmaps for these genes, as shown in Figure 3, demonstrate the differential expression patterns of those 7 genes in CT and IS. The two panels show the original data for the 7 candidate biomarkers in GSE9210 and GSE25518, in which the sample codes in red are the CT samples. A statistically significant up-regulation between IS and CT samples was found, illustrated in Figure 3.

Figure 3. Heatmaps for the SVM results based on 7 candidate genes. The x-axis samples given in red are CT samples, those in black, IS samples. Only 4 of the ITGB genes were represented on GSE9210 and ITGB1, ITGB7 and ITGB8 were clearly up-regulated and ITGB3 slightly up-regulated in the disease samples. Although the results in GSE25518 (right) were not as clear as those from GSE9210, all 7 of the putative biomarkers were upregulated in at least some of the disease samples.

Figure 3. Heatmaps for the SVM results based on 7 candidate genes. The x-axis samples given in red are CT samples, those in black, IS samples. Only 4 of the ITGB genes were represented on GSE9210 and ITGB1, ITGB7 and ITGB8 were clearly up-regulated and ITGB3 slightly up-regulated in the disease samples. Although the results in GSE25518 (right) were not as clear as those from GSE9210, all 7 of the putative biomarkers were upregulated in at least some of the disease samples.

We also used Hierarchical Clustering (HC) to analyze clustering of the samples using the expression data of the 7 candidate genes in the two sets of microarray data, and the results showed a clear separation between the IS and CT samples (Figure 4). Note that the GSE9210 database was derived from approximately twice as many samples as GSE25518 and the differences are qualitatively clearer although significant differences were found for both.

Figure 4. HC results for the 7 candidate biomarkers. Samples listed in red are CT samples, those in black are IS samples. In GSE9210, two CT samples clustered into the IS group while in GSE25518, three CT samples clustered into same group as 13 IS samples and the other CT samples clustered with the remaining 6 IS samples.

Figure 4. HC results for the 7 candidate biomarkers. Samples listed in red are CT samples, those in black are IS samples. In GSE9210, two CT samples clustered into the IS group while in GSE25518, three CT samples clustered into same group as 13 IS samples and the other CT samples clustered with the remaining 6 IS samples.

To investigate further the relevance of the proteins derived from the 78 genes to IS, we performed a computational simulation for the network using the software package Cell Illustrator (The University of Tokyo). Since, as shown in Figure 2a, the sub-network contains genes that produce several closely related proteins or isoforms of the same protein, it only represents 33 types of protein. Of these, 23 were regulated by ITGB including E2F, which represents 5 paralogs. Therefore 19 proteins or protein groups were used in the simulation, which involved changing the effective parameter of each protein type from 0-10 while keeping the parameters of the rest unchanged (Figure 5).

Only altering ITGB had an effect on the majority of other genes in its pathways (Figure 6). The most significant results of the simulation were thus that the genes in ITGB affected the expression levels of the proteins at 18 nodes in the sub-network, which is about 80% of the total (Figure 6). p53 (TP53) influenced about 24% of nodes while the CDK2/CycE and RTK genes affected 21% of nodes in the network, while CDC20 only affected the genes PTTG1 and PTTG2 (data not shown). This simulation result demonstrates that ITGB appears to be the major determinant of the development of IS, and that the higher the expression of ITGB genes, the greater the effect on the network.

Figure 5. Cell illustrator simulation file for the sub-network. This figure displays the cell illustrator simulation file including 33 paralog nodes, which contains 76 target genes and 46 process edges. The number at each node’s bottom left-hand corner is the value of each node, which was changed from 0-10 during each simulation. The top left-hand corner shows the name of each node.

Figure 5. Cell illustrator simulation file for the sub-network. This figure displays the cell illustrator simulation file including 33 paralog nodes, which contains 76 target genes and 46 process edges. The number at each node’s bottom left-hand corner is the value of each node, which was changed from 0-10 during each simulation. The top left-hand corner shows the name of each node.

Figure 6. Simulation results for the effect of ITGB on the proteins in the network whose expression it regulates. This figure displays the simulation curve distribution results for ITGB and the nodes that it affected. In order to use Cell Illustrator to simulate the effect of each gene on the network, the expressed proteins were used to represent each gene or group of paralog genes; a total of 33 proteins are produced by the 78 genes in the sub-network (see figures 2 and 2a). ITGB represents 7 susceptible genes and as shown here, affects the highest number (18) genes or groups of paralogs, six of which have been previously reported as being related to azoospermia.

Figure 6. Simulation results for the effect of ITGB on the proteins in the network whose expression it regulates. This figure displays the simulation curve distribution results for ITGB and the nodes that it affected. In order to use Cell Illustrator to simulate the effect of each gene on the network, the expressed proteins were used to represent each gene or group of paralog genes; a total of 33 proteins are produced by the 78 genes in the sub-network (see figures 2 and 2a). ITGB represents 7 susceptible genes and as shown here, affects the highest number (18) genes or groups of paralogs, six of which have been previously reported as being related to azoospermia.

4. Discussion

Conventional approaches to expression analysis with microarrays use both fold-change in expression and expression p-value [28,29]to identify significantly altered genes. However, this is based on single-gene expression, so if a particular cell type is lost or depleted in a mixed population of cells, the apparent expression level of any gene expressed in other cell types is likely to be higher than it would be in an unaffected cell population. Potentially, this greatly limits the utility of tissue-level expression profiles.

Furthermore, focusing on the expression of individual genes ignores the importance of gene-to-gene interaction. The expression of individual genes influences the expression of others, forming cascades of gene-to-gene interactions that describe biological pathways. These pathways can themselves interact, forming complex networks [30]. Our hypothesis was that the problem of expression-level distortion in organ-based microarray studies could be circumvented by analyzing expression in terms of pathways, using gene to gene co-expression relationships with computational simulation that identifies the driving power of individual target genes/proteins. Thus, although Table 4 shows there was no significant difference in average ITGB1 expression between IS and CT in the two databases, pathway analysis was able to demonstrate a significant increase in expression in IS compared with CT. In the sub-network created, >50% (49 of the 78 genes) are represented in the focal adhesion [31] pathway (Figure 2) with most of the downstream genes leading the pathways to cell proliferation, apoptosis, cell cycle regulation or survival. Furthermore, 25 of the 78 genes are characterized as being in GO groups related to male infertility, and 12 of those 25 are specifically related to azoospermia (David Bioinformatics Resources Database 6.7). To corroborate a functional relationship between any of the 78 potential target genes and azoospermia, the list of network members was compared with 77 azoospermia-related genes identified from GAD (see Materials & Methods). Four genes occurred in both sets: TP53, CCNE1, CCNE2 and CDK2. TP53 is one of the most important genes involved in cellular decision making between DNA repair, cell cycle arrest or apoptosis and has previously been suggested to be involved in idiopathic male infertility [32] . Two of the others, CCNE1 and CCNE2 are key regulators of the cell cycle to activate the fourth gene, CDK2, which forms a complex with CycE [33]. The dysregulation of CDK2/CycE leads to testicular atrophy and a reduction in male fertility because it can arrest spermatogonial mitosis [34] The pathways in the network can be divided into two sets of signal transduction cascades, each headed by a group of membrane-bound genes. One is ITGB, a group of 8 paralog genes, of which 7 members were associated with IS: ITGB1, ITGB3, ITGB4, ITGB5, ITGB6, ITGB7 & ITGB8. The other is RTK, which includes the genes EGFR, FLT1, KDR, MET, PDGFRA & IGF1R, classed as paralogs by the KEGG database. Expression changes in all of the genes in the RTK group have previously been shown to be associated with male infertility including EGFR [35], FLT1 & KDR [23], MET [36], and IGF1R [26]. None of the ITGB paralog genes has been reported in the literature as having any functional relationship with azoospermia or male infertility although DAVID functional annotation allocates ITGB1 into the gamete generation GO group. Its isoform integrin beta-1 is expressed in differentiated prostate gland epithelial cells and regulates primordial germ cell migration [37]. From the present work, the SVM results show a high correct rate in microarray data GSE9210 and the Gephi results (Table 2) for ITGB were also high with 60 in-degree value and over 4000 betweenness centrality values. In addition, the simulation result shows that altering the levels of ITGB would have the greatest effect of all genes in the pathway and around one third of the affected genes have been shown to be directly or indirectly related to azoospermia[10,12,16,18,19,26,38,39,40,41]. This includes genes downstream of ITGB (VCL, ACTB, AKT1, E2F1, TFDP1) as shown in the sub-network in Figure 1, which have all previously been found to a direct effect on azoospermia [10,12,16,18,19,26,38,39,40,41].Together, all these lines of evidence point to the ITGB genes as being key, novel mediators of IS. The work was limited to only two microarray datasets, albeit derived from two different types of pathology, so it is possible that the inclusion of more types could have revealed other pathways. In addition, the computationally derived pathways showing alterations in the microarray data do not explain the staining in interstitial or tubule-wall cells. This may be because the pathways leading to IS are controlled by 7 of the 8 ITGB family members. If fewer ITGB genes were involved in in the tubule wall and interstitial cell events, then the pathway analysis would preferentially identify testicular pathways.

The microarray data identified ITGB1 expression as a key component of the observed gene-expression profile in the IS testicular tissue in the arrays. These findings suggested that altered ITGB1 expression could be associated with testicular changes observed in the microarray dataset and provide a basis for further investigation of its potential contribution in male infertility However, the microarray results alone do not establish a functional role for ITGB1, confirm its value as a diagnostic molecular marker, or demonstrate that modulation of ITGB1 could prevent germ-cell loss or induce infertility for contraceptive purposes. Further work, using direct methods such as immunohistochemistry, is required to confirm that changes in ITGB1 are associated with perturbations of spermatogenesis.

5. Conclusions

This study presents a framework methodology for the analysis of testicular gene expression data, based on the construction of a signalling network of susceptible genes, for the systematic analysis of microarray data. The proposed approach firstly constructs a global map of signaling networks of the differentially expressed genes extracted from microarray data. The central genes of the global signaling network and their interacting genes are then used to create a susceptible signaling network for IS. The computational simulation is finally employed to identify the gene biomarkers. By using this approach, 6 statistical candidate genes for IS were revealed. The integrin-β family genes are particularly interesting as they act as membrane receptors involving cell adhesion and recognition in a variety of processes including embryogenesis and immune responses etc. Defects in the adhesion of germ cells within their micro-environment could thus be an important factor in the development of IS. The HC and SVM results confirm their importance in distinguishing IS from CT. This is the first time such an approach has been used to identify key changes in testicular gene expression associated with impaired spermatogenesis.