Transcriptome profiling of mice testes following low dose irradiation

Background Radiotherapy is used routinely to treat testicular cancer. Testicular cells vary in radio-sensitivity and the aim of this study was to investigate cellular and molecular changes caused by low dose irradiation of mice testis and to identify transcripts from different cell types in the adult testis. Methods Transcriptome profiling was performed on total RNA from testes sampled at various time points (n = 17) after 1 Gy of irradiation. Transcripts displaying large overall expression changes during the time series, but small expression changes between neighbouring time points were selected for further analysis. These transcripts were separated into clusters and their cellular origin was determined. Immunohistochemistry and in silico quantification was further used to study cellular changes post-irradiation (pi). Results We identified a subset of transcripts (n = 988) where changes in expression pi can be explained by changes in cellularity. We separated the transcripts into five unique clusters that we associated with spermatogonia, spermatocytes, early spermatids, late spermatids and somatic cells, respectively. Transcripts in the somatic cell cluster showed large changes in expression pi, mainly caused by changes in cellularity. Further investigations revealed that the low dose irradiation seemed to cause Leydig cell hyperplasia, which contributed to the detected expression changes in the somatic cell cluster. Conclusions The five clusters represent gene expression in distinct cell types of the adult testis. We observed large expression changes in the somatic cell profile, which mainly could be attributed to changes in cellularity, but hyperplasia of Leydig cells may also play a role. We speculate that the possible hyperplasia may be caused by lower testosterone production and inadequate inhibin signalling due to missing germ cells.


Background
Radiotherapy is used routinely to treat testicular cancer. Cells in the testis display different radio-sensitivity and the degree of radiation-induced damage depends on dose and fractionation of the radiation [1]. Some cells are very sensitive and undergo apoptosis after low dose irradiation, others are only partially affected and presumably recover after a period of time, while yet others seem completely unaffected. Spermatogonial stem cells (SSCs) are the most radio-resistant germ cells, probably due to their slow turnover. Contrarily, A spermatogonia have the highest turnover and are the germ cells most susceptible to irradiation-induced damage [2,3] followed by Intermediate and B spermatogonia [4,5]. The more differentiated germ cells such as spermatocytes and spermatids seem unaffected by low dose irradiation and continue their maturation toward spermatozoa [1].
The somatic cells in the testis are also to some extent affected by irradiation. Human Leydig cells are damaged by the irradiation doses used to treat testicular cancer patients (16)(17)(18)(19)(20) [1,6,7], although in many cases the testosterone production is maintained. Rodent Leydig cells have shown dysfunction after high dose irradiation (5 Gy) leading to low testosterone and increased serum luteinizing hormone (LH) levels [8]. Similarly, it has also been reported that rodent Sertoli cells are affected by high dose irradiation (5 Gy) [9].
Low dose irradiation of adult rodent testes kills A through B spermatogonia, which creates a "gap" in the germinal epithelium [10]. The gap moves through the differentiation stages of the germ cells as spermatogenesis continues and SSCs re-populate the testis. The cellular and molecular changes caused by the movement of the gap can be used to study gene expression in adult spermatogenesis. However, by studying the whole testis we investigate changes in the relative gene expression, i.e. the contribution from each cell type to the total RNA pool, which is determined by the concentration of the RNA in the specific cell type and the percentage of the total volume that the cell type occupies -the cellularity [11].
In this study, we identified five clusters of cell typespecific transcripts by transcriptome profiling. The transcripts in these clusters all had gene expression profiles that could be explained by changes in the cellularity during recovery from irradiation. We observed large expression changes in the somatic cell cluster and further investigations by immunohistochemistry (IHC) revealed that low dose irradiation seemed to cause hyperplasia of the Leydig cells, whereas peritubular myoid (PTM) and Sertoli cells seemed unaffected.

