Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

A new paradigm for profiling testicular gene expression during normal and disturbed human spermatogenesis

A new paradigm for profiling testicular gene expression during normal and disturbed human... Abstract The aim of this study was to identify gene expression patterns of the testis that correlate with the appearance of distinct stages of male germ cells. We avoided the pitfalls of mixed pathological phenotypes of the testis and circumvented the inapplicability of using the first spermatogenic wave as done previously on rodents. This was accomplished by using 28 samples showing defined and highly homogeneous pathologies selected from 578 testicular biopsies obtained from 289 men with azoospermia (two biopsies each). The molecular signature of the different developmental stages correlated with the morphological preclassification of the testicular biopsies, as shown by resampling-based hierarchical clustering using different measures of variability. By using analysis of variance (ANOVA) and extensive permutation analysis, we filtered 1181 genes that exhibit exceptional statistical significance in testicular expression and grouped subsets with transcriptional changes within the pre-meiotic (348 genes), post-meiotic (81 genes) and terminal differentiation (38 genes) phase. Several distinct molecular classes, metabolic pathways and transcription factor binding sites are involved, depending on the transcriptional profile of the gene clusters that were built using a novel clustering procedure based on not only similarity but also statistical significance. human spermatogenesis, microarray, male infertility, Johnsen score, germ cell development Introduction Human spermatogenesis is a developmental process that implies drastic changes in cellular morphology and metabolism. These severe alterations include changes in ploidy after meiotic division of diploid primary spermatocytes into haploid secondary spermatocytes and furthermore the structural and nuclear changes needed to develop into mature spermatids, also referred as ‘testicular spermatozoa’. Investigating human germ cell development on a global transcriptional scale might gain insight into the various forms of testicular pathology in terms of spermatogenic status and might also uncover those genes needed for normal and unimpaired germ cell development (reviewed in Wrobel and Primig, 2005). As yet, few studies have used microarray technology to reveal global changes of gene expression in the human testis (Sha et al., 2002; Fox et al., 2003; Rockett et al., 2004). Information on groups of genes was derived from these studies with respect to transcriptional changes in germ cell development, but a major confounding factor is in our opinion the unclear definition of the different developmental stages and thus the varying relative proportions of germ cells in the different pathological states that were used in the microarray experiments, in addition to the diverse types of somatic cells present in the testis. The classification of non-obstructive azoospermia (NOA), for example, comprises a collection of spermatogenic impairments that define the complete and/or focal absence of round spermatids, spermatocytes, spermatogonia or Sertoli cells and consequently a mixture of several phenotypes with an unclear classification. To acquire high-quality data related to testicular transcriptional changes during human spermatogenesis, we set up an experimental paradigm that would hopefully lead to a very precise description of the different testicular states. For this paradigm, two prerequisites seemed indispensable: (i) using a relatively high number of replicates to base the gene filtering procedure on solid statistical ground and (ii) a sharp definition of the spermatogenic status used for the experiments. We consider the latter point as the most crucial step because it addresses the heterogeneity of the sample input and subsequently of the derived data. Similar to the analysis of mutant mouse models, however, the presence or absence of specific germ cell stages in distinct pathologies of the human testis can be exploited to analyse their temporal programme of gene expression. Transcript levels can then be correlated with the presence or absence of distinct germ cell stages and might therefore be involved in a specific series of unique cell-associated processes. These include genetic recombination, meiosis, haploid gene expression, the formation of the acrosome and flagellum, as well as remodelling and condensation of the sperm chromatin, each of these processes requiring novel gene products to occur in a timely, precise and well-coordinated manner (reviewed in Schlecht and Primig, 2003). Thus, in our approach, the selection of histological distinct and homogeneous testicular pathologies was an essential step. Starting with the premise that the appearance of novel germ cell types will have the strongest impact on transcriptional levels, we grouped human testicular biopsies on the basis of the Johnsen score, a classical pathological scoring procedure that classifies the germ cell composition from different stages of spermatogenesis (Johnsen, 1970; Schulze et al., 1999). These different stages, as defined by score counts, imply that the seminiferous tubules progressively lose their cell contents in a specific order, that is, elongated spermatids disappear first, then round spermatids, then spermatocytes, then spermatogonia and finally the Sertoli cells. We took advantage of this observation to investigate histological subtype-specific gene expression patterns of the human testes which might be correlated with the presence or absence of specific germ cell stages in individual biopsies. Materials and methods Patients The study population comprised 28 patients selected from a cohort of 289 men presenting with azoospermia, due to either obstruction of the sperm-transporting ducts or defects in the spermatogenic tissue, who underwent an attempt for testicular sperm extraction (TESE) in the Department of Andrology, University Hospital Hamburg-Eppendorf, Germany, between June 2004 and December 2005. Informed consent and ethic committee approval (OB/X/2000) was obtained, and the studies were conducted in accordance with the guidelines of the Helsinki Declaration regarding the use of human tissues. Morning baseline serum hormone concentrations were measured before testicular biopsy. Two microdeletions of the Y chromosome had been documented for the patients. Testicular biopsy and TESE Testis tissue was obtained at the same time for therapeutic TESE and diagnostic purposes. Tissue samples were taken according to a protocol previously described as the ‘sandwich principle’ (Jezek et al., 1998). Briefly, at least four fragments per testis were processed for histological analysis, post-surgical trial-TESE and cryopreservation. An additional sample of rice grain size (∼30 mg) was submerged in 2 ml of RNAlater® (Ambion, Austin, TX, USA) for RNA extraction and subsequent gene expression profiling. Histology Tissue samples were processed as described (Jezek et al., 1998). Briefly, fragments were immersed in 5.5% glutaraldehyde at 4°C for 2 h and post-fixed in 1% OsO4 for another 2 h. Specimens were dehydrated and embedded in Epon employing standard procedures. One-micrometre semi-thin sections were obtained and stained with Toluidine Blue/pyronine for bright field microscopic evaluation. The following parameters were assessed: number of germ cell-containing tubules, score count of the seminiferous epithelium and spermatozoa obtained at trial-TESE. Cases of maldescensus testis or hypospermatogenesis with so-called mixed atrophy where score count differences between individual tubules were obvious were excluded from subsequent microarray analysis. Target preparation, hybridization, post-hybridization and scanning of oligonucleotide arrays RNAlater®-fixed tissue was homogenized in 3 ml of RNApure™ (Peqlab, Germany) using an Ultra-Turrax™ and total RNA extracted according to the manufacturer’s protocols. Extracts were post-purified on RNeasy™ columns (Qiagen, Germany). RNA integrity (28S/18S ratio) and purity was assessed on an Agilent Bioanalyzer (Model 2100; Agilent Technologies, Palo Alto, CA, USA); only samples with an RNA integrity number (RIN) >7.5 (Agilent software) were employed in the microarray analysis. Amplification and labelling of cRNA was done according to the manufacturer’s protocol (GE Healthcare, Piscataway, NJ, USA) using time-optimized amplification parameters (Spiess et al., 2003). Ten micrograms of cRNA was hybridized to a Codelink™ Human 20K Bioarray (GE Healthcare) containing 19 881 gene-specific probes. Arrays were stringency washed, stained with Cy5™-Streptavidin (GE Healthcare) and washed again according to the manufacturer’s instructions. Slides were dried by centrifugation and immediately scanned on a 428™ Array Scanner (Affymetrix) using Jaguar 2.0 Software. Images were analysed and quantified using the CodeLink™ Expression Analysis Software v4.1 (GE Healthcare) and rigorously quality controlled for hybridization artefacts. Data analysis and statistics All statistical steps were programmed in the open-source statistical scripting language R (http://www.r-project.org). The code is very stripped down and does not perform any error checking but is optimized for speed. All scripts with appropriate descriptions can be downloaded from http://www.dr-spiess.de/R-scripts. Filtering of raw signals and normalization of microarray data Signals were termed significant when the raw fluorescence was above the median + two interquartile ranges calculated from the spike-in signal of 300 negative bacterial control spots. Raw data were log2-transformed and log2-normalized by a cyclic loess approach (Bolstad et al., 2003). The raw microarray data set is found in MIAME-compliance at NCBI GEO Accession GSE4797. The filtered and normalized data set of the relevant genes used for our analysis is included as Supplementary data 1. ANOVA, identification of differentially expressed genes and permutation analysis Analysis of variance (ANOVA) was calculated with the lm function after one-sample Kolmogorov–Smirnov and Shapiro–Wilk testing confirmed a normal distribution of the logarithmized raw data. Genes were arbitrarily termed as significant when their P-value was <10–7, which was the case for 1181 genes. Using BRB Array Tools (http://linus.nci.nih.gov/BRB-ArrayTools.html), which calls compiled Fortran code for speed, extensive permutation analysis was first performed by using 20 000 f-test permutations of the complete data set to calculate the probability of getting this set of genes by chance. Additionally, the P-value was decreased to 10–3 in steps of 10–1, which is a frequent threshold value in the microarray literature. After this global permutation analysis, the significant genes were permutated 20 000 times to calculate a permutation-adjusted P-value (P*). Disclosure of the inherent grouping structure Three different measures of variability on a per-gene basis were used: (i) the Gini coefficient from the R library ineq, a variance measure for income distributions of countries, (ii) the t-statistic from a one-sample t-test over the samples, measuring the deviation from a hypothetical population mean (Tsai et al., 2003) and (iii) a principle component analysis (PCA) on a per-gene basis. The latter method is somewhat similar to the ‘gene shaving’ method described by Hastie et al. (2000), but we used the deviation of the first principle component as the measure instead of the inner product to the first principle component. For agglomerative hierarchical clustering (manhattan distance), the top 200, top 1000 and top 200 genes, respectively, were used, and cluster stability in terms of percentage stability of the tree nodes was evaluated by 2000 random permutations of the class labels via multiscale bootstrap resampling (R library pvclust; Shimodaira, 2004). Clustering of genes using multiple comparison testing We clustered the genes by their outcome from a multiple comparison analysis of the classes (ANOVA post hoc testing). In detail, we used the multcomp library in R, which computes confidence intervals and bounds for the linear combinations of the parameters in a fixed-effects model for unbalanced designs. This is done using family wise error rate protection, so that the probability that all bounds hold is at least 1-alpha, with Tukey’s method adjusted for unequal sample sizes (Tukey–Kramer). The code was written as to transform the usual confidence table output to a result frame, with the significantly differing classes on a per-gene basis as labels, after which genes with the same labels were extracted and grouped (clustering). For more details of the clustering procedure, see Figure 4. Figure 4. Open in new tabDownload slide Statistical workflow for clustering of genes with similar expression patterns. Genes with a P value < 1E–7 obtained from analysis of variance (ANOVA) analysis were subjected to multiple pairwise comparison analysis (Tukey–Kramer test) using the R multcomp package for unbalanced designs. If sample sizes within the calculated groups were high enough for subsequent gene ontology (GO) term and pathway analysis (>100 genes), then these groups were not further modified. If sample sizes were smaller, we filtered a reference gene from within the group by isolating the gene with the minimum sum of P values from all multiple comparisons that were below the threshold P value. We then merged those genes to the reference gene that follow the same expression profile by varying the Pearson correlation values (0.95–0.99). For each of these values we analysed the resulting groups in respect to over-representation and chose those with a clear tendency between the various correlation settings. The inset in the multcomp graphic is the theoretical expression profile deduced from the pairwise comparison results (red line). Note that this profile matches exactly with the isolated reference gene and merged gene group. SymAtlas database, gene ontology and promoter analysis For the bioinformatical analysis of the tissue distribution, SymATLAS batch data (http://symatlas.gnf.org/SymAtlas/) of the 1181-gene set were log2-transformed and summarized as signals across tissues and as the number of genes with their highest expression in a certain tissue. Over-representation analysis of gene ontology (GO) terms and pathways was done with the DAVID 2006 server (http://niaid.abcc.ncifcrf.gov/). Over-representation analysis of transcription factor binding sites (TFBS) was done using the OPOSSUM web tool version 1.3 (http://www.cisreg.ca/cgi-bin/oPOSSUM/opossum) with default settings. Quantitative real-time PCR Quantitative real-time PCR (qRT-PCR) was performed using LightCycler™ (Roche, Basel, Switzerland) technology. cDNAs were synthesized with Superscript™ II reverse transcriptase (Invitrogen, Carlsbad, CA, USA), and PCRs were performed using 10 pmol of each gene-specific primers (Supplementary data 2), 2 µl of dNTP mix (25 mM each, Takara Bio, Shiga, Japan), 0.5 µl of SybrGreen I [1:1000 in dimethylsulphoxide (DMSO); Molecular Probes, Leiden, the Netherlands], 0.25 µl of bovine serum albumin (BSA) (20 mg/ml; Sigma, Germany) and 0.2 µl of Ex-Taq HS (5 U/µl; Takara Bio) in a total volume of 20 µl. Cycling conditions were as follows: denaturation (95°C for 5 min), amplification and quantification (95°C for 10 s, primer-specific annealing temperature for 10 s, 72°C for 30 s with a single fluorescence measurement at the end of the segment) repeated 40 times, a melting curve programme (60–95°C with a heating rate of 0.1°C/s and continuous fluorescence measurement) and a cooling step to 40°C. The threshold cycle (crossing point) was determined by a second derivate maximum method (LightCycler™ Quantification Software). Expression levels were normalized to ribosomal protein S27a expression, showing minimum variation between individual human testis samples of differing pathologies (Figure 6). Fold differences were calculated by use of the REST© software (Pfaffl et al., 2002). The PCR efficiency was calculated to be in the range between 1.65 and 1.81 by an R script that uses a ‘sliding window’ method to estimate the slope of the amplification curve at the point of highest linearity (Ramakers et al., 2003). PCR products were electrophoretically separated on 1.3% agarose/tris-acetic acid-EDTA (TAE) gels and verified by sequence analysis. Figure 6. Open in new tabDownload slide Validation of microarray data by quantitative real-time PCR (qRT-PCR). Two housekeeping genes (S27a and L19) and eight differentially transcribed genes (ACT, HOOK1, Calicin, TSKS, CAPPA3, CROC4, PRKAR2α and AKAP4/AKAP82) were validated by qRT-PCR on the same 28 samples as used for the microarray experiments. Furthermore, three samples (two score 10 and one score 2) from an independent set of individuals were also subjected to qRT-PCR, depicted as red dots on the graph lines. Black line represents expression levels from microarray data corresponding to the left axis with log(2) fluorescence values, and red line represents expression levels from qRT-PCR data corresponding to the right axis with crossing points from the qRT-PCR. Error bars indicate standard deviation from the replicates of the four groups (n = 12, 6, 5, 5). Results Selection of testicular biopsy specimens A modification of the classical scoring procedure by Johnsen (1970) was applied to select biopsy specimens showing a homogeneous histopathology of the testis parenchyma (Schulze et al., 1999). For microarray analysis, we carefully chose 28 specimens from different individuals belonging to a collective of 578 biopsies showing maximum tubular homogeneity as extrapolated from strict histological evaluation of a parallel biopsy taken from the same testis (Figure 1). For final sample selection, results from a diagnostic TESE of the same testis (trial-TESE) were taken into account to exclude any discrepancies between the histological classification in the first biopsy and the spermatogenic activity observed in a second biopsy from another area of the same testis (Figure 2). A quantitative determination of the various parameters in one sample was found to be representative of parallel samples from the same donor testis. Also, only samples yielding high-quality total RNA and target cRNA preparations were included in the analysis. Figure 1. Open in new tabDownload slide Morphological classification used for preselection of testicular samples. Upper panel: Semi-thin section micrographs of representative testis samples that underwent modified Johnsen scoring and were subsequently used for microarray analysis. Lower panel: Diagrammatic and non-quantitative representation of cell types found in testicular biopsies as classified by the morphological criteria. Score 2: Sertoli-cell-only syndrome (SCOS); Score 5: meiotic arrest; Score 8: hypospermatogenesis; Score 10: full spermatogenesis. ES, elongated spermatids/testicular spermatozoa; LC, Leydig cells; Ma, macrophages; MC, myoid cells; PL, preleptotene spermatocytes; PS, pachytene spermatocytes; SG, spermatogonia; RS, round spermatids. White arrows indicate 50 µm scale. Figure 2. Open in new tabDownload slide Representativity of selected biopsies. Kernel density box-plot from 578 testicular biopsies obtained from 289 individuals correlating the histological score from one biopsy with number and search time of testicular spermatozoa found in an independent second biopsy used for testicular sperm extraction (TESE) from the same testis. The area of the kernel density curves correlate with the proportion of samples within the category. Red dots depict the biopsy samples, red numbers are the number versus total number (in brackets) of samples used for the microarray analysis within the category, numbers (in %) on the x-axis represent the percentage of scored samples in relation to the complete clinical population. ATF, average time to find; Sp/FOV, testicular spermatozoa per microscopic field-of-view. Both criteria (FOV, ATF) are parameters derived from microscopical detection of spermatozoa during the TESE process. On the basis of these criteria, replicates belonging to the following subgroups were analysed in parallel: biopsies obtained from patients with post-testicular obstructions showing complete spermatogenesis, i.e. >20 late spermatids (also referred as ‘testicular spermatozoa’) per seminiferous tubule cross-section and repeated occurrence of residual bodies, i.e. score count 10; biopsies showing severe hypospermatogenesis in all seminiferous tubules, i.e. score count 8; biopsies with many spermatocytes and no spermatids, i.e. score count 5; and biopsies from patients with Sertoli cell-only syndrome (SCOS), i.e. score count 2. All patient groups were similar with respect to age and testosterone levels (Table I). Patients with SCOS exhibited elevated gonadotrophin levels, which is the usual case in this group due to lack of or severely reduced levels of inhibin B (Jezek et al., 1998). One patient in the SCOS group and one patient in the score 5 group were carriers of Y chromosomal microdeletions (YCMDs). Table I. Relation between histological evaluation, clinical parameters and pathological characteristics of 28 patients and testicular biopsies used in the microarray experiments Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Open in new tab Table I. Relation between histological evaluation, clinical parameters and pathological characteristics of 28 patients and testicular biopsies used in the microarray experiments Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Open in new tab Identification of differentially expressed genes From the 15 119 genes that had significant signal after background correction and normalization, we filtered a subset for which differential testicular expression throughout the various stages of germ cell development accounts. For this purpose, we followed an ANOVA approach in which the classes were defined by the morphological score that had been acquired from the micrographed testis biopsies. Using the R anova function, we collected a high number of genes with very low P-values that showed a significant differential expression from one class to the other. One of the questions that had to be resolved because of the high number of significant genes was at which level to set the statistical threshold value to term genes as significantly differential or not. Several methods exist in defining such threshold value, such as the widely used false discovery rate (Benjamini and Hochberg, 1995) or the very conservative Bonferroni correction. One of the inherent problems in microarray data is that the assumption of normality and independence of variables does not necessarily hold, so that parametric methods are not always a good choice (Kerr et al., 2000). For this reason, computationally intensive methods such as permutation tests can circumvent the lack of normality and independence and calculate more realistic statistical estimates (Efron, 2000) as shown in the dramatically overoptimistic error rates in gene filtering for class prediction (Ambroise and McLachlan, 2002; Simon et al., 2003). Although we have found a normal distribution in our logarithmized data, using non-parametric permutation methods will always give more realistic statistical estimates. For this reason, we used the free software tool BRB-Array Tools to calculate f-tests (ANOVA) for the gene set, including 20 000 permutations to get an estimate on (i) the chance of obtaining a certain number of genes defined by a specific P-value and (ii) a permutation-adjusted P-value estimate that counts the proportion of permutations with the same or better P-values obtained from the f‐test. We chose a P-value of 10–7 as a sensible threshold value (Supplementary data 3). This value manifested in 1181 genes from the R anova and 1395 genes in the BRB-Array Tools f-test, possibly due to slight differences in algorithm or number rounding. We decided on the former number of genes because subsequent statistical analyses were also conducted in R. Some interesting aspects can be observed here. The permutation P-values for the first 1181 genes ordered by ANOVA P-value were effectively P = 0, resulting from not a single permutation being equivalent or better than the P-values derived by defining the classes. Additionally, there is the clear observation that the permutation P-values rise dramatically at the border of the Bonferroni correction at P = 3.3 × 10–6, and the unadjusted ANOVA P-values are overly optimistic. Disclosure of the grouping structure inherent in the biopsy samples Hierarchical cluster analysis has become an invaluable tool in extracting information about possible structures underlying a microarray data set, with the potential to reveal new molecular distinct cancer subtypes previously undefined (Alizadeh et al., 2000). Interpreting structural tendencies from the remaining filtered genes is often a very important step in discovering clinical and pathological subtypes. We clustered the top 200 genes (with respect to P-value) filtered by ANOVA (Figure 3D) by unsupervised hierarchical clustering. We emphasize here that after the application of a gene-filtering procedure, which is based on the predefinition of a class membership (in this case ANOVA), hierarchical clustering of the filtered data is merely a graphical representation of the preceding statistic. Thus, to confirm that the molecular signature does indeed correlate with our definition of groups based on the morphological scoring, we applied three different measures of variance that do not imply class membership. This rationale has been exemplified by microarray-based analyses from Perou et al. (2000) and Bullinger et al. (2004), who selected genes for unsupervised clustering with a high deviation from the median of all samples. Figure 3. Open in new tabDownload slide Hierarchical clustering (average linkage and manhattan metric) of gene expression profiles from human testicular biopsies. (A) Gini index (top 200 genes). (B) t-score (top 1000 genes). (C) Single-gene principle component analysis (sgPCA) (top 200 genes). (D) Clustering of the gene expression profile obtained by defining four groups [analysis of variance (ANOVA)] as obtained by morphological scoring (top 200 genes). The heat map representation of the top 200 genes for each statistic is shown below the dendrograms. Each hierarchical clustering tree was subjected to 2000 permutations with cluster stability as percent stability of the tree nodes. Colour bars on top of the heat map represent the histological subgroups (Score 2: yellow; Score 5: green; Score 8: blue; Score 10: red). We clustered the top 200 genes with a high Gini index (Figure 3A), which is a measure of inequality previously used for comparing income distributions of countries. Furthermore, hierarchical clustering was applied to the top 1000 genes derived from the t-score of a one-sample t-test, comparing the gene sample mean to a hypothetical normal distribution (Figure 3B). Finally, we conducted a single-gene PCA (sgPCA) and clustered the top 200 genes using the deviation of the first principle component as a figure of merit. All clusterings were submitted to 2000 permutations to evaluate the cluster stability and therefore the significance of the correlation to the morphological scoring. The dendrograms from the top-scoring genes of the three different measures of variability and ANOVA analysis exhibit a highly significant correlation to the morphological classification. This is underlined by the high cluster stability and the observation that all 28 samples cluster independent from the applied statistic into the corresponding morphological class, thus illustrating the molecular relatedness of the samples to the histopathological subgroups. Tissue distribution of the gene set Our filtered gene set was used to query the GNF SymAtlas, a database housing tissue-specific gene expression signatures of mouse and human tissues (Su et al., 2002). The data were collected in batch mode, log2-normalized and summarized with respect to tissue-specific expression. The specificity (Supplementary data 4, upper panel) and distribution of all signals across the various human tissues (Supplementary data 4, lower panel) show that 92% of the known transcripts (which represent 46% of all genes in our filtered data set) has predominant expression in the male gonads. The remaining 8% of genes is expressed in the testis but also has higher expression values in other tissues. Analysis of testicular expression changes during germ cell development by clustering genes on the basis of a multiple comparison analysis ANOVA testing in itself calculates P-values without the information on which classes have a significant difference, meaning that our filtered gene set showed strong differential testicular expression during germ cell development, yet information on the classes that exhibit differential expression was lacking. We applied a multiple comparison analysis (post hoc testing) on the classes on a per-gene basis of our extracted gene set and subsequently grouped those genes with the same outcome. By using this procedure, we were able to collate groups of genes that not only display a similarity in their transcriptional profile throughout the classes but also exhibit a change from one class to the other that is based on statistical significance (Figure 4, left panel). This method of grouping genes is effectively an unsupervised approach, because the numbers of clusters to be built are not selected as in classical supervised approaches of clustering genes [e.g. partitioning around medoids (PAM) or k-means clustering]. We obtained 13 different clusters by selecting a P-value of 0.01, resulting in a differing number of members within each group. The primary criterion for our analysis of testicular subtype-specific expression was that the genes within one group needed to have a statistically significant expression change from one group to another and stay unaltered within the other groups. Exemplified for the most prominent pre-meiotic group (score 2→score 5) with 348 genes, this implies a drastic change in testicular expression when spermatogonia/spermatocytes appear but static expression in the later stages of germ cell development. Proceeding this way, we obtained a second post-meiotic group (score 5→score 8) containing 81 genes and a third group (38 genes) for transcripts specifically expressed during the formation of elongated spermatids/testicular spermatozoa (score 8→score 10). Unfortunately, we made the observation that in wanting to perform over-representation analysis, groups with small number of genes led to unconvincing results. For these groups, we altered the approach by decreasing the stringency of the statistic to increase the group size to an applicable level. This modification was applied because in over-representation analysis, the number of genes needed for gaining a high statistical power increases dramatically when the two proportions (in our case, the proportions of TFBS in our data set versus the background data set) differ only slightly. To illustrate this, we have conducted a power analysis by calculating the number of genes needed in dependence on the difference of the proportions, which can be found in Supplementary data 5. From the existing groups we extracted a ‘reference gene’ that possessed the best overall statistic from the multiple comparison analysis (Figure 4, right panel). We then merged genes to this gene by varying the Pearson correlation coefficient from 0.99 to 0.95 and obtained different group sizes, which we analysed with respect to over-representation. We then kept the biggest group with the same analysis outcome concerning over-representation than more stringent groups, by making the premise that over-representation has an either stabilizing or destabilizing tendency when using higher number of genes. Using this approach, we could increase the group sizes to 100–250 genes that gave reasonable results within the over-representation analysis (Supplementary data 6). This procedure is similar to the analysis by Mansson et al. (2004) for the identification of genes that co-regulate with the Early-B-Cell-Factor (EBF) by filtering genes with a high Pearson correlation. Testicular gene expression patterns specific to the appearance of germ cell subtypes By using the described procedure of gene clustering, we obtained three clusters, with genes discriminating the testicular expression patterns during the pre-meiotic, post-meiotic and terminal differentiation stages (transcriptional profile insets in Figure 5). Over-representation analysis for these three groups was performed with respect to GO terms, pathways and TFBS. Although for approximately half of the genes no biological process is as yet known, distinct functional categories and cellular components could nevertheless be correlated with the different expression profiles. For ∼95% of all differentially expressed genes, the profiles reflected an increase of transcript levels with ongoing germ cell development. Figure 5. Open in new tabDownload slide Circular diagram summarizing differences in protein function, pathways and promoters between the various testicular subtypes. Clustered gene groups were subjected to extensive analysis in respect to gene ontology (GO) terms (biological process), over-represented pathways and over-represented transcription factor binding sites (TFBS) to discriminate the functional differences. Over-representation and axis scales are shown as P values. Circular lines represent gene ontology (GO) terms, blue and red text symbols in the graph are over-represented pathways and transcription factor binding sites, respectively. The most frequent GO categories correlated with germ cell development were ‘mitosis’, ‘spermatogenesis’, ‘transcription’, ‘kinase activity’ and ‘ubiquitin cycle’, depending on the subset of genes examined (Figure 5). ‘Mitosis’ (GO term), ‘Cell cycle’ (pathway) and an over-representation of NF-Y, HLF and Pax6 binding sites were the prevalent characteristics of the pre-meiotic stage. Post-meiotic development was marked by an increase in ‘kinase activity’, ‘cytoskeleton organization’ and ‘transforming growth factor (TGF)-beta signalling’ (GO terms and pathways). Two other different over-represented binding sites were found in this group, namely retinoic acid related orphan receptor (ROR)-alpha and B-cell lineage specific activator protein (BSAP). For the terminal differentiation stage (i.e. spermiogenesis, differentiation of round spermatids to testicular spermatozoa), GO terms and pathways for ‘apoptosis’, ‘ubiquitin cycle’ and ‘glycerophospholipid metabolism’ were markedly increased compared with the previous stages. C-Fos and Freac-4 were the most prevalent TFBS found in this group. One large group was characterized by significantly different expression levels between all scores, meaning that the corresponding transcript levels increased steadily with germ cell development (140 genes, data not shown). No characteristic over-representation could be associated with this profile. Rather, considerable overlap was observed with all other expression profiles. Finally, regarding a small gene subset that was distinct by an expression profile showing significantly up-regulated transcript levels with increasing testicular pathology (i.e. decreasing scores; 70 genes, data not shown), ‘apoptosis’ was the prevailing biological process. Additional characteristic processes were ‘G-protein-coupled receptor protein signalling’ and ‘cell-surface receptor-linked signal transduction’. Also different from all other subsets, for which the nucleus was the major cellular component, ∼50% of the down-regulated gene products had functions associated with membranes. We have additionally validated our microarray data by qRT-PCR on eight genes that showed significant expression changes, of which all showed a high correlation in the obtained transcriptional profiles (Figure 6). To evaluate whether the high concordance between microarray data and qRT-PCR results also holds for independent test samples not included in our reference data set, we also quantified three additional samples (one score 2, two score 10) by qRT-PCR. All three samples had expression levels that correlated with the levels from the according group. There is also high concordance to expression levels found in the literature, but by this validation, we are able to extend the existing data on rodents to the human model. Discussion To describe distinct human testicular pathologies in molecular terms and also to unravel the complex expression patterns characteristic of the normal adult human testis, we categorized 28 testicular biopsies by their ‘expression phenotypes’ which reflected a parallel histological classification of the same testis based on a morphological score count procedure. The paradigm to which we adhered was that the use of testicular biopsies with a clear categorization with respect to cellular composition will give high-quality data on differential testicular expression during germ cell development. Previous investigations on rodents with respect to postnatal and pubertal development of the testis have already gained insight into the expression programme of murine male germ cells (Almstrup et al., 2004; Ellis et al., 2004; Schlecht et al., 2004; Shima et al., 2004). In addition, the isolation of purified germ cells from whole rodent testis revealed interesting genes involved in the various differentiation steps of spermatogenesis (Pang et al., 2006). These experimental setups are evidently not applicable to the analysis of the human, added to the fact that rodent and human spermatogenesis are considerably different concerning the spermatogenic cycle (Heller and Clermont, 1964), organization of seminiferous epithelium (Schulze and Rehder, 1984; Schulze et al., 1986), sperm morphology (Gage, 1998) and endocrine control of spermatogenesis (reviewed in Luetjens et al., 2005). Consequently, in the human, distinct testicular pathologies were exploited for the acquisition of information on germ cell development. Recent studies assessed transcriptional differences between infertile men with NOA and those with obstructive azoospermia (Rockett et al., 2004), SCOS versus OA (Fox et al., 2003) and fetal compared with adult testis (Sha et al., 2002). The inherent problem in such studies is a mixture of several different pathological conditions. By describing testicular pathologies with NOA/OA, it is virtually impossible to pinpoint more detailed pre-meiotic and post-meiotic transcriptional events. In the case of OA, this term comprises three different pathological states (score 10, score 8 and score 10/8), whereas NOA consists of at least 92 different score mixtures. We selected biopsies with no discrepancy between the morphological scoring of the first biopsy and the number of spermatozoa found in a second independent biopsy from the same testis. There is debate whether a biopsy taken from the testis at random is representative (Meng et al., 2000; Turek et al., 2000). From the pool of 578 samples that we investigated, a clear relationship between the histological classification and the trial-TESE finding was observed. There is a slight heterogeneity that prevails in some groups (score 2, score 6 and score 8), but the percentage of total disagreement between both criteria is usually in the region of 10–15%, underlining previous observations (Schulze et al., 1999). Indeed, the tendency for testicular biopsies in our selected cases to be representative for the whole organ is further corroborated by the third biopsy which underwent transcriptional profiling and that clustered exclusively into the histological classification from the first biopsy, even when using measures of variance without preclassification. In our case, all samples clustered with high significance into the according preclassified score group, with score 2 being the class with the highest distance from the other groups. Scores 8 and 10 tended to cluster in near vicinity, which underlines the validity of our filtering approach as these two groups differ only quantitatively in germ cell composition (hypospermatogenesis versus full spermatogenesis). With the analysis of testicular gene expression in the different pathological states, we observed a reduction in the number of scores in comparison with the morphological classification (i.e. no molecular score 9, 7, 6, 4; data not shown) such that five distinct scores can be distinguished by clustering (scores 1, 2, 5, 8 and 10). We believe that by a reduction on a molecular level, this can serve as a more noise-free way of classification, which is less prone to subjective criteria from the investigator. Our approach, in using samples with normal testicular morphology and defined pathological conditions, resulted in a set of 1181 genes with P-values below the Bonferroni correction. This is a characteristic that is seldom found within the microarray literature and only when comparing the transcriptional profiles of different tissues (Rinn et al., 2004). This feature is what we believe the outcome from several factors adding to highly significant statistics: (i) the appearance of new cell types within our classified samples, (ii) the homogeneity of our sample classification leading to a very low intragroup variance in comparison with intergroup variance, (iii) the relatively high number of replicates and (iv) the reproducibility and specificity of our microarray system that exhibits a technical variance, as measured by independent hybridizations with the same RNA, of 2.5% (data not shown). The sensitivity of our system might also be the reason for a slightly higher number of transcripts found in the testis than previously described (see Wrobel and Primig, 2005). Interestingly, the permutation analysis revealed a dramatic underestimation of the ANOVA P-values when higher than the Bonferroni border, an effect that we believe has previously not been shown for an experimental data set. Several features are inherent in our gene set. Forty-six percent of all genes have undefined biological function. The majority (92% of known genes) have their highest expression in testicular tissue, the main functional categories being ‘spermatogenesis’, ‘mitosis’, ‘kinase activity’ and ‘ubiquitin cycle’. Involved pathways include cell cycle, TGF-beta signalling and glycerophospholipid metabolism. As we aimed to analyse the functional differences in various stages of germ cell development, we clustered our gene set by using multiple comparison analysis. Hereby, we obtained 13 groups (P = 0.01) with transcriptional profiles that were based on statistical significant changes from one group to the other. We focused on three groups that characterized the pre-meiotic, post-meiotic and terminal differentiation stage. The numbers of members in the latter two groups were increased by successively decreasing the stringency of the applied statistic to obtain group sizes that gave reasonable results in the over-representation analysis. The terminal differentiation stage is characterized by the smallest number of differential transcripts, which underscores the observation that the spermatozoal transcriptome provides mainly the history of past-transcriptional events. The three groups are very distinct with regard to biological processes, pathways and TFBS. Over-representation analysis revealed biological processes present in all groups, yet with very different weighting, i.e. number of genes with similar function. In the pre-meiotic group, NF-Y binding sites dominate, in this case with the corresponding transcription factor subunit NF-Y alpha (Sinha et al., 1995) being in the same class. Although the androgen receptor is expressed constitutively in our microarray data, two different co-regulators, namely insulin-degrading enzyme (Kupfer et al., 1994) and ras related nuclear protein (RAN) binding protein 9 (Rao et al., 2002), also have their highest expression change in the pre-meiotic stage. The post-meiotic class is marked again by the over-representation of HLF and Pax6 binding sites. Additionally, the orphan receptor ROR-alpha is over-represented in this stage, for which ligands remain to be found; yet, it has been shown that the highly related retinoic acid receptors (RARs) play a crucial role in late germ cell development (Chung et al., 2004). The terminal differentiation stage is characterized by a prevalence of FREAC4 and c-fos binding sites. FREAC4/FOXD1 has a slight increase in expression in testes with late spermatids/testicular spermatozoa but does not fall into our strict statistical criteria (P = 1.1E–4). From the plethora of co-regulators that have been found for c-fos, differential testicular expression during germ cell development as deduced from our 1181 gene list accounts for ELK1, nuclear factor of activated T-cells (NFAT), activating transcription factor 2 (ATF2) and protein kinase C alpha subunit (PRKCa). In summary, we have made the observation that there is a strong concordance in our data set between the over-representation of TFBS and the corresponding transcription factors, or even more often their co-regulators. We think that this undermines our clustering procedure to build the different transcriptional profiles with high statistical stringency. A set of genes (data not shown) had different transcription levels in the samples carrying YCMDs as compared with samples without this chromosomal aberration. These affected genes, however, did not influence the molecular classification and clustering outcomes of the biopsies, because the overlap with our set of 1181 developmental genes was only minor. In both cases included in our study, the chromosomal aberration caused a characteristic testicular histology which was by far the dominating factor. To re-evaluate the effects of these two samples on the data set, we conducted the same ANOVA approach as already applied on the data set but left the two samples out. If the samples with the Y chromosomal deletions had a confounding effect due to aberrant gene expression, one would assume an increase in statistical strength. We observed, however, that by leaving the two samples out, the statistical strength decreased dramatically, that is, P-values increased 2–3 orders of magnitude, an observation that can be induced by leaving any two samples out from the analysis. From the 1181 genes included in our analysis, this was the case for 1173 genes. Additionally, we analysed the two groups (score 2 and score 5) with outlier tests and found no prevalence of aberrant gene expression in the two samples with YCMD. The complete analysis can be found in Supplementary data 1. One must emphasize here that an inherent problem in investigating testicular expression changes is the cellular complexity of the organ. In our experimental setup, we analysed the transcriptional changes found in a complete organ because of the presence/absence of distinct germ cell types. One advantage is that we revealed complex transcriptional changes related to the whole testis during germ cell differentiation. The same point has the inherent disadvantage, compared with isolated cell fractions, that we cannot directly pinpoint the locus of expression change and thus cannot directly differentiate between primary expression change (due to the appearance of novel germ cell subtypes) and secondary effects (transcriptional changes in non-germ cell types, i.e. Sertoli cells or Leydig cells) as a consequence of the appearance. We believe that the majority of differentially expressed genes pertains to the appearance of germ cells because (i) most of the genes have a near background expression in SCOS and (ii) by literature search of the top 100 genes, we found a high occurrence of germ-cell-specific transcripts. Our qRT-PCR validation of several genes on the original ‘training’ set and three samples from an independent ‘test’ set confirmed and extended the expression data shown in rodents and human, with additional information on novel genes being expressed during spermatogenesis (i.e. CROC4; Jeffrey et al. 2000). We also observed transcriptional changes in our validated qRT-PCR data that are evident in the earlier stages of spermatogenesis than previously described. This supports the observation that we and others (Hild et al., 2003) have made, that in terms of sensitivity, microarrays are superior to in situ hybridizations. We hope that our high-quality data set on transcriptional changes under normal and defined pathological conditions of the human testis will lead to the identification of common and species-specific mechanisms in mammalian spermatogenesis and contribute to the elucidation of human male fertility on a molecular level. Supplementary material Supplementary data are available at http://molehr.oxfordjournals.org/. Supplementary data 1: Filtered, log-transformed and cyclic-loess normalized gene expression data sorted by ANOVA P-values and used for the statistical analysis and interpretation in this study. This also includes the results from the statistical testing and clustering procedures. Microsoft Excel spreadsheet. Description of column labels at end of data. Supplementary data 2: Characteristics of the primers used for validation by qRT-PCR. Microsoft Excel spreadsheet. Supplementary data 3: Permutation f-test (ANOVA) results from BRB-Array Tools output. HTML file. Supplementary data 4: Tissue distribution summary of the obtained differential gene expression set. Upper panel: histogram showing the number of genes with their highest expression in a specific tissue. Lower panel: box-plot of the relative expression strength of the gene set throughout the various human tissues. Data for both representations were extracted and summarized from the SymAtlas database (http://symatlas.gnf.org/SymAtlas/) using the batch function. Supplementary data 5: Power analysis based on 2000 random data sets for Fisher exact tests used in over-representation analysis in dependence of the number of observations in the test set. Analysis was based on a proportion of 1% of a certain TFBS in the background data set. Power analysis was done with different proportions (1–16%) of the same TFBS in the test set. Dashed line at P = 0.8 is the common minimal power to be accepted. Supplementary data 6: Overrepresentation analysis for GO terms, KEGG pathways and transcription binding sites (TFBS). Data are P-values of the results in dependence of the number of genes obtained from the Pearson correlation gene merging procedure. Numbers in red were used for Figure 5. Acknowledgements We thank Ms Stefanie Reinhardt for assistance in RNA isolation and Drs Stefan Hartung, Heike Obermann and Martina Behnen for helpful discussions. We also thank Dr Wilma Bilger, Serono, Germany, for generous financial support. This research was supported by the Deutsche Forschungsgemeinschaft grant Sp721/1-1. References Alizadeh AA , Eisen MB, Davis RE, Ma C, Lossos IS, Rosenwald A, Boldrick JC, Sabet H, Tran T, Yu X et al. ( 2000 ) Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature 403 , 503 –511. Almstrup K , Nielsen JE, Hansen MA, Tanaka M, Skakkebaek NE and Leffers H ( 2004 ) Analysis of cell-type-specific gene expression during mouse spermatogenesis. Biol Reprod 70 , 1751 –1761. Ambroise C and McLachlan GJ ( 2002 ) Selection bias in gene extraction on the basis of microarray gene-expression data. Proc Natl Acad Sci USA 99 , 6562 –6566. Benjamini Y and Hochberg Y ( 1995 ) Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 57 , 289 –300. Bolstad BM , Irizarry RA, Astrand M and Speed TP ( 2003 ) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19 , 185 –193. Bullinger L , Dohner K, Bair E, Frohling S, Schlenk RF, Tibshirani R, Dohner H and Pollack JR ( 2004 ) Use of gene-expression profiling to identify prognostic subclasses in adult acute myeloid leukemia. N Engl J Med 350 , 1605 –1616. Chung SS , Sung W, Wang X and Wolgemuth DJ ( 2004 ) Retinoic acid receptor alpha is required for synchronization of spermatogenic cycles and its absence results in progressive breakdown of the spermatogenic process. Dev Dyn 230 , 754 –766. Efron B ( 2000 ) The bootstrap and modern statistics. J Am Stat Assoc 95 , 1293 –1296. Ellis PJ , Furlong RA, Wilson A, Morris S, Carter D, Oliver G, Print C, Burgoyne PS, Loveland KL and Affara NA ( 2004 ) Modulation of the mouse testis transcriptome during postnatal development and in selected models of male infertility. Mol Hum Reprod 10 , 271 –281. Fox MS , Ares VX, Turek PJ, Haqq C and Reijo-Pera RA ( 2003 ) Feasibility of global gene expression analysis in testicular biopsies from infertile men. Mol Reprod Dev 66 , 403 –421. Gage MJ ( 1998 ) Mammalian sperm morphometry. Proc Biol Sci 265 , 97 –103. Hastie T , Tibshirani R, Eisen MB, Alizadeh A, Levy R, Staudt L, Chan WC, Botstein D and Brown P ( 2000 ) ‘Gene shaving’ as a method for identifying distinct sets of genes with similar expression patterns. Genome Biol 1 , RESEARCH0003. Heller CH and Clermont Y ( 1964 ) Kinetics of the germinal epithelium in man. Recent Prog Horm Res 20 , 545 –575. Hild M , Beckmann B, Haas SA, Koch B, Solovyev V, Busold C, Fellenberg K, Boutros M, Vingron M, Sauer F et al. ( 2003 ) An integrated gene annotation and transcriptional profiling approach towards the full gene content of the Drosophila genome. Genome Biol 5 , R3 . Jeffrey PL , Capes-Davis A, Dunn JM, Tolhurst O, Seeto G, Hannan AJ and Lin SL ( 2000 ) CROC-4: a novel brain specific transcriptional activator of c-fos expressed from proliferation through to maturation of multiple neuronal cell types. Mol Cell Neurosci 16 , 185 –196. Jezek D , Knuth UA and Schulze W ( 1998 ) Successful testicular sperm extraction (TESE) in spite of high serum follicle stimulating hormone and azoospermia: correlation between testicular morphology, TESE results, semen analysis and serum hormone values in 103 infertile men. Hum Reprod 13 , 1230 –1234. Johnsen SG ( 1970 ) Testicular biopsy score count—a method for registration of spermatogenesis in human testes: normal values and results in 335 hypogonadal males. Hormones 1 , 2 –25. Kerr MK , Martin M and Churchill GA ( 2000 ) Analysis of variance for gene expression microarray data. J Comput Biol 7 , 819 –837. Kupfer SR , Wilson EM and French FS ( 1994 ) Androgen and glucocorticoid receptors interact with insulin degrading enzyme. J Biol Chem 269 , 20622 –20628. Luetjens CM , Weinbauer GF and Wistuba J ( 2005 ) Primate spermatogenesis: new insights into comparative testicular organisation, spermatogenic efficiency and endocrine control. Biol Rev Camb Philos Soc 80 , 475 –488. Mansson R , Tsapogas P, Akerlund M, Lagergren A, Gisler R and Sigvardsson M ( 2004 ) Pearson correlation analysis of microarray data allows for the identification of genetic targets for early B-cell factor. J Biol Chem 279 , 17905 –17913. Meng MV , Cha I, Ljung BM and Turek PJ ( 2000 ) Relationship between classic histological pattern and sperm findings on fine needle aspiration map in infertile men. Hum Reprod 15 , 1973 –1977. Pang AL , Johnson W, Ravindranath N, Dym M, Rennert OM and Chan WY ( 2006 ) Expression profiling of purified male germ cells: stage-specific expression patterns related to meiosis and postmeiotic development. Physiol Genomics 24 , 75 –85. Perou CM , Sorlie T, Eisen MB, van de Rijn M, Jeffrey SS, Rees CA, Pollack JR, Ross DT, Johnsen H, Akslen LA et al. ( 2000 ) Molecular portraits of human breast tumours. Nature 406 , 747 –752. Pfaffl MW , Horgan GW and Dempfle L ( 2002 ) Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res 30 , e36 . Ramakers C , Ruijter JM, Deprez RH and Moorman AF ( 2003 ) Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data. Neurosci Lett 339 , 62 –66. Rao MA , Cheng H, Quayle AN, Nishitani H, Nelson CC and Rennie PS ( 2002 ) RanBPM, a nuclear protein that interacts with and regulates transcriptional activity of androgen receptor and glucocorticoid receptor. J Biol Chem 277 , 48020 –48027. Rinn JL , Rozowsky JS, Laurenzi IJ, Petersen PH, Zou K, Zhong W, Gerstein M and Snyder M ( 2004 ) Major molecular differences between mammalian sexes are involved in drug metabolism and renal function. Dev Cell 6 , 791 –800. Rockett JC , Patrizio P, Schmid JE, Hecht NB and Dix DJ ( 2004 ) Gene expression patterns associated with infertility in humans and rodent models. Mutat Res 549 , 225 –240. Schlecht U and Primig M ( 2003 ) Mining meiosis and gametogenesis with DNA microarrays. Reproduction 125 , 447 –456. Schlecht U , Demougin P, Koch R, Hermida L, Wiederkehr C, Descombes P, Pineau C, Jegou B and Primig M ( 2004 ) Expression profiling of mammalian male meiosis and gametogenesis identifies novel candidate genes for roles in the regulation of fertility. Mol Biol Cell 15 , 1031 –1043. Schulze W and Rehder U ( 1984 ) Organization and morphogenesis of the human seminiferous epithelium. Cell Tissue Res 237 , 395 –407. Schulze W , Riemer M, Rehder U and Hoehne KH ( 1986 ) Computer-aided three-dimensional reconstructions of the arrangement of primary spermatocytes in human seminiferous epithelium. Cell Tissue Res 244 , 1 –8. Schulze W , Thoms F and Knuth UA ( 1999 ) Testicular sperm extraction: comprehensive analysis with simultaneously performed histology in 1418 biopsies from 766 subfertile men. Hum Reprod 14 ( Suppl 1 ), 82 –96. Sha J , Zhou Z, Li J, Yin L, Yang H, Hu G, Luo M, Chan HC and Zhou K ( 2002 ) Identification of testis development and spermatogenesis-related genes in human and mouse testes using cDNA arrays. Mol Hum Reprod 8 , 511 –517. Shima JE , McLean DJ, McCarrey JR and Griswold MD ( 2004 ) The murine testicular transcriptome: characterizing gene expression in the testis during the progression of spermatogenesis. Biol Reprod 71 , 319 –330. Shimodaira H ( 2004 ) Approximately unbiased tests of regions using multistep-multiscale bootstrap resampling. Ann Stat 32 , 2616 –2641. Simon R , Radmacher MD, Dobbin K and McShane LM ( 2003 ) Pitfalls in the use of DNA microarray data for diagnostic and prognostic classification. J Natl Cancer Inst 95 , 14 –18. Sinha S , Maity SN, Lu J and de Crombrugghe B ( 1995 ) Recombinant rat CBF-C, the third subunit of CBF/NFY, allows formation of a protein-DNA complex with CBF-A and CBF-B and with yeast HAP2 and HAP3. Proc Natl Acad Sci USA 92 , 1624 –1628. Spiess AN , Mueller N and Ivell R ( 2003 ) Amplified RNA degradation in T7-amplification methods results in biased microarray hybridizations. BMC Genomics 4 , 44 . Su AI , Cooke MP, Ching KA, Hakak Y, Walker JR, Wiltshire T, Orth AP, Vega RG, Sapinoso LM, Moqrich A et al. ( 2002 ) Large-scale analysis of the human and mouse transcriptomes. Proc Natl Acad Sci USA 99 , 4465 –4470. Tsai CA , Chen YJ and Chen JJ ( 2003 ) Testing for differentially expressed genes with microarray data. Nucleic Acids Res 31 , e52 . Turek PJ , Ljung BM, Cha I and Conaghan J ( 2000 ) Diagnostic findings from testis fine needle aspiration mapping in obstructed and nonobstructed azoospermic men. J Urol 163 , 1709 –1716. Wrobel G and Primig M ( 2005 ) Mammalian male germ cells are fertile ground for expression profiling of sexual reproduction. Reproduction 129 , 1 –7. Author notes 1Department of Andrology, University Hospital Hamburg-Eppendorf, Hamburg, Germany, 2School of Molecular and Biomedical Science, University of Adelaide, Adelaide, Australia and 3Fertility Center Hamburg, Hamburg, Germany © The Author 2006. Published by Oxford University Press on behalf of the European Society of Human Reproduction and Embryology. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org
 The online version of this article has been published under an open access model. Users are entitled to use, reproduce, disseminate, or display the open access version of this article for non-commercial purposes provided that: the original authorship is properly and fully attributed; the Journal and Oxford University Press are attributed as the original place of publication with the correct citation details given; if an article is subsequently reproduced or disseminated not in its entirety but only in part or as a derivative work this must be clearly indicated. For commercial re-use, please contact journals.permissions@oxfordjournals.org © The Author 2006. Published by Oxford University Press on behalf of the European Society of Human Reproduction and Embryology. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org
 The online version of this article has been published under an open access model. Users are entitled to use, reproduce, disseminate, or display the open access version of this article for non-commercial purposes provided that: the original authorship is properly and fully attributed; the Journal and Oxford University Press are attributed as the original place of publication with the correct citation details given; if an article is subsequently reproduced or disseminated not in its entirety but only in part or as a derivative work this must be clearly indicated. For commercial re-use, please contact journals.permissions@oxfordjournals.org http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Molecular Human Reproduction Oxford University Press

A new paradigm for profiling testicular gene expression during normal and disturbed human spermatogenesis

Download PDF
11 pages

Loading next page...
 
/lp/oxford-university-press/a-new-paradigm-for-profiling-testicular-gene-expression-during-normal-oLnrIicJJr

References (46)

Publisher
Oxford University Press
Copyright
Copyright © 2022 European Society of Human Reproduction and Embryology
ISSN
1360-9947
eISSN
1460-2407
DOI
10.1093/molehr/gal097
pmid
17114209
Publisher site
See Article on Publisher Site

Abstract

Abstract The aim of this study was to identify gene expression patterns of the testis that correlate with the appearance of distinct stages of male germ cells. We avoided the pitfalls of mixed pathological phenotypes of the testis and circumvented the inapplicability of using the first spermatogenic wave as done previously on rodents. This was accomplished by using 28 samples showing defined and highly homogeneous pathologies selected from 578 testicular biopsies obtained from 289 men with azoospermia (two biopsies each). The molecular signature of the different developmental stages correlated with the morphological preclassification of the testicular biopsies, as shown by resampling-based hierarchical clustering using different measures of variability. By using analysis of variance (ANOVA) and extensive permutation analysis, we filtered 1181 genes that exhibit exceptional statistical significance in testicular expression and grouped subsets with transcriptional changes within the pre-meiotic (348 genes), post-meiotic (81 genes) and terminal differentiation (38 genes) phase. Several distinct molecular classes, metabolic pathways and transcription factor binding sites are involved, depending on the transcriptional profile of the gene clusters that were built using a novel clustering procedure based on not only similarity but also statistical significance. human spermatogenesis, microarray, male infertility, Johnsen score, germ cell development Introduction Human spermatogenesis is a developmental process that implies drastic changes in cellular morphology and metabolism. These severe alterations include changes in ploidy after meiotic division of diploid primary spermatocytes into haploid secondary spermatocytes and furthermore the structural and nuclear changes needed to develop into mature spermatids, also referred as ‘testicular spermatozoa’. Investigating human germ cell development on a global transcriptional scale might gain insight into the various forms of testicular pathology in terms of spermatogenic status and might also uncover those genes needed for normal and unimpaired germ cell development (reviewed in Wrobel and Primig, 2005). As yet, few studies have used microarray technology to reveal global changes of gene expression in the human testis (Sha et al., 2002; Fox et al., 2003; Rockett et al., 2004). Information on groups of genes was derived from these studies with respect to transcriptional changes in germ cell development, but a major confounding factor is in our opinion the unclear definition of the different developmental stages and thus the varying relative proportions of germ cells in the different pathological states that were used in the microarray experiments, in addition to the diverse types of somatic cells present in the testis. The classification of non-obstructive azoospermia (NOA), for example, comprises a collection of spermatogenic impairments that define the complete and/or focal absence of round spermatids, spermatocytes, spermatogonia or Sertoli cells and consequently a mixture of several phenotypes with an unclear classification. To acquire high-quality data related to testicular transcriptional changes during human spermatogenesis, we set up an experimental paradigm that would hopefully lead to a very precise description of the different testicular states. For this paradigm, two prerequisites seemed indispensable: (i) using a relatively high number of replicates to base the gene filtering procedure on solid statistical ground and (ii) a sharp definition of the spermatogenic status used for the experiments. We consider the latter point as the most crucial step because it addresses the heterogeneity of the sample input and subsequently of the derived data. Similar to the analysis of mutant mouse models, however, the presence or absence of specific germ cell stages in distinct pathologies of the human testis can be exploited to analyse their temporal programme of gene expression. Transcript levels can then be correlated with the presence or absence of distinct germ cell stages and might therefore be involved in a specific series of unique cell-associated processes. These include genetic recombination, meiosis, haploid gene expression, the formation of the acrosome and flagellum, as well as remodelling and condensation of the sperm chromatin, each of these processes requiring novel gene products to occur in a timely, precise and well-coordinated manner (reviewed in Schlecht and Primig, 2003). Thus, in our approach, the selection of histological distinct and homogeneous testicular pathologies was an essential step. Starting with the premise that the appearance of novel germ cell types will have the strongest impact on transcriptional levels, we grouped human testicular biopsies on the basis of the Johnsen score, a classical pathological scoring procedure that classifies the germ cell composition from different stages of spermatogenesis (Johnsen, 1970; Schulze et al., 1999). These different stages, as defined by score counts, imply that the seminiferous tubules progressively lose their cell contents in a specific order, that is, elongated spermatids disappear first, then round spermatids, then spermatocytes, then spermatogonia and finally the Sertoli cells. We took advantage of this observation to investigate histological subtype-specific gene expression patterns of the human testes which might be correlated with the presence or absence of specific germ cell stages in individual biopsies. Materials and methods Patients The study population comprised 28 patients selected from a cohort of 289 men presenting with azoospermia, due to either obstruction of the sperm-transporting ducts or defects in the spermatogenic tissue, who underwent an attempt for testicular sperm extraction (TESE) in the Department of Andrology, University Hospital Hamburg-Eppendorf, Germany, between June 2004 and December 2005. Informed consent and ethic committee approval (OB/X/2000) was obtained, and the studies were conducted in accordance with the guidelines of the Helsinki Declaration regarding the use of human tissues. Morning baseline serum hormone concentrations were measured before testicular biopsy. Two microdeletions of the Y chromosome had been documented for the patients. Testicular biopsy and TESE Testis tissue was obtained at the same time for therapeutic TESE and diagnostic purposes. Tissue samples were taken according to a protocol previously described as the ‘sandwich principle’ (Jezek et al., 1998). Briefly, at least four fragments per testis were processed for histological analysis, post-surgical trial-TESE and cryopreservation. An additional sample of rice grain size (∼30 mg) was submerged in 2 ml of RNAlater® (Ambion, Austin, TX, USA) for RNA extraction and subsequent gene expression profiling. Histology Tissue samples were processed as described (Jezek et al., 1998). Briefly, fragments were immersed in 5.5% glutaraldehyde at 4°C for 2 h and post-fixed in 1% OsO4 for another 2 h. Specimens were dehydrated and embedded in Epon employing standard procedures. One-micrometre semi-thin sections were obtained and stained with Toluidine Blue/pyronine for bright field microscopic evaluation. The following parameters were assessed: number of germ cell-containing tubules, score count of the seminiferous epithelium and spermatozoa obtained at trial-TESE. Cases of maldescensus testis or hypospermatogenesis with so-called mixed atrophy where score count differences between individual tubules were obvious were excluded from subsequent microarray analysis. Target preparation, hybridization, post-hybridization and scanning of oligonucleotide arrays RNAlater®-fixed tissue was homogenized in 3 ml of RNApure™ (Peqlab, Germany) using an Ultra-Turrax™ and total RNA extracted according to the manufacturer’s protocols. Extracts were post-purified on RNeasy™ columns (Qiagen, Germany). RNA integrity (28S/18S ratio) and purity was assessed on an Agilent Bioanalyzer (Model 2100; Agilent Technologies, Palo Alto, CA, USA); only samples with an RNA integrity number (RIN) >7.5 (Agilent software) were employed in the microarray analysis. Amplification and labelling of cRNA was done according to the manufacturer’s protocol (GE Healthcare, Piscataway, NJ, USA) using time-optimized amplification parameters (Spiess et al., 2003). Ten micrograms of cRNA was hybridized to a Codelink™ Human 20K Bioarray (GE Healthcare) containing 19 881 gene-specific probes. Arrays were stringency washed, stained with Cy5™-Streptavidin (GE Healthcare) and washed again according to the manufacturer’s instructions. Slides were dried by centrifugation and immediately scanned on a 428™ Array Scanner (Affymetrix) using Jaguar 2.0 Software. Images were analysed and quantified using the CodeLink™ Expression Analysis Software v4.1 (GE Healthcare) and rigorously quality controlled for hybridization artefacts. Data analysis and statistics All statistical steps were programmed in the open-source statistical scripting language R (http://www.r-project.org). The code is very stripped down and does not perform any error checking but is optimized for speed. All scripts with appropriate descriptions can be downloaded from http://www.dr-spiess.de/R-scripts. Filtering of raw signals and normalization of microarray data Signals were termed significant when the raw fluorescence was above the median + two interquartile ranges calculated from the spike-in signal of 300 negative bacterial control spots. Raw data were log2-transformed and log2-normalized by a cyclic loess approach (Bolstad et al., 2003). The raw microarray data set is found in MIAME-compliance at NCBI GEO Accession GSE4797. The filtered and normalized data set of the relevant genes used for our analysis is included as Supplementary data 1. ANOVA, identification of differentially expressed genes and permutation analysis Analysis of variance (ANOVA) was calculated with the lm function after one-sample Kolmogorov–Smirnov and Shapiro–Wilk testing confirmed a normal distribution of the logarithmized raw data. Genes were arbitrarily termed as significant when their P-value was <10–7, which was the case for 1181 genes. Using BRB Array Tools (http://linus.nci.nih.gov/BRB-ArrayTools.html), which calls compiled Fortran code for speed, extensive permutation analysis was first performed by using 20 000 f-test permutations of the complete data set to calculate the probability of getting this set of genes by chance. Additionally, the P-value was decreased to 10–3 in steps of 10–1, which is a frequent threshold value in the microarray literature. After this global permutation analysis, the significant genes were permutated 20 000 times to calculate a permutation-adjusted P-value (P*). Disclosure of the inherent grouping structure Three different measures of variability on a per-gene basis were used: (i) the Gini coefficient from the R library ineq, a variance measure for income distributions of countries, (ii) the t-statistic from a one-sample t-test over the samples, measuring the deviation from a hypothetical population mean (Tsai et al., 2003) and (iii) a principle component analysis (PCA) on a per-gene basis. The latter method is somewhat similar to the ‘gene shaving’ method described by Hastie et al. (2000), but we used the deviation of the first principle component as the measure instead of the inner product to the first principle component. For agglomerative hierarchical clustering (manhattan distance), the top 200, top 1000 and top 200 genes, respectively, were used, and cluster stability in terms of percentage stability of the tree nodes was evaluated by 2000 random permutations of the class labels via multiscale bootstrap resampling (R library pvclust; Shimodaira, 2004). Clustering of genes using multiple comparison testing We clustered the genes by their outcome from a multiple comparison analysis of the classes (ANOVA post hoc testing). In detail, we used the multcomp library in R, which computes confidence intervals and bounds for the linear combinations of the parameters in a fixed-effects model for unbalanced designs. This is done using family wise error rate protection, so that the probability that all bounds hold is at least 1-alpha, with Tukey’s method adjusted for unequal sample sizes (Tukey–Kramer). The code was written as to transform the usual confidence table output to a result frame, with the significantly differing classes on a per-gene basis as labels, after which genes with the same labels were extracted and grouped (clustering). For more details of the clustering procedure, see Figure 4. Figure 4. Open in new tabDownload slide Statistical workflow for clustering of genes with similar expression patterns. Genes with a P value < 1E–7 obtained from analysis of variance (ANOVA) analysis were subjected to multiple pairwise comparison analysis (Tukey–Kramer test) using the R multcomp package for unbalanced designs. If sample sizes within the calculated groups were high enough for subsequent gene ontology (GO) term and pathway analysis (>100 genes), then these groups were not further modified. If sample sizes were smaller, we filtered a reference gene from within the group by isolating the gene with the minimum sum of P values from all multiple comparisons that were below the threshold P value. We then merged those genes to the reference gene that follow the same expression profile by varying the Pearson correlation values (0.95–0.99). For each of these values we analysed the resulting groups in respect to over-representation and chose those with a clear tendency between the various correlation settings. The inset in the multcomp graphic is the theoretical expression profile deduced from the pairwise comparison results (red line). Note that this profile matches exactly with the isolated reference gene and merged gene group. SymAtlas database, gene ontology and promoter analysis For the bioinformatical analysis of the tissue distribution, SymATLAS batch data (http://symatlas.gnf.org/SymAtlas/) of the 1181-gene set were log2-transformed and summarized as signals across tissues and as the number of genes with their highest expression in a certain tissue. Over-representation analysis of gene ontology (GO) terms and pathways was done with the DAVID 2006 server (http://niaid.abcc.ncifcrf.gov/). Over-representation analysis of transcription factor binding sites (TFBS) was done using the OPOSSUM web tool version 1.3 (http://www.cisreg.ca/cgi-bin/oPOSSUM/opossum) with default settings. Quantitative real-time PCR Quantitative real-time PCR (qRT-PCR) was performed using LightCycler™ (Roche, Basel, Switzerland) technology. cDNAs were synthesized with Superscript™ II reverse transcriptase (Invitrogen, Carlsbad, CA, USA), and PCRs were performed using 10 pmol of each gene-specific primers (Supplementary data 2), 2 µl of dNTP mix (25 mM each, Takara Bio, Shiga, Japan), 0.5 µl of SybrGreen I [1:1000 in dimethylsulphoxide (DMSO); Molecular Probes, Leiden, the Netherlands], 0.25 µl of bovine serum albumin (BSA) (20 mg/ml; Sigma, Germany) and 0.2 µl of Ex-Taq HS (5 U/µl; Takara Bio) in a total volume of 20 µl. Cycling conditions were as follows: denaturation (95°C for 5 min), amplification and quantification (95°C for 10 s, primer-specific annealing temperature for 10 s, 72°C for 30 s with a single fluorescence measurement at the end of the segment) repeated 40 times, a melting curve programme (60–95°C with a heating rate of 0.1°C/s and continuous fluorescence measurement) and a cooling step to 40°C. The threshold cycle (crossing point) was determined by a second derivate maximum method (LightCycler™ Quantification Software). Expression levels were normalized to ribosomal protein S27a expression, showing minimum variation between individual human testis samples of differing pathologies (Figure 6). Fold differences were calculated by use of the REST© software (Pfaffl et al., 2002). The PCR efficiency was calculated to be in the range between 1.65 and 1.81 by an R script that uses a ‘sliding window’ method to estimate the slope of the amplification curve at the point of highest linearity (Ramakers et al., 2003). PCR products were electrophoretically separated on 1.3% agarose/tris-acetic acid-EDTA (TAE) gels and verified by sequence analysis. Figure 6. Open in new tabDownload slide Validation of microarray data by quantitative real-time PCR (qRT-PCR). Two housekeeping genes (S27a and L19) and eight differentially transcribed genes (ACT, HOOK1, Calicin, TSKS, CAPPA3, CROC4, PRKAR2α and AKAP4/AKAP82) were validated by qRT-PCR on the same 28 samples as used for the microarray experiments. Furthermore, three samples (two score 10 and one score 2) from an independent set of individuals were also subjected to qRT-PCR, depicted as red dots on the graph lines. Black line represents expression levels from microarray data corresponding to the left axis with log(2) fluorescence values, and red line represents expression levels from qRT-PCR data corresponding to the right axis with crossing points from the qRT-PCR. Error bars indicate standard deviation from the replicates of the four groups (n = 12, 6, 5, 5). Results Selection of testicular biopsy specimens A modification of the classical scoring procedure by Johnsen (1970) was applied to select biopsy specimens showing a homogeneous histopathology of the testis parenchyma (Schulze et al., 1999). For microarray analysis, we carefully chose 28 specimens from different individuals belonging to a collective of 578 biopsies showing maximum tubular homogeneity as extrapolated from strict histological evaluation of a parallel biopsy taken from the same testis (Figure 1). For final sample selection, results from a diagnostic TESE of the same testis (trial-TESE) were taken into account to exclude any discrepancies between the histological classification in the first biopsy and the spermatogenic activity observed in a second biopsy from another area of the same testis (Figure 2). A quantitative determination of the various parameters in one sample was found to be representative of parallel samples from the same donor testis. Also, only samples yielding high-quality total RNA and target cRNA preparations were included in the analysis. Figure 1. Open in new tabDownload slide Morphological classification used for preselection of testicular samples. Upper panel: Semi-thin section micrographs of representative testis samples that underwent modified Johnsen scoring and were subsequently used for microarray analysis. Lower panel: Diagrammatic and non-quantitative representation of cell types found in testicular biopsies as classified by the morphological criteria. Score 2: Sertoli-cell-only syndrome (SCOS); Score 5: meiotic arrest; Score 8: hypospermatogenesis; Score 10: full spermatogenesis. ES, elongated spermatids/testicular spermatozoa; LC, Leydig cells; Ma, macrophages; MC, myoid cells; PL, preleptotene spermatocytes; PS, pachytene spermatocytes; SG, spermatogonia; RS, round spermatids. White arrows indicate 50 µm scale. Figure 2. Open in new tabDownload slide Representativity of selected biopsies. Kernel density box-plot from 578 testicular biopsies obtained from 289 individuals correlating the histological score from one biopsy with number and search time of testicular spermatozoa found in an independent second biopsy used for testicular sperm extraction (TESE) from the same testis. The area of the kernel density curves correlate with the proportion of samples within the category. Red dots depict the biopsy samples, red numbers are the number versus total number (in brackets) of samples used for the microarray analysis within the category, numbers (in %) on the x-axis represent the percentage of scored samples in relation to the complete clinical population. ATF, average time to find; Sp/FOV, testicular spermatozoa per microscopic field-of-view. Both criteria (FOV, ATF) are parameters derived from microscopical detection of spermatozoa during the TESE process. On the basis of these criteria, replicates belonging to the following subgroups were analysed in parallel: biopsies obtained from patients with post-testicular obstructions showing complete spermatogenesis, i.e. >20 late spermatids (also referred as ‘testicular spermatozoa’) per seminiferous tubule cross-section and repeated occurrence of residual bodies, i.e. score count 10; biopsies showing severe hypospermatogenesis in all seminiferous tubules, i.e. score count 8; biopsies with many spermatocytes and no spermatids, i.e. score count 5; and biopsies from patients with Sertoli cell-only syndrome (SCOS), i.e. score count 2. All patient groups were similar with respect to age and testosterone levels (Table I). Patients with SCOS exhibited elevated gonadotrophin levels, which is the usual case in this group due to lack of or severely reduced levels of inhibin B (Jezek et al., 1998). One patient in the SCOS group and one patient in the score 5 group were carriers of Y chromosomal microdeletions (YCMDs). Table I. Relation between histological evaluation, clinical parameters and pathological characteristics of 28 patients and testicular biopsies used in the microarray experiments Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Open in new tab Table I. Relation between histological evaluation, clinical parameters and pathological characteristics of 28 patients and testicular biopsies used in the microarray experiments Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Histological evaluation of testis tissue (Johnsen score count) . 10 full spermatogenesis . 8 hypospermatogenesis . 5 meiotic arrests . 2 Sertoli-cell only . Number of samples n = 12 n = 6 n = 5 n = 5 Age of patients [years (mean ± SD)] 39.1 ± 7.0 34.7 ± 3.9 32.5 ± 7.6 30.3 ± 7.0 Serum FSH concentration [IU/ml (mean ± SD)] 5.0 ± 3.2 5.7 ± 3.9 5.6 ± 0.8 10.6 ± 4.2 Serum LH concentration [IU/ml (mean ± SD)] 3.5 ± 1.4 3.9 ± 2.2 3.7 ± 1.8 8.2 ± 2.3 Serum testosterone concentration [ng/ml (mean ± SD)] 5.0 ± 3.2 6.4 ± 4.0 3.4 ± 1.5 4.8 ± 1.6 Y chromosomal microdeletions None None 1 (AZF-c) 1 (AZF-b) Open in new tab Identification of differentially expressed genes From the 15 119 genes that had significant signal after background correction and normalization, we filtered a subset for which differential testicular expression throughout the various stages of germ cell development accounts. For this purpose, we followed an ANOVA approach in which the classes were defined by the morphological score that had been acquired from the micrographed testis biopsies. Using the R anova function, we collected a high number of genes with very low P-values that showed a significant differential expression from one class to the other. One of the questions that had to be resolved because of the high number of significant genes was at which level to set the statistical threshold value to term genes as significantly differential or not. Several methods exist in defining such threshold value, such as the widely used false discovery rate (Benjamini and Hochberg, 1995) or the very conservative Bonferroni correction. One of the inherent problems in microarray data is that the assumption of normality and independence of variables does not necessarily hold, so that parametric methods are not always a good choice (Kerr et al., 2000). For this reason, computationally intensive methods such as permutation tests can circumvent the lack of normality and independence and calculate more realistic statistical estimates (Efron, 2000) as shown in the dramatically overoptimistic error rates in gene filtering for class prediction (Ambroise and McLachlan, 2002; Simon et al., 2003). Although we have found a normal distribution in our logarithmized data, using non-parametric permutation methods will always give more realistic statistical estimates. For this reason, we used the free software tool BRB-Array Tools to calculate f-tests (ANOVA) for the gene set, including 20 000 permutations to get an estimate on (i) the chance of obtaining a certain number of genes defined by a specific P-value and (ii) a permutation-adjusted P-value estimate that counts the proportion of permutations with the same or better P-values obtained from the f‐test. We chose a P-value of 10–7 as a sensible threshold value (Supplementary data 3). This value manifested in 1181 genes from the R anova and 1395 genes in the BRB-Array Tools f-test, possibly due to slight differences in algorithm or number rounding. We decided on the former number of genes because subsequent statistical analyses were also conducted in R. Some interesting aspects can be observed here. The permutation P-values for the first 1181 genes ordered by ANOVA P-value were effectively P = 0, resulting from not a single permutation being equivalent or better than the P-values derived by defining the classes. Additionally, there is the clear observation that the permutation P-values rise dramatically at the border of the Bonferroni correction at P = 3.3 × 10–6, and the unadjusted ANOVA P-values are overly optimistic. Disclosure of the grouping structure inherent in the biopsy samples Hierarchical cluster analysis has become an invaluable tool in extracting information about possible structures underlying a microarray data set, with the potential to reveal new molecular distinct cancer subtypes previously undefined (Alizadeh et al., 2000). Interpreting structural tendencies from the remaining filtered genes is often a very important step in discovering clinical and pathological subtypes. We clustered the top 200 genes (with respect to P-value) filtered by ANOVA (Figure 3D) by unsupervised hierarchical clustering. We emphasize here that after the application of a gene-filtering procedure, which is based on the predefinition of a class membership (in this case ANOVA), hierarchical clustering of the filtered data is merely a graphical representation of the preceding statistic. Thus, to confirm that the molecular signature does indeed correlate with our definition of groups based on the morphological scoring, we applied three different measures of variance that do not imply class membership. This rationale has been exemplified by microarray-based analyses from Perou et al. (2000) and Bullinger et al. (2004), who selected genes for unsupervised clustering with a high deviation from the median of all samples. Figure 3. Open in new tabDownload slide Hierarchical clustering (average linkage and manhattan metric) of gene expression profiles from human testicular biopsies. (A) Gini index (top 200 genes). (B) t-score (top 1000 genes). (C) Single-gene principle component analysis (sgPCA) (top 200 genes). (D) Clustering of the gene expression profile obtained by defining four groups [analysis of variance (ANOVA)] as obtained by morphological scoring (top 200 genes). The heat map representation of the top 200 genes for each statistic is shown below the dendrograms. Each hierarchical clustering tree was subjected to 2000 permutations with cluster stability as percent stability of the tree nodes. Colour bars on top of the heat map represent the histological subgroups (Score 2: yellow; Score 5: green; Score 8: blue; Score 10: red). We clustered the top 200 genes with a high Gini index (Figure 3A), which is a measure of inequality previously used for comparing income distributions of countries. Furthermore, hierarchical clustering was applied to the top 1000 genes derived from the t-score of a one-sample t-test, comparing the gene sample mean to a hypothetical normal distribution (Figure 3B). Finally, we conducted a single-gene PCA (sgPCA) and clustered the top 200 genes using the deviation of the first principle component as a figure of merit. All clusterings were submitted to 2000 permutations to evaluate the cluster stability and therefore the significance of the correlation to the morphological scoring. The dendrograms from the top-scoring genes of the three different measures of variability and ANOVA analysis exhibit a highly significant correlation to the morphological classification. This is underlined by the high cluster stability and the observation that all 28 samples cluster independent from the applied statistic into the corresponding morphological class, thus illustrating the molecular relatedness of the samples to the histopathological subgroups. Tissue distribution of the gene set Our filtered gene set was used to query the GNF SymAtlas, a database housing tissue-specific gene expression signatures of mouse and human tissues (Su et al., 2002). The data were collected in batch mode, log2-normalized and summarized with respect to tissue-specific expression. The specificity (Supplementary data 4, upper panel) and distribution of all signals across the various human tissues (Supplementary data 4, lower panel) show that 92% of the known transcripts (which represent 46% of all genes in our filtered data set) has predominant expression in the male gonads. The remaining 8% of genes is expressed in the testis but also has higher expression values in other tissues. Analysis of testicular expression changes during germ cell development by clustering genes on the basis of a multiple comparison analysis ANOVA testing in itself calculates P-values without the information on which classes have a significant difference, meaning that our filtered gene set showed strong differential testicular expression during germ cell development, yet information on the classes that exhibit differential expression was lacking. We applied a multiple comparison analysis (post hoc testing) on the classes on a per-gene basis of our extracted gene set and subsequently grouped those genes with the same outcome. By using this procedure, we were able to collate groups of genes that not only display a similarity in their transcriptional profile throughout the classes but also exhibit a change from one class to the other that is based on statistical significance (Figure 4, left panel). This method of grouping genes is effectively an unsupervised approach, because the numbers of clusters to be built are not selected as in classical supervised approaches of clustering genes [e.g. partitioning around medoids (PAM) or k-means clustering]. We obtained 13 different clusters by selecting a P-value of 0.01, resulting in a differing number of members within each group. The primary criterion for our analysis of testicular subtype-specific expression was that the genes within one group needed to have a statistically significant expression change from one group to another and stay unaltered within the other groups. Exemplified for the most prominent pre-meiotic group (score 2→score 5) with 348 genes, this implies a drastic change in testicular expression when spermatogonia/spermatocytes appear but static expression in the later stages of germ cell development. Proceeding this way, we obtained a second post-meiotic group (score 5→score 8) containing 81 genes and a third group (38 genes) for transcripts specifically expressed during the formation of elongated spermatids/testicular spermatozoa (score 8→score 10). Unfortunately, we made the observation that in wanting to perform over-representation analysis, groups with small number of genes led to unconvincing results. For these groups, we altered the approach by decreasing the stringency of the statistic to increase the group size to an applicable level. This modification was applied because in over-representation analysis, the number of genes needed for gaining a high statistical power increases dramatically when the two proportions (in our case, the proportions of TFBS in our data set versus the background data set) differ only slightly. To illustrate this, we have conducted a power analysis by calculating the number of genes needed in dependence on the difference of the proportions, which can be found in Supplementary data 5. From the existing groups we extracted a ‘reference gene’ that possessed the best overall statistic from the multiple comparison analysis (Figure 4, right panel). We then merged genes to this gene by varying the Pearson correlation coefficient from 0.99 to 0.95 and obtained different group sizes, which we analysed with respect to over-representation. We then kept the biggest group with the same analysis outcome concerning over-representation than more stringent groups, by making the premise that over-representation has an either stabilizing or destabilizing tendency when using higher number of genes. Using this approach, we could increase the group sizes to 100–250 genes that gave reasonable results within the over-representation analysis (Supplementary data 6). This procedure is similar to the analysis by Mansson et al. (2004) for the identification of genes that co-regulate with the Early-B-Cell-Factor (EBF) by filtering genes with a high Pearson correlation. Testicular gene expression patterns specific to the appearance of germ cell subtypes By using the described procedure of gene clustering, we obtained three clusters, with genes discriminating the testicular expression patterns during the pre-meiotic, post-meiotic and terminal differentiation stages (transcriptional profile insets in Figure 5). Over-representation analysis for these three groups was performed with respect to GO terms, pathways and TFBS. Although for approximately half of the genes no biological process is as yet known, distinct functional categories and cellular components could nevertheless be correlated with the different expression profiles. For ∼95% of all differentially expressed genes, the profiles reflected an increase of transcript levels with ongoing germ cell development. Figure 5. Open in new tabDownload slide Circular diagram summarizing differences in protein function, pathways and promoters between the various testicular subtypes. Clustered gene groups were subjected to extensive analysis in respect to gene ontology (GO) terms (biological process), over-represented pathways and over-represented transcription factor binding sites (TFBS) to discriminate the functional differences. Over-representation and axis scales are shown as P values. Circular lines represent gene ontology (GO) terms, blue and red text symbols in the graph are over-represented pathways and transcription factor binding sites, respectively. The most frequent GO categories correlated with germ cell development were ‘mitosis’, ‘spermatogenesis’, ‘transcription’, ‘kinase activity’ and ‘ubiquitin cycle’, depending on the subset of genes examined (Figure 5). ‘Mitosis’ (GO term), ‘Cell cycle’ (pathway) and an over-representation of NF-Y, HLF and Pax6 binding sites were the prevalent characteristics of the pre-meiotic stage. Post-meiotic development was marked by an increase in ‘kinase activity’, ‘cytoskeleton organization’ and ‘transforming growth factor (TGF)-beta signalling’ (GO terms and pathways). Two other different over-represented binding sites were found in this group, namely retinoic acid related orphan receptor (ROR)-alpha and B-cell lineage specific activator protein (BSAP). For the terminal differentiation stage (i.e. spermiogenesis, differentiation of round spermatids to testicular spermatozoa), GO terms and pathways for ‘apoptosis’, ‘ubiquitin cycle’ and ‘glycerophospholipid metabolism’ were markedly increased compared with the previous stages. C-Fos and Freac-4 were the most prevalent TFBS found in this group. One large group was characterized by significantly different expression levels between all scores, meaning that the corresponding transcript levels increased steadily with germ cell development (140 genes, data not shown). No characteristic over-representation could be associated with this profile. Rather, considerable overlap was observed with all other expression profiles. Finally, regarding a small gene subset that was distinct by an expression profile showing significantly up-regulated transcript levels with increasing testicular pathology (i.e. decreasing scores; 70 genes, data not shown), ‘apoptosis’ was the prevailing biological process. Additional characteristic processes were ‘G-protein-coupled receptor protein signalling’ and ‘cell-surface receptor-linked signal transduction’. Also different from all other subsets, for which the nucleus was the major cellular component, ∼50% of the down-regulated gene products had functions associated with membranes. We have additionally validated our microarray data by qRT-PCR on eight genes that showed significant expression changes, of which all showed a high correlation in the obtained transcriptional profiles (Figure 6). To evaluate whether the high concordance between microarray data and qRT-PCR results also holds for independent test samples not included in our reference data set, we also quantified three additional samples (one score 2, two score 10) by qRT-PCR. All three samples had expression levels that correlated with the levels from the according group. There is also high concordance to expression levels found in the literature, but by this validation, we are able to extend the existing data on rodents to the human model. Discussion To describe distinct human testicular pathologies in molecular terms and also to unravel the complex expression patterns characteristic of the normal adult human testis, we categorized 28 testicular biopsies by their ‘expression phenotypes’ which reflected a parallel histological classification of the same testis based on a morphological score count procedure. The paradigm to which we adhered was that the use of testicular biopsies with a clear categorization with respect to cellular composition will give high-quality data on differential testicular expression during germ cell development. Previous investigations on rodents with respect to postnatal and pubertal development of the testis have already gained insight into the expression programme of murine male germ cells (Almstrup et al., 2004; Ellis et al., 2004; Schlecht et al., 2004; Shima et al., 2004). In addition, the isolation of purified germ cells from whole rodent testis revealed interesting genes involved in the various differentiation steps of spermatogenesis (Pang et al., 2006). These experimental setups are evidently not applicable to the analysis of the human, added to the fact that rodent and human spermatogenesis are considerably different concerning the spermatogenic cycle (Heller and Clermont, 1964), organization of seminiferous epithelium (Schulze and Rehder, 1984; Schulze et al., 1986), sperm morphology (Gage, 1998) and endocrine control of spermatogenesis (reviewed in Luetjens et al., 2005). Consequently, in the human, distinct testicular pathologies were exploited for the acquisition of information on germ cell development. Recent studies assessed transcriptional differences between infertile men with NOA and those with obstructive azoospermia (Rockett et al., 2004), SCOS versus OA (Fox et al., 2003) and fetal compared with adult testis (Sha et al., 2002). The inherent problem in such studies is a mixture of several different pathological conditions. By describing testicular pathologies with NOA/OA, it is virtually impossible to pinpoint more detailed pre-meiotic and post-meiotic transcriptional events. In the case of OA, this term comprises three different pathological states (score 10, score 8 and score 10/8), whereas NOA consists of at least 92 different score mixtures. We selected biopsies with no discrepancy between the morphological scoring of the first biopsy and the number of spermatozoa found in a second independent biopsy from the same testis. There is debate whether a biopsy taken from the testis at random is representative (Meng et al., 2000; Turek et al., 2000). From the pool of 578 samples that we investigated, a clear relationship between the histological classification and the trial-TESE finding was observed. There is a slight heterogeneity that prevails in some groups (score 2, score 6 and score 8), but the percentage of total disagreement between both criteria is usually in the region of 10–15%, underlining previous observations (Schulze et al., 1999). Indeed, the tendency for testicular biopsies in our selected cases to be representative for the whole organ is further corroborated by the third biopsy which underwent transcriptional profiling and that clustered exclusively into the histological classification from the first biopsy, even when using measures of variance without preclassification. In our case, all samples clustered with high significance into the according preclassified score group, with score 2 being the class with the highest distance from the other groups. Scores 8 and 10 tended to cluster in near vicinity, which underlines the validity of our filtering approach as these two groups differ only quantitatively in germ cell composition (hypospermatogenesis versus full spermatogenesis). With the analysis of testicular gene expression in the different pathological states, we observed a reduction in the number of scores in comparison with the morphological classification (i.e. no molecular score 9, 7, 6, 4; data not shown) such that five distinct scores can be distinguished by clustering (scores 1, 2, 5, 8 and 10). We believe that by a reduction on a molecular level, this can serve as a more noise-free way of classification, which is less prone to subjective criteria from the investigator. Our approach, in using samples with normal testicular morphology and defined pathological conditions, resulted in a set of 1181 genes with P-values below the Bonferroni correction. This is a characteristic that is seldom found within the microarray literature and only when comparing the transcriptional profiles of different tissues (Rinn et al., 2004). This feature is what we believe the outcome from several factors adding to highly significant statistics: (i) the appearance of new cell types within our classified samples, (ii) the homogeneity of our sample classification leading to a very low intragroup variance in comparison with intergroup variance, (iii) the relatively high number of replicates and (iv) the reproducibility and specificity of our microarray system that exhibits a technical variance, as measured by independent hybridizations with the same RNA, of 2.5% (data not shown). The sensitivity of our system might also be the reason for a slightly higher number of transcripts found in the testis than previously described (see Wrobel and Primig, 2005). Interestingly, the permutation analysis revealed a dramatic underestimation of the ANOVA P-values when higher than the Bonferroni border, an effect that we believe has previously not been shown for an experimental data set. Several features are inherent in our gene set. Forty-six percent of all genes have undefined biological function. The majority (92% of known genes) have their highest expression in testicular tissue, the main functional categories being ‘spermatogenesis’, ‘mitosis’, ‘kinase activity’ and ‘ubiquitin cycle’. Involved pathways include cell cycle, TGF-beta signalling and glycerophospholipid metabolism. As we aimed to analyse the functional differences in various stages of germ cell development, we clustered our gene set by using multiple comparison analysis. Hereby, we obtained 13 groups (P = 0.01) with transcriptional profiles that were based on statistical significant changes from one group to the other. We focused on three groups that characterized the pre-meiotic, post-meiotic and terminal differentiation stage. The numbers of members in the latter two groups were increased by successively decreasing the stringency of the applied statistic to obtain group sizes that gave reasonable results in the over-representation analysis. The terminal differentiation stage is characterized by the smallest number of differential transcripts, which underscores the observation that the spermatozoal transcriptome provides mainly the history of past-transcriptional events. The three groups are very distinct with regard to biological processes, pathways and TFBS. Over-representation analysis revealed biological processes present in all groups, yet with very different weighting, i.e. number of genes with similar function. In the pre-meiotic group, NF-Y binding sites dominate, in this case with the corresponding transcription factor subunit NF-Y alpha (Sinha et al., 1995) being in the same class. Although the androgen receptor is expressed constitutively in our microarray data, two different co-regulators, namely insulin-degrading enzyme (Kupfer et al., 1994) and ras related nuclear protein (RAN) binding protein 9 (Rao et al., 2002), also have their highest expression change in the pre-meiotic stage. The post-meiotic class is marked again by the over-representation of HLF and Pax6 binding sites. Additionally, the orphan receptor ROR-alpha is over-represented in this stage, for which ligands remain to be found; yet, it has been shown that the highly related retinoic acid receptors (RARs) play a crucial role in late germ cell development (Chung et al., 2004). The terminal differentiation stage is characterized by a prevalence of FREAC4 and c-fos binding sites. FREAC4/FOXD1 has a slight increase in expression in testes with late spermatids/testicular spermatozoa but does not fall into our strict statistical criteria (P = 1.1E–4). From the plethora of co-regulators that have been found for c-fos, differential testicular expression during germ cell development as deduced from our 1181 gene list accounts for ELK1, nuclear factor of activated T-cells (NFAT), activating transcription factor 2 (ATF2) and protein kinase C alpha subunit (PRKCa). In summary, we have made the observation that there is a strong concordance in our data set between the over-representation of TFBS and the corresponding transcription factors, or even more often their co-regulators. We think that this undermines our clustering procedure to build the different transcriptional profiles with high statistical stringency. A set of genes (data not shown) had different transcription levels in the samples carrying YCMDs as compared with samples without this chromosomal aberration. These affected genes, however, did not influence the molecular classification and clustering outcomes of the biopsies, because the overlap with our set of 1181 developmental genes was only minor. In both cases included in our study, the chromosomal aberration caused a characteristic testicular histology which was by far the dominating factor. To re-evaluate the effects of these two samples on the data set, we conducted the same ANOVA approach as already applied on the data set but left the two samples out. If the samples with the Y chromosomal deletions had a confounding effect due to aberrant gene expression, one would assume an increase in statistical strength. We observed, however, that by leaving the two samples out, the statistical strength decreased dramatically, that is, P-values increased 2–3 orders of magnitude, an observation that can be induced by leaving any two samples out from the analysis. From the 1181 genes included in our analysis, this was the case for 1173 genes. Additionally, we analysed the two groups (score 2 and score 5) with outlier tests and found no prevalence of aberrant gene expression in the two samples with YCMD. The complete analysis can be found in Supplementary data 1. One must emphasize here that an inherent problem in investigating testicular expression changes is the cellular complexity of the organ. In our experimental setup, we analysed the transcriptional changes found in a complete organ because of the presence/absence of distinct germ cell types. One advantage is that we revealed complex transcriptional changes related to the whole testis during germ cell differentiation. The same point has the inherent disadvantage, compared with isolated cell fractions, that we cannot directly pinpoint the locus of expression change and thus cannot directly differentiate between primary expression change (due to the appearance of novel germ cell subtypes) and secondary effects (transcriptional changes in non-germ cell types, i.e. Sertoli cells or Leydig cells) as a consequence of the appearance. We believe that the majority of differentially expressed genes pertains to the appearance of germ cells because (i) most of the genes have a near background expression in SCOS and (ii) by literature search of the top 100 genes, we found a high occurrence of germ-cell-specific transcripts. Our qRT-PCR validation of several genes on the original ‘training’ set and three samples from an independent ‘test’ set confirmed and extended the expression data shown in rodents and human, with additional information on novel genes being expressed during spermatogenesis (i.e. CROC4; Jeffrey et al. 2000). We also observed transcriptional changes in our validated qRT-PCR data that are evident in the earlier stages of spermatogenesis than previously described. This supports the observation that we and others (Hild et al., 2003) have made, that in terms of sensitivity, microarrays are superior to in situ hybridizations. We hope that our high-quality data set on transcriptional changes under normal and defined pathological conditions of the human testis will lead to the identification of common and species-specific mechanisms in mammalian spermatogenesis and contribute to the elucidation of human male fertility on a molecular level. Supplementary material Supplementary data are available at http://molehr.oxfordjournals.org/. Supplementary data 1: Filtered, log-transformed and cyclic-loess normalized gene expression data sorted by ANOVA P-values and used for the statistical analysis and interpretation in this study. This also includes the results from the statistical testing and clustering procedures. Microsoft Excel spreadsheet. Description of column labels at end of data. Supplementary data 2: Characteristics of the primers used for validation by qRT-PCR. Microsoft Excel spreadsheet. Supplementary data 3: Permutation f-test (ANOVA) results from BRB-Array Tools output. HTML file. Supplementary data 4: Tissue distribution summary of the obtained differential gene expression set. Upper panel: histogram showing the number of genes with their highest expression in a specific tissue. Lower panel: box-plot of the relative expression strength of the gene set throughout the various human tissues. Data for both representations were extracted and summarized from the SymAtlas database (http://symatlas.gnf.org/SymAtlas/) using the batch function. Supplementary data 5: Power analysis based on 2000 random data sets for Fisher exact tests used in over-representation analysis in dependence of the number of observations in the test set. Analysis was based on a proportion of 1% of a certain TFBS in the background data set. Power analysis was done with different proportions (1–16%) of the same TFBS in the test set. Dashed line at P = 0.8 is the common minimal power to be accepted. Supplementary data 6: Overrepresentation analysis for GO terms, KEGG pathways and transcription binding sites (TFBS). Data are P-values of the results in dependence of the number of genes obtained from the Pearson correlation gene merging procedure. Numbers in red were used for Figure 5. Acknowledgements We thank Ms Stefanie Reinhardt for assistance in RNA isolation and Drs Stefan Hartung, Heike Obermann and Martina Behnen for helpful discussions. We also thank Dr Wilma Bilger, Serono, Germany, for generous financial support. This research was supported by the Deutsche Forschungsgemeinschaft grant Sp721/1-1. References Alizadeh AA , Eisen MB, Davis RE, Ma C, Lossos IS, Rosenwald A, Boldrick JC, Sabet H, Tran T, Yu X et al. ( 2000 ) Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling. Nature 403 , 503 –511. Almstrup K , Nielsen JE, Hansen MA, Tanaka M, Skakkebaek NE and Leffers H ( 2004 ) Analysis of cell-type-specific gene expression during mouse spermatogenesis. Biol Reprod 70 , 1751 –1761. Ambroise C and McLachlan GJ ( 2002 ) Selection bias in gene extraction on the basis of microarray gene-expression data. Proc Natl Acad Sci USA 99 , 6562 –6566. Benjamini Y and Hochberg Y ( 1995 ) Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B 57 , 289 –300. Bolstad BM , Irizarry RA, Astrand M and Speed TP ( 2003 ) A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics 19 , 185 –193. Bullinger L , Dohner K, Bair E, Frohling S, Schlenk RF, Tibshirani R, Dohner H and Pollack JR ( 2004 ) Use of gene-expression profiling to identify prognostic subclasses in adult acute myeloid leukemia. N Engl J Med 350 , 1605 –1616. Chung SS , Sung W, Wang X and Wolgemuth DJ ( 2004 ) Retinoic acid receptor alpha is required for synchronization of spermatogenic cycles and its absence results in progressive breakdown of the spermatogenic process. Dev Dyn 230 , 754 –766. Efron B ( 2000 ) The bootstrap and modern statistics. J Am Stat Assoc 95 , 1293 –1296. Ellis PJ , Furlong RA, Wilson A, Morris S, Carter D, Oliver G, Print C, Burgoyne PS, Loveland KL and Affara NA ( 2004 ) Modulation of the mouse testis transcriptome during postnatal development and in selected models of male infertility. Mol Hum Reprod 10 , 271 –281. Fox MS , Ares VX, Turek PJ, Haqq C and Reijo-Pera RA ( 2003 ) Feasibility of global gene expression analysis in testicular biopsies from infertile men. Mol Reprod Dev 66 , 403 –421. Gage MJ ( 1998 ) Mammalian sperm morphometry. Proc Biol Sci 265 , 97 –103. Hastie T , Tibshirani R, Eisen MB, Alizadeh A, Levy R, Staudt L, Chan WC, Botstein D and Brown P ( 2000 ) ‘Gene shaving’ as a method for identifying distinct sets of genes with similar expression patterns. Genome Biol 1 , RESEARCH0003. Heller CH and Clermont Y ( 1964 ) Kinetics of the germinal epithelium in man. Recent Prog Horm Res 20 , 545 –575. Hild M , Beckmann B, Haas SA, Koch B, Solovyev V, Busold C, Fellenberg K, Boutros M, Vingron M, Sauer F et al. ( 2003 ) An integrated gene annotation and transcriptional profiling approach towards the full gene content of the Drosophila genome. Genome Biol 5 , R3 . Jeffrey PL , Capes-Davis A, Dunn JM, Tolhurst O, Seeto G, Hannan AJ and Lin SL ( 2000 ) CROC-4: a novel brain specific transcriptional activator of c-fos expressed from proliferation through to maturation of multiple neuronal cell types. Mol Cell Neurosci 16 , 185 –196. Jezek D , Knuth UA and Schulze W ( 1998 ) Successful testicular sperm extraction (TESE) in spite of high serum follicle stimulating hormone and azoospermia: correlation between testicular morphology, TESE results, semen analysis and serum hormone values in 103 infertile men. Hum Reprod 13 , 1230 –1234. Johnsen SG ( 1970 ) Testicular biopsy score count—a method for registration of spermatogenesis in human testes: normal values and results in 335 hypogonadal males. Hormones 1 , 2 –25. Kerr MK , Martin M and Churchill GA ( 2000 ) Analysis of variance for gene expression microarray data. J Comput Biol 7 , 819 –837. Kupfer SR , Wilson EM and French FS ( 1994 ) Androgen and glucocorticoid receptors interact with insulin degrading enzyme. J Biol Chem 269 , 20622 –20628. Luetjens CM , Weinbauer GF and Wistuba J ( 2005 ) Primate spermatogenesis: new insights into comparative testicular organisation, spermatogenic efficiency and endocrine control. Biol Rev Camb Philos Soc 80 , 475 –488. Mansson R , Tsapogas P, Akerlund M, Lagergren A, Gisler R and Sigvardsson M ( 2004 ) Pearson correlation analysis of microarray data allows for the identification of genetic targets for early B-cell factor. J Biol Chem 279 , 17905 –17913. Meng MV , Cha I, Ljung BM and Turek PJ ( 2000 ) Relationship between classic histological pattern and sperm findings on fine needle aspiration map in infertile men. Hum Reprod 15 , 1973 –1977. Pang AL , Johnson W, Ravindranath N, Dym M, Rennert OM and Chan WY ( 2006 ) Expression profiling of purified male germ cells: stage-specific expression patterns related to meiosis and postmeiotic development. Physiol Genomics 24 , 75 –85. Perou CM , Sorlie T, Eisen MB, van de Rijn M, Jeffrey SS, Rees CA, Pollack JR, Ross DT, Johnsen H, Akslen LA et al. ( 2000 ) Molecular portraits of human breast tumours. Nature 406 , 747 –752. Pfaffl MW , Horgan GW and Dempfle L ( 2002 ) Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res 30 , e36 . Ramakers C , Ruijter JM, Deprez RH and Moorman AF ( 2003 ) Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data. Neurosci Lett 339 , 62 –66. Rao MA , Cheng H, Quayle AN, Nishitani H, Nelson CC and Rennie PS ( 2002 ) RanBPM, a nuclear protein that interacts with and regulates transcriptional activity of androgen receptor and glucocorticoid receptor. J Biol Chem 277 , 48020 –48027. Rinn JL , Rozowsky JS, Laurenzi IJ, Petersen PH, Zou K, Zhong W, Gerstein M and Snyder M ( 2004 ) Major molecular differences between mammalian sexes are involved in drug metabolism and renal function. Dev Cell 6 , 791 –800. Rockett JC , Patrizio P, Schmid JE, Hecht NB and Dix DJ ( 2004 ) Gene expression patterns associated with infertility in humans and rodent models. Mutat Res 549 , 225 –240. Schlecht U and Primig M ( 2003 ) Mining meiosis and gametogenesis with DNA microarrays. Reproduction 125 , 447 –456. Schlecht U , Demougin P, Koch R, Hermida L, Wiederkehr C, Descombes P, Pineau C, Jegou B and Primig M ( 2004 ) Expression profiling of mammalian male meiosis and gametogenesis identifies novel candidate genes for roles in the regulation of fertility. Mol Biol Cell 15 , 1031 –1043. Schulze W and Rehder U ( 1984 ) Organization and morphogenesis of the human seminiferous epithelium. Cell Tissue Res 237 , 395 –407. Schulze W , Riemer M, Rehder U and Hoehne KH ( 1986 ) Computer-aided three-dimensional reconstructions of the arrangement of primary spermatocytes in human seminiferous epithelium. Cell Tissue Res 244 , 1 –8. Schulze W , Thoms F and Knuth UA ( 1999 ) Testicular sperm extraction: comprehensive analysis with simultaneously performed histology in 1418 biopsies from 766 subfertile men. Hum Reprod 14 ( Suppl 1 ), 82 –96. Sha J , Zhou Z, Li J, Yin L, Yang H, Hu G, Luo M, Chan HC and Zhou K ( 2002 ) Identification of testis development and spermatogenesis-related genes in human and mouse testes using cDNA arrays. Mol Hum Reprod 8 , 511 –517. Shima JE , McLean DJ, McCarrey JR and Griswold MD ( 2004 ) The murine testicular transcriptome: characterizing gene expression in the testis during the progression of spermatogenesis. Biol Reprod 71 , 319 –330. Shimodaira H ( 2004 ) Approximately unbiased tests of regions using multistep-multiscale bootstrap resampling. Ann Stat 32 , 2616 –2641. Simon R , Radmacher MD, Dobbin K and McShane LM ( 2003 ) Pitfalls in the use of DNA microarray data for diagnostic and prognostic classification. J Natl Cancer Inst 95 , 14 –18. Sinha S , Maity SN, Lu J and de Crombrugghe B ( 1995 ) Recombinant rat CBF-C, the third subunit of CBF/NFY, allows formation of a protein-DNA complex with CBF-A and CBF-B and with yeast HAP2 and HAP3. Proc Natl Acad Sci USA 92 , 1624 –1628. Spiess AN , Mueller N and Ivell R ( 2003 ) Amplified RNA degradation in T7-amplification methods results in biased microarray hybridizations. BMC Genomics 4 , 44 . Su AI , Cooke MP, Ching KA, Hakak Y, Walker JR, Wiltshire T, Orth AP, Vega RG, Sapinoso LM, Moqrich A et al. ( 2002 ) Large-scale analysis of the human and mouse transcriptomes. Proc Natl Acad Sci USA 99 , 4465 –4470. Tsai CA , Chen YJ and Chen JJ ( 2003 ) Testing for differentially expressed genes with microarray data. Nucleic Acids Res 31 , e52 . Turek PJ , Ljung BM, Cha I and Conaghan J ( 2000 ) Diagnostic findings from testis fine needle aspiration mapping in obstructed and nonobstructed azoospermic men. J Urol 163 , 1709 –1716. Wrobel G and Primig M ( 2005 ) Mammalian male germ cells are fertile ground for expression profiling of sexual reproduction. Reproduction 129 , 1 –7. Author notes 1Department of Andrology, University Hospital Hamburg-Eppendorf, Hamburg, Germany, 2School of Molecular and Biomedical Science, University of Adelaide, Adelaide, Australia and 3Fertility Center Hamburg, Hamburg, Germany © The Author 2006. Published by Oxford University Press on behalf of the European Society of Human Reproduction and Embryology. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org
 The online version of this article has been published under an open access model. Users are entitled to use, reproduce, disseminate, or display the open access version of this article for non-commercial purposes provided that: the original authorship is properly and fully attributed; the Journal and Oxford University Press are attributed as the original place of publication with the correct citation details given; if an article is subsequently reproduced or disseminated not in its entirety but only in part or as a derivative work this must be clearly indicated. For commercial re-use, please contact journals.permissions@oxfordjournals.org © The Author 2006. Published by Oxford University Press on behalf of the European Society of Human Reproduction and Embryology. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org
 The online version of this article has been published under an open access model. Users are entitled to use, reproduce, disseminate, or display the open access version of this article for non-commercial purposes provided that: the original authorship is properly and fully attributed; the Journal and Oxford University Press are attributed as the original place of publication with the correct citation details given; if an article is subsequently reproduced or disseminated not in its entirety but only in part or as a derivative work this must be clearly indicated. For commercial re-use, please contact journals.permissions@oxfordjournals.org

Journal

Molecular Human Reproduction – Oxford University Press

Published: Jan 1, 2007

Keywords: gene expression; germ cells

There are no references for this article.