A major issue to resolve in this work was the discovery and selection of genes that are affected by AZFc/AZFb + c deletions in a background of transcriptional changes that are dominated by the depletion of specific germ cell types in the different spermatogenic stages when the deletions are present. To eliminate the dominating effect of germ cell-specific transcript levels within the comparison of two (or more) biopsies, it is an indispensable prerequisite to start with samples that present a homogeneous germ cell phenotype. Figure exemplifies this problem using samples from our initial set. When comparing biopsies presenting continuously hypospermatogenesis with those that have an AZFc deletion and a ‘mixed atrophy’ phenotype with predominantly Sertoli-cell-only tubules (upper panel), any statistical filtering method will select genes specific to all germ cells that differ between the two states, i.e., spermatogonia, primary/secondary spermatocytes, and round/elongated spermatids. Furthermore, any clustering method will result in groups defining the subtypes of spermatogenic failure. In contrast, if one employs samples with similar germ cell composition (lower panel), the dominating ‘germ cell effect’ is minimized, and a ‘deletion effect’ might be more apparent.
Fig. 1 Definition of the ‘germ-cell effect’ problem involved in filtering differential genes influenced by Y-chromosomal microdeletions. Based on the semi-thin section micrographs of three samples from our collection (A: MA_8; B: MA_11_c; C: (more ...)
Our initial set of 34 biopsies (Supplemental Data 1, ‘initial set’) had been classified on the basis of germ cell types that are present in the seminiferous tubules. We grouped the top 200 variant genes by hierarchical clustering (Fig. , upper panel). Additional information was acquired about the relative quantity of germ cell composition found in the histological examination of the biopsies (lower panel), illustrating the paradigm shift from the most differentiated germ cell type to that which is most frequent: the separation of groups in the dendrogram correlates not so much with initial histological classification but significantly more with the dominant germ cell type present. For instance, sample SCO_10_c clusters within the MA samples, which tallies with the observation of tubules containing spermatocytes (MA) being the most frequent observation in this sample. A comparison of the dendrogram with the histology for all samples clearly shows (and without exception) that all samples cluster according to their quantitatively prevailing germ cell type. This is an important observation as many of the initial samples with YCMD cluster in groups distinct from their most differentiated germ cell type. To complement these results, we employed principle component analysis (PCA) on the 34 samples using all genes. As PCA is a data reduction technique, one can visualize the first three main variances (principle components) within the gene expression data by three-dimensional display (Supplemental Data 2). Corroborating the dominant effect of relative germ cell composition on the grouping outcome, the main variance of gene expression within the 34 samples is again defined by the most frequent germ cell type (spread along the first principle component PCA1). A similar observation had been made by Ellis et al. (2007
) using a mixture of arrest types with variable etiology. The second and third principle components (PCA2 and PCA3) separate a mixture of germ cell composition and YCMD status. It must be noted, however, that PCA1-3 only explained 65% of the total variance found in the data, and therefore other, yet undescribed effects seem likely. We subjected the top 200 variant genes to GO (gene ontology) term analysis and found a highly significant over-representation of terms such as spermatogenesis (p
= 2E−14), mitosis (p
= 5E−6) and meiosis (p
= 1E−3). These results strengthen once again the fact that germ cell transcription is the dominating effect found in our initial sample set. Furthermore, the findings above indicate that the clustering of gene expression profiles from testicular biopsies is extremely sensitive to the relative germ cell composition, a property that enables us to select homogeneous samples from the clustering outcome.
Fig. 2 Visualizing the appropriate pre-selection of testicular biopsies based on gene expression similarity. Testicular biopsies were classified based on their germ cell composition. The gene expression data of all 34 samples were grouped by hierarchical clustering (more ...)
Based on the grouping of samples in Figs. and we selected 26 samples for further analysis (Supplemental Data 1, ‘selected set’). We clustered the top 200 variant genes from the selected subset of 26 samples by bootstrap hierarchical clustering to acquire some additional information about cluster stabilities (Fig. ). Similar to the previous grouping structure of the initial set of 34 samples, we obtained three main clusters defining the different pathological conditions HYS, MA, and SCO as the most variant genes correlate also here with the germ cell state. In contrast to the initial clustering, samples with YCMD now grouped with the corresponding germ cell phenotypes without deletion and this with relatively high cluster stability. The procedure therefore resulted in a subset in which the difference in germ cell composition between samples of the same arrest type but with absence/presence of YCMD was minimized.
Fig. 3 Hierarchical clustering of the selected sample subset. Biopsy samples with names as indicated above were selected from Fig. based on their homogeneous germ cell composition. The top 200 variant genes were clustered by bootstrap hierarchical (more ...)
We conducted t-tests with subsequent Benjamini–Hochberg multiple testing correction on the MA and SCO phenotypes of this subset between normal/AZFc samples. For the HYS group, we were only able to calculate a ratio between the averaged normal samples and one sample with AZFc due to the fact that possible replicates had been sorted out by the sample selection procedure. We filtered 214 genes with fourfold up/down-regulation for the HYS/HYS-AZFc group, 495 genes with pcorr < 0.05 for the MA/MA-AZFc group, 22 genes with pcorr < 0.05 for the SCO/SCO-AZFc group, and 18 genes with praw < 0.05 for the SCO/SCO-AZFb + c group (Supplemental Data 1, ‘HYS AZFc’, ‘MA AZFc’, ‘SCO AZFc’, and ‘SCO AZFbc’). This way, we obtained lists of differential gene expression between normal samples and those with either AZFc or AZFbc deletions specific to each of the three different spermatogenic failure types and putatively enriched for the ‘deletion effect’. In terms of gene numbers, the largest difference could be attributed to the MA phenotype with 495 genes that account for roughly 3% of genes that were under investigation. A further interest was to obtain a set of genes that show differential expression irrespective of the spermatogenic state. Interesting to note here is that for both HYS and MA, DAZ2 was found to be the transcript with the highest ratio of downregulation (~20-fold), which demonstrates that our gene expression data mirrors the depletion of DAZ transcripts as a consequence of the AZFc deletion. However, the probes used on our array platform were not able to discriminate between the four isoforms of DAZ (DAZ1-4). Furthermore, no overlap was observed between all three failure types, suggesting that the majority of affected genes correlate with specific germ cell expression.
Gene ontology analysis of differential genes within HYS and MA types in the presence of an AZFc deletion (SCO failed to be analyzed due to the small gene number) indicates functional differences in dependence of the phenotype (Table ). Nevertheless, terms such as ‘spermatogenesis’, ‘sexual reproduction’, and ‘meiosis’, although having been decreased by the selection procedure by 2–10 orders of magnitude in p
values, could still be attributed to the fact that our selection procedure was not optimal. We added another step of minimizing the ‘germ cell effect’ by eliminating all genes by an intersection approach from the obtained lists that are differentially expressed throughout human spermatogenesis, based on 2,568 genes obtained from t
-tests between HYS and SCO normal samples, in analogy to Feig et al. (2007
) (Supplemental Data 1, ‘Germ cell genes’). Intersections are given in Fig. . We removed these intersections from the initial HYS and MA gene lists, and reanalyzed in respect to GO term over-representation (Table , terms in bold). This resulted in the elimination of the term ‘sexual reproduction’ from the HYS group and a significant reduction of genes that correlate with the terms ‘spermatogenesis’ and ‘meiosis’ as shown on the increase of p
values at 4–8 orders of magnitude. Several other GO terms such as ‘negative regulation of glucocorticoid synthesis’ and ‘killer cell-mediated cytotoxicity’ appeared after elimination of intersecting genes. An overlap analysis of the reduced sets resulted in seven genes that are up-/downregulated in the presence of AZFc in HYS as well as MA phenotype (Supplemental Data 1, ‘HYS AZFc vs. MA AZFc’).
Over-representation analysis of genes affected by the presence of AZFc or AZFb + c microdeletions
Fig. 4 Venn diagram of gene expression overlap between the filtered groups. Genes affected by the AZFc deletion between the different spermatogenic states (HYS hypospermatogenesis, MA meiotic arrest, SCO Sertoli-cell only) were analyzed in respect to overlap (more ...)
In a last step, we queried the ‘germ cell-depleted’ gene lists on two gene expression datasets of human spermatogenesis that were obtained from a completely different cohort of samples as the analysis of microarray data has the inherent problem of ‘overtraining’ and thus can lack generality (Zervakis et al. 2009
). The remaining 178 genes of the reduced HYS dataset (Supplemental Data 1, ‘HYS wo GSG’) showed no significant trend of expression changes throughout spermatogenesis (Supplemental Data 3, upper panel). The analysis of the reduced MA dataset of 232 genes (Supplemental Data 1, ‘MA wo GSG’) revealed a small increase from SCO to MA that was not statistically significant on the averaged values (Supplemental Data 3, lower panel; Codelink platform: 3.8-fold, p
= 0.27; Affymetrix platform: 2.23-fold; p
= 0.08), which may be interpreted as a remnant and non-removable ‘germ cell effect’ although other interpretations are possible (see “Discussion
A transcriptional deletion map for Y-chromosomal genes residing within the deleted loci was built, demonstrating that the deletions were mirrored almost completely on the transcriptional level, and that our microarray analysis was sensitive enough to reveal the absence of transcripts as a consequence of the deleted genes (Fig. a, b). However, the general impression of Y-chromosomal transcript levels was relatively heterogeneous in the different samples. We also noticed a complete absence of any signal for NLGN4Y transcripts in sample SCO_13_bc in contrast to the second sample with AZFb + c deletion, MA_14_bc (Supplemental Data 1, ‘SCO AZFbc’). To confirm absence of this gene, which is located proximal to the classic AZFb boundary defined by Vogt et al. (1996
), we mapped a possible deletion in this region by using gene-specific primers (Supplemental Data 4) for NLGN4Y, its proximal genes VCY and TMSB4Y, and two STS markers proximal and distal of the NLGN4Y gene (sY95 and sY98). For this sample, the complete chromosomal region for AZFb + c reaching proximally to sY95 was deleted (Fig. c–e), and therefore, NLGN4Y transcripts are undetectable due to the presence of a relatively rare massive AZFb + c deletion (Silber et al. 1998
), also referred to as ‘cytogenetic abnormality’ (Krausz et al. 1999
Fig. 5 Compilation of AZFc/AZFb + c deletion-induced transcriptional changes of Y-chromosomal genes. a Heatmap representation of Y-chromosomal gene expression for the complete subset of 26 samples. b Localization (‘Position’ given (more ...)
To validate our microarray data, we conducted quantitative real-time PCR on selected genes using 13 independent samples with/without deletions. All transcripts were normalized to the S27a transcript, which was shown in recent microarray investigations (compare GEO Database Accessions GSE4797, GSE6023, GSE18997) and also this study (GSE21613) to be highly stable expressed in testicular tissue with different states of spermatogenic impairment. It was therefore selected as a good normalization gene to prevent wrong interpretations of expression changes due to differences in the amount of input cDNA in the real-time PCR reactions (Supplemental Data 5A + B). We confirmed a strong downregulation for HIV-1 Rev-binding protein (HRB; 7.5-fold array, 2.8-fold qPCR) in samples with hypospermatogenesis (HYS_7_c). For this sample, there were several other transcripts with more than tenfold upregulation in the presence of an AZFc deletion (marked in red in Supplemental Data 1), but we failed to verify several of these by qPCR (data not shown). Here, the results for genes with a very high upregulation must be taken critically as we cannot identify hybridization artefacts when replicates are lacking. In contrast to this, other genes with interesting functions in respect to AZF deletions that had prominent down-/upregulation in samples with meiotic arrest or Sertoli-cell-only syndrome (meiotic arrest: breast cancer antigen NY-BR-1, apolipoprotein E; Sertoli-cell only: ring finger protein 126; bradykinin receptor B1) could be validated in respect to the direction of expression changes by qPCR, most probably because replicates were available. The ratios obtained from microarray data and qPCR analysis differed slightly, a common finding that may be attributed to the effect of ‘ratio compression’ for microarray data (Wang et al. 2006