Mice testis preparation
Treatment of mice and preparation of testicular RNA and sections were performed as in Shah et al., [10]. In brief, male C3H/He mice from Japan SLC (Shizuoka, Japan) were maintained under controlled conditions (22 ± 2°C, 55 ± 5% humidity, 12 h light/dark cycle, lights on 0600 h) with laboratory chow (CE-2, Japan Crea, Tokyo, Japan) and water ad libitum. Eleven-week old mice were anesthetised with pentobarbital and covered with lead sheeting except for the scrotum. The testes were locally exposed to X-ray irradiation with 1 Gy. Testes from one to four mice were sampled and weighted regularly during recovery on days 0, 3,7,10,14,17,21,24,28,31,35,38,42,45,48,52,56, and 59. One testis was fixed in 4% paraformaldehyde in 0.1 M phosphate buffer, pH 7.4, overnight at 4°C and subsequently dehydrated in graded series of ethanol and embedded in paraffin for IHC. The contralateral testis was snap-frozen in liquid nitrogen for later preparation of total RNA.
The Japanese Pharmacological Society approved the animal study. The animals were treated according to generally accepted guidelines for animal experimentation at St. Marianna University Graduate School of Medicine and guiding principles for the care and use of laboratory animals.

Gene expression microarrays
Total RNA was purified from whole testis samples collected pi day 3 to day 59 (n = 17) with NucleoSpin RNA II (Macherey-Nagel, Dueren, Germany) according to the manufacturer's protocol. RNA quality was determined using Bioanalyzer nano kit (Agilent Technologies, Santa Clara, California, US). The samples were amplified (one round) using the MessageAmp II aRNA Amplification Kit (Applied Biosystems, Carlsbad, California, US) and the aRNA was applied to Agilent whole mouse genome 4 × 44K oligo microarrays. Hybridisation and scanning of the one-colour arrays were done as described by the manufacturer (Agilent Technologies, Santa Clara, California, US). The microarray data analysis was performed in the R software [12] where the gProcessedSignals were loaded into the limma R/Bioconductor package [13,14] and normalised between arrays using the quantile normalisation method [15]. The probes were collapsed for each systematic transcript ID by the median expression value.

Enrichment score
Transcripts with expression profiles potentially explained by changes in cellularity due to the recovery from irradiation were identified and selected for further analysis based on the enrichment score (ES) calculated as follows: where GE i and GE (i + 1) is the expression value for a transcript at time point i and the neighbouring time point i + 1, respectively, and -GE is the mean expression value for a transcript through the time series.
The false discovery rate (FDR) of the ES was calculated for each transcript based on ten calculations on shuffled time points. A subset of transcripts was chosen for further analysis all with a cumulative minimum of the FDR ≤ 30% and a standard deviation (SD) of the most extreme expression value in the time series ≤ 15%.

Cluster analysis
Cluster analyses were performed on the subset of transcripts selected as described above. A distance matrix was generated with the pair-wise correlation variances of the gene expression values during the time series centred around 1. The transcripts were clustered according to the distance matrix using partitioning around medoids (PAM) clustering, which was performed multiple times with different number of clusters. The number of clusters that separated most unique transcript clusters was chosen as the best separation.

Association of transcript clusters to specific testis cells
To determine the cellular origin of the transcript clusters, we compared their gene expression patterns to the expression profiles of cell markers determined in our previous study of the same testis samples where we used differential display and in situ hybridisation [10]. We also mapped testicular cell-specific markers to the clusters to confirm their cellular origin [16][17][18].

Gene set enrichment analysis
Gene set enrichment analysis (GSEA) was performed on the transcripts in each cluster separately using DAVID [19,20]. All genes represented on the array were used as background. The cut off for statistical significance was set to the Bonferroni corrected p-value ≤ 0.01.

Immunohistochemistry (IHC)
The following primary antibodies were used: Vimentin/ HRP (Dako, Glostrup, Denmark; U7034), Transforming growth factor β receptor III (Tgfbr3 Vimentin was used according to the manufacturer's protocol. In short, the sections were deparaffinised, rehydrated and blocked for endogen peroxidase with H 2 O 2 , washed in tap water, placed 5 min in TBS (0.5 M Tris/HCl, 0.15 M NaCl, pH 7.6) at 37°C, incubated with 1:10 Trypsin (Gibson 15400) in TBS 15 min at 37°C, and exposed to the antibody for 1 h at room temperature. Development was performed with 3-amino-9-ethylcarbazole (AEC). Thorough washing with TBS was performed after each individual step and finally the sections were washed in water before a short staining with Meyers haematoxylin.
The three remaining antibodies, against Tgfbr3, Hsd3b and SMA, were used as in a protocol based on a Zymed histostain kit (Invitrogen, Carlsbad, CA, USA). In short, the sections were deparaffinised, rehydrated and blocked for endogen peroxidase as described above, followed by microwave treatment for 15 min in TEG buffer (10 mM Tris, 0.5 mM EGTA, pH 9.0) for Tgfbr3 and Hsd3b, whereas Citrate buffer (10 mM, pH 6) was used for SMA. Cross reactivity of the antibodies was minimized by treatment with 0.5% milk powder diluted in TBS. Sections were exposed to the primary antibodies over night at 5°C and 1 h at room temperature, then incubated with biotinylated goat anti-rabbit IgG or with biotinylated donkey anti-goat IgG 1:400 in TBS (The binding site Ltd., Birmingham, UK; AB360) for Tgfbr3 before exposure to a peroxidase-conjugated streptavidin complex. Finally, the sections were developed with AEC and counter stained with Meyers haematoxylin. Thorough washing with TBS was performed after each individual step. Slides without addition of primary antibodies were used as controls.

In silico quantification of Leydig cells
In order to obtain a measure of the amount of Leydig cells relative to other cells in the testis, we quantified the area covered by Leydig cells relative to the area covered by other testicular cells. Three different areas, each representing more than ten seminiferous tubules, from sections of the testis of representative days pi were stained with the Leydig cell marker Hsd3b, scanned and evaluated in silico with the VisioMorph add-in in the VIS program (Visiopharm A/S, Hørsholm, Denmark). A red colour was assigned to the red Hsd3b staining, blue to other cells (Meyers haematoxylin) and white to the background (see Additional file 1: Figure S1). A Baysian classification of the colour space from the stained sections was applied and subsequently red areas smaller than 20 μm 2 assigned as blue. The areas of red, blue and white were measured and the ratio of red to blue used as a measure of the amount of Leydig cells relative to the amount of other testis cells.

Clusters of transcripts affected by low dose irradiation
Gene expression profiling was performed on total RNA from mice testis (n = 17) from a time series following low dose irradiation (1 Gy). In total, 988 transcripts were selected for further analysis. These transcripts all had expression profiles during the time series that potentially were caused by the recovery from irradiation, i.e. transcripts with large overall expression changes during the time series, but small expression changes between neighbouring time points. We performed PAM clustering of the distance matrix of the pair-wise correlation variances and identified five clusters with unique expression patterns. This was the maximum number of clusters separable in this data set. The names of the transcripts in each cluster are available in additional material (see Additional file 2: Table S1, Additional file 3: Table S2, Additional file 4: Table S3, Additional file 5: Table S4, Additional file 6: Table S5).

Assignment of transcript clusters to specific testicular cells
To identify the cellular origin of the clusters, we compared the expression patterns of the five clusters with the expression patterns of cell markers investigated in our previous study of the same testicular material [10]. We found that four of the clusters identified in this study independently had similar expression profiles as four cell-specific markers identified in our previous study. The four markers were in our previous study associated with spermatogonia (cluster 1 in the current study; n = 109), spermatocytes (cluster 2; n = 164), spermatids (cluster 3; n = 246) and Sertoli cells (cluster 5; n = 196). Since transcripts in the early spermatid cluster showed a decrease in expression with a trough already at pi day 24, we assigned the fifth novel cluster, with transcripts that showed a decrease with a trough in expression at pi day 27, to late spermatids (cluster 4; n = 273). The expression level of most The table lists the 39 previously published cell-specific markers used to identify the cellular origin of the transcripts clusters [16][17][18].
transcripts returned to pre-irradiation levels around pi day 40 where a complete round of spermatogenesis had taken place [21].
To add further evidence of the association of cell types to the clusters, we identified previously published cellspecific markers in the five clusters [16][17][18]. In total, 39 markers were identified in the five clusters representing spermatogonia (n = 1), spermatocytes (n = 19), spermatids (n = 11), Leydig cells (n = 3), PTM cells (n = 1), and Sertoli cells (n = 4) (Table 1, Figure 1 and Figure 2). The single marker associated with spermatogonia was found in the corresponding cluster (cluster 1). The spermatocyte markers were evenly distributed in all five clusters, but comparing the expression patterns to the profiles identified earlier [10,16] strongly indicated that this cluster (cluster 2) with a reduced expression at pi day 17, mainly contained transcripts from spermatocytes. Most spermatid markers were present in the two clusters that we associated with early and late spermatids (cluster 3 and 4). The transcripts in the somatic cell cluster contained all somatic cell markers (cluster 5), which confirmed that the cluster represented all somatic cells and not only transcripts from Sertoli cells as suggested in our earlier study [10]. Gene set enrichment analysis GSEA was performed on the genes from each cluster separately to gain knowledge of overrepresented functional categories in each cluster. Statistical significant categories are listed in Table 2. Of particular interest, this analysis revealed overrepresentation of genes involved in mRNA processing, spliceosome and nucleus for the transcripts in the spermatogonia cluster. The transcripts in the spermatocyte cluster were enriched in male gamete generation, spermatogenesis and cell cycle process. The early spermatid cluster differed from the late spermatid cluster in enrichment for acrosomal vesicle. The somatic cell cluster was enriched in processes such as lipid metabolic process, response to oxidative stress and lipid biosynthetic process. All enriched Gene Ontology categories are available in the additional material (see Additional file 7: Table S6, Additional file 8: Table S7, Additional file 9: Table  S8, Additional file 10: Table S9 and Additional file 11: Table S10).

IHC suggested irradiation-induced Leydig cell hyperplasia
Large expression changes were observed in the somatic cell cluster. We expected this to be caused mainly by changes in cellularity, but we decided to use IHC to investigate whether histological changes contributed to the observed changes in expression. We stained sections of the irradiated testis tissue with the following somatic cell-specific markers: two Leydig cell markers: Tgfbr3 and Hsd3b; a PTM cell marker: SMA; and a Sertoli cell marker: Vimentin (Figure 3).
In control material (day 0; Figure 3), Tgfbr3 was expressed in the cell membrane of the majority of Leydig cells. Already at pi day 7 the Leydig cell membrane reaction appeared disturbed, incomplete or became cytoplasmic. In addition, the shape of the Leydig cells became more round and the close configuration apparent in the control material seemed to be lost. At pi day 21 some membrane reaction reappeared, but was lost again day 28-31 where hyperplasia of the Leydig cells was pronounced. During pi day 49-59 some of the Leydig cell membrane reactions was restored, but Tgfbr3 was only expressed at the membrane in about half of the Leydig cells at pi day 59 while the rest either showed a cytoplasmic reaction or no reaction at all.
Hsd3b was expressed in the Leydig cell cytoplasm in the control (day 0; Figure 3). Only a few Leydig cells reacted strongly, whereas the majority reacted moderately to faint or were negative. Again, already from pi day 7, this picture changed with the vast majority of the Leydig cells now reacting strongly to the antibody and although the intensity of the reaction gradually declined, it was still elevated at pi day 59. Staining with Hsd3b indicated an increase in Leydig cell numbers at pi day 21-28, but the number gradually returned to almost pre-irradiation levels at pi day 59.
Vimentin showed a specific staining in Sertoli cells in the control material (day 0; Figure 3). The staining remained in Sertoli cells throughout recovery, but the intensity and amount of staining seemed to decrease as specific germ cells disappeared. Already from pi day 7 where only  spermatogonia are absent from the testis, the amount of staining was reduced, and it remained reduced until pi day 49 where the intensity and amount of staining was back to control levels. In addition, the classical extension of Sertoli cell cytoplasm towards the centre of the tubuli was mainly observed at day 0 and after pi day 49. Although the staining was reduced as distinct germ cell populations were missing, we did not find any indications of cell death among the Sertoli cells (deduced from the Vimentin staining). Neither did we find any histological changes pi in the PTM cell population by staining with SMA (results not shown). Nevertheless, IHC staining of the testis tissue pi with the Leydig cell markers Tgfbr3 and Hsd3b suggested that some Leydig cell hyperplasia occurred after low dose irradiation. To confirm the hyperplasia pi, we quantified the ratio of the areas occupied by Leydig cells relative to the areas covered by other testicular cells. We found that the area of Leydig cells increased after irradiation and gradually decreased towards pre-irradiation levels as the testis was repopulated by germ cells (Figure 4). This further underlined that the number of Leydig cells seemed to be affected by low dose irradiation that may have caused Leydig cell hyperplasia.

Cell-specific transcript clusters
We identified five unique transcript clusters with an unsupervised cluster analysis of the transcripts that changed in expression during recovery from irradiation. The expression changes observed in this study were caused by a combination of changes in cellularity and cells changing transcription pattern due to their differentiation towards spermatids and spermatozoa -both events change the total RNA pool. It further complicates the separation of the clusters that cells have a mixed cellular expression of many of the genes that are specific for spermatogenesis e.g. a relatively low expression that initiates in spermatogonia that is followed by a large increase in expression in spermatocytes (for example Deleted in azoospermia-like (Dazl) [22]). We identified two spermatid clusters, one associated with early spermatids and a novel cluster associated with late spermatids. We anticipate that the transcripts in the early spermatid cluster initially were expressed in late pachytene or diplotene spermatocytes just before they entered meiosis, but with their main expression in early spermatids. It has previously been hypothesized that DNA is inaccessible to transcription during meiosis but also that a large group of mRNAs needed in early spermatids are expressed in the late spermatocyte stages [16,23]. Transcription is reinitiated in early spermatids [16] resulting in a massive wave of transcriptional activity immediately after meiosis before the chromatin is condensed and transcription is silenced during spermiogenesis [24,25]. Thus, the late spermatid cluster that we identified likely includes transcripts that are initially expressed during this wave of transcription in early spermatids.
The somatic cell cluster included markers specific for Leydig, PTM and Sertoli cells, suggesting that this cluster contained transcripts expressed by all the somatic cells in the testis. The somatic cells further comprise the microvasculature, lymph vessels, nerve fibres and connective tissue [26] and transcripts expressed by these cells are properly also present in the somatic cell cluster. From inspection of the expression profile of the somatic cell cluster it seemed that the somatic cells underwent large expression changes pi. However, the majority of the somatic cells did not seem to be affected by the irradiation and the gene expression changes detected in the somatic cell cluster were most likely caused by changes in cellularity. Since pachytene spermatocytes physically are the largest and most RNA-rich germ cells in the testis [10,26] their disappearing leads to the apparent expressional upregulation of genes from the somatic cell cluster where especially RNA from Sertoli cells comprise a larger part of the RNA from the testis.
We tried to separate more than five cell clusters from this data set. Yet, some of the original clusters just separated into subclusters with very similar expression patterns. The matured spermatozoa have a highly condensed transcriptionally inactive chromatin [27] and it is therefore reasonable that we did not identify any cluster(s) originating from the later germ cell stages.

Gene set enrichment analysis
We performed GSEA of the subsets of the transcripts in each cluster to gain knowledge of their protein function. The analyses supported the association of the clusters with the cell types. However, the data set was complex and other cells than those specifically associated to the cluster may contribute to the pool of transcripts in each cluster.
The transcripts in the spermatogonia cluster were enriched for transcripts involved in mRNA processing, spliceosome and nucleus, in agreement with the transcriptional regulation involved during the numerous divisions and differentiations undergone by spermatogonia [28]. The transcripts in the spermatocyte cluster were enriched in male gamete generation, spermatogenesis and cell cycle process, which fit well with the meiotic processes in spermatocytes. Both spermatid clusters were enriched in spermatogenesis. The early spermatid cluster was additionally enriched in acrosomal vesicle, which agrees with the acrosome being formed during the maturation of early to late spermatids and further to spermatozoa [29]. The transcripts in the late spermatid cluster were in addition enriched in serine-type peptidase activity and especially one serine protease, acrosin, has been identified as an important enzyme in the acrosome [30].
The very diverse functions of the somatic cells were also observed in the GSEA where the enriched categories among others included lipid metabolic process, response to oxidative stress and lipid biosynthetic process. Since cholesterol is the basis for steroidogenesis in Leydig cells, lipid metabolism is essential for steroid production [31][32][33]. Both the many cell divisions in the germinal tissue and Leydig cell steroidogenesis require large amounts Figure 4 In silico determination of Leydig cells relative to other testicular cells pi. Areas covered by Leydig cells highly increased following irradiation compared to areas covered by other testicular cells. The peak is observed at pi day 14 where the ratio decreased again. of oxygen. The testis is poorly vascularized and therefore vulnerable to oxidative stress and consequently have several defence mechanisms where Sertoli cells superoxide dismutase is important [34]. Thus, the results from the GSEA on the transcripts in the somatic cell cluster represented both Leydig and Sertoli cell functions.
Low dose irradiation may cause Leydig cell hyperplasia IHC validation of the two Leydig cell markers, Tgfbr3 and Hsd3b, showed that low dose irradiation affected the appearance of the Leydig cells with indications of Leydig cell hyperplasia around pi day 21-31. Previously published studies reported that human Leydig cells are preserved even after high dose irradiation up to 20 Gy, but that the endocrine function of the cells are affected also at lower doses [35][36][37][38]. In this study, the testes were irradiated with 1 Gy in a single dose, which is a low dose even in mice compared to the clinical doses of 16-20 Gy used to treat human testicular cancer [39].
Leydig cell hyperplasia has also been observed in men with fertility problems and testicular cancer [40]. The proliferation might arise due to a decreased testosterone production by the Leydig cells pi, which will induce increased LH production that may stimulate proliferation of pre-Leydig cells and their differentiation to Leydig cells. This is also the clinical observation upon irradiation of human testis harbouring carcinoma in situ [41] and in rats with Sertoli cell-only testes [42]. However, it is not clear if a dose of 1 Gy leads to reduced testosterone production and increased LH level, although higher irradiation doses lead to changes in testosterone and LH levels [43]. Also, the presence of Follicle-stimulating hormone (FSH) receptors in Leydig cells could suggest that the proliferation could be a response to decreased inhibin production caused by the absence of germ cells, which leads to increased FSH levels that may induce proliferation of pre-Leydig cells [44]. Further, the GSEA showed that the transcripts in the somatic cell cluster were highly enriched in lipid and steroid biosynthetic processes (see Additional file 11: Table S10). This enrichment can either be caused by an increased activity of the remaining Leydig cells or by an increased number of Leydig cells. From our results it seems to be caused by an increased number of Leydig cells. Alternatively, the reduced number of germ cells in the tubuli will likely result in a reduced size of the tubuli, and this might alter the ratio between tubuli and interstitial areas.
IHC staining of Vimentin specific to Sertoli cells did not show any histologically changes during recovery from low dose irradiation. Earlier studies found Sertoli cell death following irradiation, but these were all studies that used higher irradiation doses [45,46]. SMA staining of PTM cells did not suggest any histological changes in the PTM cells.

Conclusion
In this study, we investigated the molecular and cellular changes following re-establishment of spermatogenesis pi by analysis of the testicular transcriptome. We identified five clusters each representing a profile of a specific type of cell(s) during recovery from irradiation. We associated the five clusters with: spermatogonia, spermatocytes, early spermatids, late spermatids and somatic cells, respectively. We observed large expression changes in the somatic cell cluster during the time series. Therefore, we investigated if changes in histology, in addition to changes in cellularity, contributed to this observation. IHC staining of Leydig cell-specific markers suggested Leydig cell hyperplasia at pi day 21-31, which was further supported by in silico quantification of the relative area occupied by Leydig cells.