- Research article
- Open Access
Functional changes in mRNA expression and alternative pre-mRNA splicing associated with the effects of nutrition on apoptosis and spermatogenesis in the adult testis
BMC Genomics volume 18, Article number: 64 (2017)
The effects of nutrition on testis mass in the sexually mature male have long been known, however, the cellular and molecular processes of the testis response to nutrition was not fully understood.
We tested whether the defects in spermatogenesis and increases in germ cell apoptosis in the testis that are induced by under-nutrition are associated with changes in mRNA expression and pre-mRNA alternative splicing using groups of 8 male sheep fed for a 10% increase or 10% decrease in body mass over 65 days.
We identified 2,243 mRNAs, including TP53 and Claudin 11, that were differentially expressed in testis from underfed and well-fed sheep (FDR < 0.1), and found that their expression changed in parallel with variations in germ cell numbers, testis size, and spermatogenesis. Furthermore, pairs of 269 mRNAs and 48 miRNAs were identified on the basis of target prediction. The regulatory effect of miRNAs on mRNA expression, in combination with functional analysis, suggests that these miRNAs are involved in abnormal reproductive morphology, apoptosis and male infertility. Nutrition did not affect the total number of alternative splicing events, but affected 206 alternative splicing events. A total of 159 genes, including CREM, SPATA6, and DDX4, were differentially spliced between dietary treatments, with functions related to RNA splicing and spermatogenesis. In addition, three gene modules were positively correlated with spermatogenesis-related phenotypic traits and negatively related to apoptosis-related phenotypic traits. Among these gene modules, seven (CFLAR, PTPRC, F2R, MAP3K1, EPHA7, APP, BCAP31) were also differentially expressed between nutritional treatments, indicating their potential as markers of spermatogenesis or apoptosis.
Our findings on significant changes in mRNAs and pre-mRNA alternative splicing under-nutrition suggest that they may partly explain the disruption of spermatogenesis and the increase germ cell apoptosis. However, more research is required to verify their causal effects in regulating spermatogenesis and germ cell apoptosis.
The development of mature haploid spermatozoa from diploid spermatogonial cells  can be affected by many factors, including photoperiod, hormones, temperature and nutrition. The effects of nutrition on testis mass in the sexually mature male have long been known, as has the direct relationship between testicular mass and sperm production . In addition, with change in testicular size, the efficiency of sperm production also changes . We have been investigating the cellular and molecular processes of the testis response to nutrition, and we have found that under-nutrition despaired spermatogenesis in adult sheep [4, 5].
Within the testis, spermatogenesis is a strictly regulated process, at both the transcriptional and the post-transcriptional level . In recent years, a novel mechanism of post-transcriptional control, mediated by microRNAs (miRNAs), has emerged as an important regulator of spermatogenesis [6, 7]. miRNAs are small (~22 nucleotides) endogenous RNAs that negatively regulate gene expression by targeting the 3′untranslated region (3′UTR)  and/or coding region  of mRNAs. We have recently found that the expression of a number of miRNAs is affected by nutrition in sexually mature male sheep, and most of the predicted targets of the differentially expressed miRNAs were mainly involved in reproductive system development and function . However, the regulatory relationship between these miRNAs and their corresponding mRNAs targets in testis remains to be verified. We therefore decided to profile mRNA expression in the testes of well-fed and underfed male sheep using RNA-seq so we could explore the relationships between the miRNAs we had identified and their putative targets.
In addition to the disruption of spermatogenesis, under-nutrition of sexually mature male sheep increased apoptosis in germ cells . Our recent study has revealed higher levels of expression of miR-98 in underfed sheep than in well-fed sheep  and this miRNA has been reported to play a critical role in apoptosis . Since the molecular mechanisms through which miRNAs regulate the expression of apoptosis-related genes are still controversial, we decided to explore these processes further using our nutritional model.
It has also been reported that spermatogenesis and a large number of apoptotic factors are regulated by alternative pre-mRNA splicing that generates multiple transcript species from a common mRNA precursor and thus raises protein diversity and allows the system to cope with the increasingly broad spectrum of functional and behavioural complexity [12, 13]. To date, eight types of alternative splicing have been reported: cassette exon, alternative 5′ splice site, alternative 3′ splice site, mutually exclusive exon, coordinates cassette exons, alternative first exon, alternative last exon and intron retention . We therefore also tested the hypothesis that nutritional treatment will induce differences in alternative splicing, and these changes will be related to the regulation of spermatogenesis and germ cell apoptosis in the testis.
Overall, this study used testicular tissue from under-fed and well-fed sexually mature sheep to pursue four objectives: 1) To investigate the differences of the expression of mRNAs; 2) To evaluate the influence of miRNAs on spermatogenesis and the expression of apoptosis-related genes; 3) To explore the relationships between alternative pre-mRNA splicing and spermatogenesis and apoptosis; 4) To investigate the relationships between the gene modules and spermatogenesis and apoptosis related phenotypic traits.
The experimental protocol was approved by the Animal Ethics Committee of the CSIRO Centre for Environment and Life Sciences, Floreat, Western Australia (Project No.1202).
Animals and treatments
From May to July (autumn-winter), 24 Merino male sheep (age 24 months, body mass 65.7 ± 4.7 kg, scrotal circumference 31.8 ± 2.5 cm) were housed in individual pens in a building with windows that allowed good penetration of natural light at Floreat, Western Australia. During the 3-week acclimatization period, all sheep were fed daily with 750 g oaten chaff (8.4% crude protein; 8.0 MJ/Kg metabolisable energy) and 200 g lupin grain (35.8% crude protein; 13.0 MJ/Kg metabolisable energy). At the start of the treatment period (end of May; mid-autumn), the animals were allocated into three dietary treatment groups (high, maintenance and low) balanced for training success to semen collection, body mass, scrotal circumference, temperament, poll-horn type, and sperm quality (the percentage of live and motile sperm, sperm concentration). Animals fed their maintenance requirements were expected to maintain constant body mass. The high diet was designed to allow the animals to gain 10% live weight over 65 days whereas the low diet was designed to allow 10% loss in weight. At the start of the treatment period, individual daily allowance was 1.2 kg oaten chaff plus 0.3 kg lupin grain for the rams in the high-diet group, 0.7 kg chaff and 0.18 kg lupin grain for the maintenance group, and 0.51 kg chaff and 0.13 kg lupin grain for the low-diet group. Every week, the animals were weighed and the amount of feed offered to each individual was adjusted to ensure achievement of target live weight.
The data on the effects of nutritional treatments on body mass, scrotal circumference, semen quality and spermatozoal quality, Sertoli cell number and quality were published before [4, 5]. Our previous work showed significant differences between High diet and Low diet in terms of the above parameters, with the values for maintenance-fed rams being generally similar to those for the High diet. For purposes of the current study, we decided to perform RNA-Seq only with the two extreme groups (testis growing versus testis shrinking) to determine the factors that contribute to the huge differences in phenotype.
Tissue Collection and preservation
After 65 days, all male sheep were killed with intravenous overdose of sodium pentobarbitone, and the testes were immediately removed, dissected and weighed. Three samples were chosen from top, middle and bottom parts of both testes (~1 cm3 for each); those from the right testis were snap-frozen in liquid nitrogen and stored at −80 °C for RNA isolation.
Isolation of RNA
About 1 cm3 tissues from top, middle and bottom parts of the right testes were mixed and grinded to powder in liquid nitrogen for RNA isolation. The trizol method was used to isolate total RNA  from testis samples. The quality and quantity of the RNA were determined by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA) and Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA). Only RNA with an integrity number (RIN) greater than 7.0 was used for further analysis.
Small RNA library sequencing and Identification of miRNAs
The methodology details and outcomes for miRNAs have been reported elsewhere .
Construction and sequencing of the RNA-seq library
Total RNA (1.0 μg each) from each sample was used to construct RNA-seq libraries with a unique index using the TruSeq mRNA Sample Preparation kit (Illumina, San Diego, CA) according to the manufacturer’s instruction. Quantitative real time PCR (qPCR) was performed for library quantification using the StepOnePlus™ Real-Time PCR System (Applied Biosystems, Carlsbad, CA) and KAPA SYBR Fast ABI Prism qPCR kit (Kapa Biosystems, Woburn, MA). Individual libraries were then pooled for sequencing at Génome Québec (Montréal, Canada) using the HiSeq 2000 system (Illumina). Sequencing was performed as 100 bp paired-end reads. All reads were de-multiplexed according to their index sequences with CASAVA version 1.8 (Illumina) and reads that did not pass the Illumina chastity filter were discarded.
Mapping and annotation of RNA-seq reads
RNA-seq reads were aligned to the ovine genome (OAR 3.1) using Tophat 2.0.10 with default parameters . Each BAM file obtained from TopHat2 and the GTF file obtained from ENSEMBL (v75.30) were used in the htseq-count (http://www-huber.embl.de/users/anders/HTSeq/) to determine the number of reads mapped to each gene.
Identification of differentially expressed (DE) mRNAs
Differentially expressed (DE) mRNAs were investigated with the bioinformatics tool, edgeR that utilizes a negative binomial distribution to model sequencing data . The expression of mRNAs in each library was normalized to counts per million reads (CPM) with the formula: CPM = (reads number/total reads number per library) × 1,000,000. mRNAs with CPM > 5 in at least 50% of the samples were subjected to DE analysis. Fold changes were defined as ratios of arithmetic means of CPM within each comparison group. Significant DE mRNAs were determined by an adjusted P value (False discovery rate, FDR < 0.1) based on Benjamini and Hochberg multiple testing correction  as well as fold change > 1.5 .
Validation of mRNA expression using RT-qPCR
RT-qPCR was performed using SYBR Green (Fast SYBR® Green Master Mix; Applied Biosystems) to validate mRNA expression of six differentially expressed genes: PIWIL1, SPATA4, INHBA, FOXO3, PTEN, CYP51A1. Oligonucleotide primer sequences for these genes were designed using NCBI primer blast (http://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome) and the primer for glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was obtained from a published source (Table 1). Total RNA (1 μg) from each sample was treated with DNAase I (Invitrogen), and reverse-transcribed to cDNA using SuperScript II reverse transcriptase following the manufacturer’s protocol (Invitrogen). Fluorescence signal was detected with StepOnePlus™ Real-Time PCR System (Applied Biosystems). In total, each reaction contained 10 μl Fast SYBR Green Master Mix (Applied Biosystems), 1 μl of forward primer (20 pmol/μl), 1 μl of reverse primer (20 pmol/μl), 7 μl nuclease-free water, and 1 μl cDNA template. Samples were measured in triplicate using the following protocol: 95 °C for 10 min for initial denaturation and then 40 cycles of 95 °C for 20 s, followed by annealing/extension for 30 s at 60 °C. Analysis of melting curves was used to monitor PCR product purity. Amplification of a single PCR product was confirmed by agarose gel electrophoresis and DNA sequencing (data not shown). One-way ANOVA was used to compare the groups, and P < 0.05 was considered significant. Data are expressed as Mean ± SEM.
Construction of miRNA-mRNA regulatory relationships
The results for miRNAs were all obtained from a previous study using the same samples . The predicted regulatory relationships between differentially expressed miRNAs and differentially expressed mRNAs were identified on the basis of two criteria as suggested by previous studies [20, 21]: negative correlation and computational target prediction. Genes targeted by miRNAs were predicted by TargetScan Release 6.0 (http://www.targetscan.org/) and miRanda (http://www.microrna.org/microrna/home.do). Target genes predicted by both TargetScan (default parameters) and miRanda (total score > 145, total energy < −10 kcal/mol) were used for further analysis . Pairwise Pearson correlation coefficient (R) was computed for each miRNA and their predicted target genes, and multiple testing corrections were done by calculating FDR. Significant miRNA-mRNA pairs were defined as R < −0.9 and FDR < 0.1.
Functional analysis for DE mRNAs
The identities of the DE mRNAs in the miRNA-mRNA regulatory relationships were uploaded into IPA software (Ingenuity Systems, www.ingenuity.com) to detect the top functions. A threshold of P < 0.01 was applied to enrich significant biological functions. The IPA regulation z-score algorithm was used to predict the direction of change for a given function (increase or decrease), with a z-score > 2 suggesting a significant increase whereas a z-score < −2 suggesting a significant decrease. The GO terms were defined and the KEGG pathways were enriched using Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov) . For each analysis, the functional annotation clustering option was used and significant GO terms and KEGG pathways were declared at P < 0.05 and molecule number > 2.
Identification and annotation of alternative splicing (AS) events
TopHat2 was used to predict the splice junctions (options: −a 19, −g 1, −-max-intron-length 17325, and --min-intron-length 81). A total of 15 million reads were randomly selected from each sample for analysis to make sure that the comparison was at the same level . JuncBASE was used to annotate all AS events (cassette exons, alternative 5′ splice site, alternative 3′ splice site, mutually exclusive exons, coordinate cassette exons, alternative first exons, alternative last exons, and intron retention) and calculate the Percentage Spliced Index (PSI) . Splicing analysis was performed for events that had at least 20 reads and the PSI differences (ΔPSI) are higher than 10% .
The relationship between gene modules with sheep phenotypic traits
A weighted gene co-expression network was constructed for all samples using the WGCNA package in R to analyse the expressed mRNAs (CPM > 1 in at least 8 samples) . Briefly, a matrix of pairwise correlations between all pairs of genes across all samples was constructed. An adjacency matrix was then calculated, using the correlation matrix of the expression sets, and transformed into a topological overlap matrix that was then used to derive a distance matrix of hierarchical clustering. Finally, the mRNAs were assigned into different modules based on hierarchical clustering. Modules with eigengenes (defined as the first principal component of each module and considered as a representative of the gene expression profiles in that module) that were highly correlated were then merged. All these steps were performed using the ‘blockwiseModules’ function in the WGCNA package, using the major parameters described previously . The correlations were calculated for the relationships between module eigengenes and phenotypes (testis weight, sperm number per testis, tubule diameter, seminiferous epithelium volume, change of scrotal circumference, apoptotic cells/tubule; Additional file 1: Table S1).
High Quality RNA-seq data were obtained from all samples
More than 350 million sequenced paired-end reads were obtained from 16 libraries, of which an average of 76% could be mapped to OAR3.1 (http://www.livestockgenomics.csiro.au/). The genomics region of reads, the RNA-seq 3′/5′ bias and the sequencing depth were analysed to evaluate the quality of the RNA-seq data. Around 50% of the reads were derived from exonic regions, while around 20% were derived from intergenic or intronic regions (Fig. 1a and b). In general, the coverage of reads along each transcript revealed no obvious 3′/5′ bias (Fig. 1c). As can be seen in Fig. 1d, the number of transcripts detected increased as the number of the sequencing reads increased. Finally, analysis of the sequencing coverage along each chromosome showed extensive transcriptional activity for the entire genome (Additional file 2: Figure S1).
Profile of mRNAs in sheep testis
An average of 13,980,416 (SD = 2,930,788) reads from high diet and 11,014,809 (SD = 2,524,631) from low diet were mapped to Ensembl gene annotation database (P < 0.05). A total of 13,859 genes were detected in testicular tissue from the low diet group, compared to 14,561 from the high diet group. In total, 11,748 genes were expressed in all 16 animals. The most abundant transcript (~2% of total reads) was from the 7SK gene, a small nuclear RNA involved in pre-mRNA splicing and processing. Functional analysis with DAVID software revealed that the most highly expressed 3,000 genes were mainly related to cell cycles, protein catabolic processes, and spermatogenesis (Table 2). Only genes (14,385) that were expressed in at least 8 libraries were used for further analysis.
Effects of nutrition on mRNA expression
In total, 2,243 mRNAs were found to be differentially expressed (DE) when comparing underfed with well-fed male sheep (Additional file 3: Table S2), of which 1,081 were expressed more in underfed than well-fed sheep (eg, TP53 and Claudin 11) and 1,162 were expressed less in underfed than well-fed sheep (eg, CYP51A1 and SPATA4). IPA analysis revealed that the functions of most of the DE mRNAs are related to quantity of germ cells, testis size, quantity of Sertoli cells, and quantity of connective tissue cells (Fig. 2). We considered genes that were related to more than one function may be more important, such as FOXO3, PTEN, CYP51A1, INHBA and SPATA4. Therefore, they were selected for further RT-qPCR validation analysis. Further functional analysis using DAVID, produced largely the same outcome, indicating that most common functions of DE mRNAs to be in the cell cycle (n = 116), spermatid development (n = 16), spermatogenesis (n = 48), and DNA replication (n = 18) (Additional file 4: Table S3). Importantly, one gene, PIWIL1 (MIWI), were associated with all of these functions.
RT-qPCR validation of differentially expressed genes
Six mRNAs were selected from the DE mRNAs and the RT-qPCR results were consistent with the sequencing data for all of them. For example, both the sequencing data and the RT-qPCR results showed that PIWIL1 was expressed at a lower level in underfed males than in well-fed males (Additional file 5: Figure S2). In addition, the expression of INHBA was down regulated in well-fed male sheep (Additional file 5: Figure S2).
The relationship between gene modules with sheep phenotypic traits
A total of 15 modules were obtained using WGCNA analysis, of which three (Modules 7, 9 and 10) were of interest because the relationships were strong (correlation coefficient > ± 0.5; P < 0.01). These modules were negatively correlated with testis weight, tubule diameter, seminiferous epithelium volume and change in scrotal circumference, and positively correlated with apoptotic cells/tubule. Among the 15 modules, only Module 11 was correlated with sperm number per testis (r = 0.56, P = 0.02; Fig. 3).
Functional analysis for the genes in Modules 7, 9 and 10
Functional analysis suggested that genes in Modules 7 and 9 were associated with both spermatogenesis and apoptosis, but genes in Module 10 were only related to apoptosis. Specifically, in Module 7, 35 genes were related to spermatogenesis whereas 46 genes were associated with apoptosis. Six genes (NFKBIL1, XRCC5, ERCC1, APP, BCAP31, RRAGA) were related to both spermatogenesis and apoptosis. For Module 9, eight genes (RXFP1, ITCH, ITGB1, XHD, TAF7L, WNT2, SPIN1, LNPEP) were associated with spermatogenesis, whereas 16 genes (CFLAR, TAF9B, CCK, VAV3, NR3C1, CDH13, CROP, CDKN1B, PTPRC, ATP7A, RTN3, HSPD1, ITM2B, F2R, MAP3K1, RAD21) were associated with apoptosis, but no genes were common to spermatogenesis and apoptosis. For Module 10, four genes (EPHA7, SCIN, NGFRAP1 and MAP3K7) were related to apoptosis. Interestingly, of all the genes mentioned above, seven (APP, BCAP31, CFLAR, PTPRC, F2R, MAP3K1and EPHA7) were differentially expressed between nutritional treatments, indicating their pivotal roles in reproduction or apoptosis.
miRNA-mRNA regulatory relationships
Putative miRNA-mRNA pairs were identified on the basis of target prediction and the negative regulatory effect of miRNAs on the expression levels of their target genes. A total of 940 miRNA-mRNA pairs (48 miRNAs, 269 mRNAs) were identified (Additional file 6: Table S4). Among these pairs, oar-novel-miR-33 and oar-novel-miR-31 paired with the highest number of mRNAs: oar-novel-miR-33 paired with 68 mRNAs and oar-novel-miR-31 paired with 52 mRNAs (Additional file 7: Figure S3). IPA analysis indicated that the mRNAs in the negative pairs were mainly involved with organization of cytoplasm, cell morphology, abnormal morphology of the reproductive system, cell death and male infertility (Fig. 4a). In addition, these mRNAs were also involved in 76 signalling pathways, of which Sertoli cell-Sertoli cell junction signalling, germ cell-Sertoli cell junction signalling, androgen signalling, and apoptosis signalling were among the 15 most relevant (Fig. 4b). FOXO3 and PTEN were related to more than 8 functions out of 15 most common functions, we assumed these two genes may be crucial for testes function, therefore, they were selected for RT-qPCR validation. The expression of FOXO3 was higher in well-fed sheep than underfed sheep, while PTEN expression was lower in well-fed than underfed sheep.
Identification of alternative splicing events
We initially obtained 940,607 junctions from the 16 RNA-seq libraries with the Tophat2 software. Totally, 6376 alternative splicing events (∆PSI > 10% in at least one library and at least 20 reads mapped) were identified from these junctions (Additional file 8: Table S5). Among these alternative splicing events, 4,820 (75.6%) were previously annotated as known alternative splicing events in the Ensembl Database, which can be mapped to 2288 unique genes. Eight types of alternative splicing events were identified, including 1,131 cassette exons, 645 alternative 5′ splice sites, 578 alternative 3′ splice sites, 17 mutually exclusive exons, 86 coordinate cassette exons, 578 alternative first exons, 247 alternative last exons, and 2,810 intron retentions.
Effects of nutrition on alternative splicing
We found 2551 ± 189 (Mean ± SEM) alternative splicing events in the High Diet group and 2455 ± 126 alternative splicing events in the Low Diet group (not significant). There was no difference as for the total number of each type of alternative splicing event between two groups (Fig. 5). PSI values from each diet group were used to test the differentially spliced genes between treatments, resulting 206 differentially spliced isoforms (21 alternative 3′ splice sites, 23 alternative 5′ splice sites, 34 alternative first exon, 17 alternative last exon, 86 cassette, 7 coordinate cassette exon, 8 intron retention events, P < 0.05, Wilcoxon test). These differentially spliced isoforms were mapped to 159 known unique genes, and DAVID functional analysis revealed the most common functions of these genes were related to RNA splicing and spermatogenesis (Fig. 6).
This study appears to be the first to profile the whole transcriptome in sheep testis, to construct the regulatory relationships between miRNAs and mRNAs in testis, to explore the relationships between pre-mRNA alternative splicing and testis function, and to link the gene modules with phenotypic traits related to spermatogenesis and apoptosis. In the context of an experimental model of reversible testis growth in the sexually mature male, we have been able to identify mRNAs that are associated with testis function and, more importantly, apoptosis in germ cells. These findings strongly support the hypothesis that the decline in spermatogenesis and increase in germ cell apoptosis induced by under-nutrition in the sexually mature male sheep are, at least, partially due to changes in mRNA expression and pre-mRNA alternative splicing.
We found over 2,000 mRNAs that were differentially expressed between treatments, with over 1,000 mRNAs, including TP53 and Claudin 11, that were more highly expressed in underfed than in well-fed sheep. This result supports our previous observations based on qPCR . A high level of TP53 indicates more cells going through apoptosis , and this result is in agreement with our previous finding-more TUNEL positive cells were observed in testes from underfed than in well-fed sheep , so we conclude that under-nutrition increases apoptosis in germ cells. On the other hand, Claudin-11 is a tight junction protein expressed in Sertoli cells and rarely in other cell types in the testis  and plays a central role in the formation of tight junctions [30, 31]. In testicular tissue from underfed sheep, increased expression of Claudin 11 and disorganization of Claudin 11 protein strongly indicate the impairment of tight junctions . In addition, in the present study, over 1000 mRNAs showed lower expression in underfed than in well-fed sheep, including CYP51A1 and SPATA4. CYP51A1 is a member of the cytochrome P450 family and is expressed strongly by germ cells , illustrating its crucial role in spermatogenesis. Therefore, the lower level of CYP51A1 expression in underfed sheep is coherent with the decrease in numbers of germ cells and defective spermatogenesis caused by undernutrition . SPATA4 has also been reported to be testis-specific and associated with spermatogenesis, and involved in maintaining spermatogenesis . Therefore, the reduced expression of SPATA4 in underfed sheep is also consistent with compromised spermatogenesis with under-nutrition.
The functional analysis using IPA revealed that, for mRNAs that are differentially expressed between nutritional treatments, the most common functions are quantity of germ cells, testis size, quantity of Sertoli cells and quantity of connective tissue cells. Thus, the differentially expressed transcriptomes are consistent with the reductions of testis mass and sperm production in underfed rams . Among all the genes that were related to these functions, we considered the genes that were related to more than one function may be more important, such as IGF1R, INHBA, TP53. In the testes of adult mice lacking IGF1R in their Sertoli cells, there is a reduction in testis size and daily sperm production , indicating a role for IGF1R in control of sperm production by Sertoli cells. A protein product of the INHBA gene, activin A, is an important regulator of testicular cell proliferation . As indicated above, TP53 regulates spermatogenesis by inducing apoptosis . We propose that, the expression of IGF1R, INHBA and TP53 may be used as biomarkers of sperm production.
To further investigate the crucial genes in spermatogenesis and apoptosis, we looked at the relationships between these phenotypes and the gene modules that we discovered in the testis. We found five genes (CFLAR, PTPRC, F2R, MAP3K1, EPHA7) that appear to be crucial for apoptosis and two genes (APP and BCAP31) related to both apoptosis and spermatogenesis. More importantly, all seven of these genes were differentially expressed between nutritional treatments, suggesting that they played pivotal roles in the control of testis function. These conclusions are supported by previous studies. For example, CFLAR is involved in inhibition of the death receptor-activated pathway ; MAP3K1 has both anti- and pro-apoptotic functions ; EPHRINA5-EPHA7 complex induces apoptosis through TNFR1 . Therefore, these seven genes are potential biomarkers of spermatogenesis and apoptosis. Interestingly, the changes in these genes were associated with change in testis mass, raising the possibility that factors associated with change in testicular tissue, rather than direct effects of nutritional treatments on testicular tissue, are responsible for changes in spermatogenesis and apoptosis [5, 10]. If this were to be the case, then the relationships between gene modules and phenotypes observed in the present study could be applied more generally to other factors that can cause changes in the testis mass, such as photoperiod, stress and temperament, or physical fitness.
In addition, this study further identified miRNA-RNA relationships which may regulate the above altered expression events. Of particular importance are oar-novel-miR-33 and oar-novel-miR-31 because they paired with the greatest number of mRNAs, indicating a crucial role in testis function. Novel-miR-33 is homologous to miR-296 that is specific to embryonic stem cells and has been reported to be highly conserved between species . So far, there is no direct evidence for a role for miR-296 in testis function. However, one study proved that miR-296 was more highly expressed in mature testis than in immature testis, indicating a pivotal role in spermatogenesis in the adult. In addition, miR-296 was also defined as anti-apoptotic . Therefore, the reduced expression of novel-miR-33 (miR-296) in underfed sheep  illustrates decreased testis function and increased cell apoptosis. By contrast, novel-miR-31 is homologous with miR-34 which has been shown to enhance germ cell phenotype during the late stages of spermatogenesis in other species . In the current study, therefore, the lower level of novel-miR-31 in underfed sheep is coherent with the loss of germ cell function .
Nutritional treatment did not affect the total number of alternative splicing junctions, in contrast with some previous reports of nutritional effects on other biological processes [43, 44]. The lack of effect of nutritional treatment on the total number of alternative splicing junctions suggests that alternative splicing is a fine-tuner in the testis that stabilizes testis function, as suggested for other tissues [45, 46]. However, with respect to specific genes, we found 159 that were differentially spliced between high diet and low diet groups. Functional analysis of these genes indicated roles in RNA splicing and spermatogenesis, and suggests that nutrition affects spermatogenesis by changing pre-mRNA splicing. For example, the alternative splicing event for CREM (cAMP response element modulator) is alternative last exon. It has been reported that CREM mRNA exhibits a remarkable array of alternative splice variants . For example, during male meiosis, the inactive CREM variant is switched to active CREM variant (by incorporation of transactivating domains) directed by alternative splicing. Therefore, in the mature male sheep, it is possible that nutrition affects the number of active CREM isoforms, potentially explaining the disruption of spermatogenesis in underfed sheep. Future studies on verification of alternative splicing events genes related to spermatogenesis detected by RNA-seq using qPCR and how such changes affect the activity, possibly involving construction of a shortened ‘minigene’ containing the regulated exons and splicing signals , will be required to better understand this process. In addition, it is essential to determine whether the changes in splicing can affect protein expression.
In conclusion, we have identified two molecular mechanisms that could explain the effect of nutrition on spermatogenesis and germ cell apoptosis in the adult male: 1) nutrition-induced changes in the expression of mRNAs in sheep testis, the functions of the differentially expressed mRNAs are mainly related to spermatogenesis and germ cell apoptosis, an important regulator in these processes are the regulatory relationships between miRNAs and mRNAs; 2) nutritional treatment causes differences in pre-mRNA alternative splicing, and these changes are closely involved in spermatogenesis and germ cell apoptosis in testis. Some differentially spliced genes (CREM and DDX4), and testis phenotype related genes (CFLAR, PTPRC, F2R, MAP3K1, EPHA7, APP and BCAP31) should be able to work as potential biomarkers for spermatogenesis and apoptosis. To make the study move forward, confirming the predicted functions of these genes using in vivo and in vitro experiments are required.
Hecht NB. Molecular mechanisms of male germ cell differentiation. Bioessays. 1998;20(7):555–61.
Oldham CM, Adams NR, Gherardi PB, Lindsay DR, Mackintosh JB. The influence of level of feed intake on sperm-producing capacity of testicular tissue in the ram. Aust J Agr Res. 1978;29(1):173–9.
Walkden-Brown SW, Restall BJ, Taylor WA. Testicular and epididymal sperm content in grazing Cashmere bucks: seasonal variation and prediction from measurements in vivo. Reprod Fertil Dev. 1994;6(6):727–36.
Guan Y, Malecki IA, Hawken PA, Linden MD, Martin GB. Under-nutrition reduces spermatogenic efficiency and sperm velocity, and increases sperm DNA damage in sexually mature male sheep. Anim Reprod Sci. 2014;149(3–4):163–72.
Guan Y, Liang G, Hawken PA, Meachem SJ, Malecki IA, Ham S, Stewart T, Guan le L, Martin GB. Nutrition affects Sertoli cell function but not Sertoli cell numbers in sexually mature male sheep. Reprod Fertil Dev. 2014;28(8):1152–63.
Papaioannou MD, Nef S. microRNAs in the testis: building up male fertility. J Androl. 2010;31(1):26–33.
Liang G, Malmuthuge N, le Guan L, Griebel P. Model systems to analyze the role of miRNAs and commensal microflora in bovine mucosal immune system development. Mol Immunol. 2015;66(1):57–67.
Krutzfeldt J, Stoffel M. MicroRNAs: a new class of regulatory genes affecting metabolism. Cell Metab. 2006;4(1):9–12.
Hausser J, Syed AP, Bilen B, Zavolan M. Analysis of CDS-located miRNA target sites suggests that they can effectively inhibit translation. Genome Res. 2013;23(4):604–15.
Guan Y, Liang G, Hawken PA, Malecki IA, Cozens G, Vercoe PE, Martin GB, Guan LL. Roles of small RNAs in the effects of nutrition on apoptosis and spermatogenesis in the adult testis. Sci Rep. 2015;5:10372.
Wang S, Tang Y, Cui H, Zhao X, Luo X, Pan W, Huang X, Shen N. Let-7/miR-98 regulate Fas and Fas-mediated apoptosis. Genes Immun. 2011;12(2):149–54.
Schwerk C, Schulze-Osthoff K. Regulation of apoptosis by alternative pre-mRNA splicing. Mol Cell. 2005;19(1):1–13.
Liang G, Malmuthuge N, Guan Y, Ren Y, Griebel PJ, Guan le L. Altered microRNA expression and pre-mRNA splicing events reveal new mechanisms associated with early stage Mycobacterium avium subspecies paratuberculosis infection. Sci Rep. 2016;6:24964.
Ding F, Cui P, Wang Z, Zhang S, Ali S, Xiong L. Genome-wide analysis of alternative splicing of pre-mRNA under salt stress in Arabidopsis. BMC Genomics. 2014;15:431.
Hellani A, Ji J, Mauduit C, Deschildre C, Tabone E, Benahmed M. Developmental and hormonal regulation of the expression of oligodendrocyte-specific protein/claudin 11 in mouse testis. Endocrinology. 2000;141(8):3012–9.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. Controlling the false discovery rate in behavior genetics research. Behav Brain Res. 2001;125(1–2):279–84.
McCarthy DJ, Smyth GK. Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics. 2009;25(6):765–71.
Lim EL, Trinh DL, Scott DW, Chu A, Krzywinski M, Zhao Y, Robertson AG, Mungall AJ, Schein J, Boyle M, et al. Comprehensive miRNA sequence analysis reveals survival differences in diffuse large B-cell lymphoma patients. Genome Biol. 2015;16:18.
Liang G, Malmuthuge N, Bao H, Stothard P, Griebel PJ, le Guan L. Transcriptome analysis reveals regional and temporal differences in mucosal immune system development in the small intestine of neonatal calves. BMC Genomics. 2016;17(1):602.
Liang G, Malmuthuge N, McFadden TB, Bao H, Griebel PJ, Stothard P, le Guan L. Potential regulatory role of microRNAs in the development of bovine gastrointestinal tract during early life. PLoS One. 2014;9(3), e92592.
da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.
Brooks AN, Yang L, Duff MO, Hansen KD, Park JW, Dudoit S, Brenner SE, Graveley BR. Conservation of an RNA regulatory map between Drosophila and mammals. Genome Res. 2011;21(2):193–202.
Wang ET, Sandberg R, Luo S, Khrebtukova I, Zhang L, Mayr C, Kingsmore SF, Schroth GP, Burge CB. Alternative isoform regulation in human tissue transcriptomes. Nature. 2008;456(7221):470–6.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Kommadath A, Bao H, Arantes AS, Plastow GS, Tuggle CK, Bearson SM, Guan le L, Stothard P. Gene co-expression network analysis identifies porcine genes associated with variation in Salmonella shedding. BMC Genomics. 2014;15:452.
Li H, Jogl G. Structural and biochemical studies of TIGAR (TP53-induced glycolysis and apoptosis regulator). J Biol Chem. 2009;284(3):1748–54.
Morita K, Sasaki H, Fujimoto K, Furuse M, Tsukita S. Claudin-11/OSP-based tight junctions of myelin sheaths in brain and Sertoli cells in testis. J Cell Biol. 1999;145(3):579–88.
Gow A, Southwood CM, Li JS, Pariali M, Riordan GP, Brodie SE, Danias J, Bronstein JM, Kachar B, Lazzarini RA. CNS myelin and sertoli cell tight junction strands are absent in Osp/claudin-11 null mice. Cell. 1999;99(6):649–59.
Wolburg H, Wolburg-Buchholz K, Liebner S, Engelhardt B. Claudin-1, claudin-2 and claudin-11 are present in tight junctions of choroid plexus epithelium of the mouse. Neurosci Lett. 2001;307(2):77–80.
Rozman D, Waterman MR. Lanosterol 14alpha-demethylase (CYP51) and spermatogenesis. Drug Metab Dispos. 1998;26(12):1199–201.
Liu SF, He S, Liu BW, Zhao Y, Wang Z. Cloning and characterization of testis-specific spermatogenesis associated gene homologous to human SPATA4 in rat. Biol Pharm Bull. 2004;27(11):1867–70.
Pitetti JL, Calvel P, Zimmermann C, Conne B, Papaioannou MD, Aubry F, Cederroth CR, Urner F, Fumel B, Crausaz M, et al. An essential role for insulin and IGF1 receptors in regulating sertoli cell proliferation, testis size, and FSH action in mice. Mol Endocrinol. 2013;27(5):814–27.
Archambeault DR, Yao HH. Activin A, a product of fetal Leydig cells, is a unique paracrine regulator of Sertoli cell proliferation and fetal testis cord expansion. Proc Natl Acad Sci U S A. 2010;107(23):10526–31.
Chen D, Zheng W, Lin A, Uyhazi K, Zhao H, Lin H. Pumilio 1 suppresses multiple activators of p53 to safeguard spermatogenesis. Curr Biol. 2012;22(5):420–5.
Irmler M, Thome M, Hahne M, Schneider P, Hofmann K, Steiner V, Bodmer JL, Schroter M, Burns K, Mattmann C, et al. Inhibition of death receptor signals by cellular FLIP. Nature. 1997;388(6638):190–5.
Pham TT, Angus SP, Johnson GL. MAP3K1: Genomic Alterations in Cancer and Function in Promoting Cell Survival or Apoptosis. Genes Cancer. 2013;4(11–12):419–26.
Lee H, Park E, Kim Y, Park S. EphrinA5-EphA7 complex induces apoptotic cell death via TNFR1. Mol Cells. 2013;35(5):450–5.
Gangaraju VK, Lin H. MicroRNAs: key regulators of stem cells. Nat Rev Mol Cell Biol. 2009;10(2):116–25.
Cheng AM, Byrom MW, Shelton J, Ford LP. Antisense inhibition of human miRNAs and indications for an involvement of miRNA in cell growth and apoptosis. Nucleic Acids Res. 2005;33(4):1290–7.
Bouhallier F, Allioli N, Lavial F, Chalmel F, Perrard MH, Durand P, Samarut J, Pain B, Rouault JP. Role of miR-34c microRNA in the late steps of spermatogenesis. RNA. 2010;16(4):720–31.
Salati LM, Szeszel-Fedorowicz W, Tao H, Gibson MA, Amir-Ahmady B, Stabile LP, Hodge DL. Nutritional regulation of mRNA processing. J Nutr. 2004;134(9):2437S–43S.
Kulseth MA, Berge KE, Bogsrud MP, Leren TP. Analysis of LDLR mRNA in patients with familial hypercholesterolemia revealed a novel mutation in intron 14, which activates a cryptic splice site. J Hum Genet. 2010;55(10):676–80.
Wu CT, Chiou CY, Chiu HC, Yang UC. Fine-tuning of microRNA-mediated repression of mRNA by splicing-regulated and highly repressive microRNA recognition element. BMC Genomics. 2013;14:438.
Elton TS, Martin MM. Alternative splicing: a novel mechanism to fine-tune the expression and function of the human AT1 receptor. Trends Endocrinol Metab. 2003;14(2):66–71.
Sanborn BM, Millan JL, Meistrich ML, Moore LC. Alternative splicing of CREB and CREM mRNAs in an immortalized germ cell line. J Androl. 1997;18(1):62–70.
Stoss O, Stoilov P, Hartmann AM, Nayler O, Stamm S. The in vivo minigene approach to analyze tissue-specific splicing. Brain Res Brain Res Protoc. 1999;4(3):383–94.
The authors are very grateful to Dr John Milton (School of Animal Biology, University of Western Australia) for designing the nutrition treatments, and to Mr Matt Wilmot (CSIRO, Perth, Western Australia), Dr César Rosales Nieto and Mr Gary Cass (School of Animal Biology, University of Western Australia) for their assistance in handling the sheep. Mr Sheng Zhang (Zhejiang University, China), Mrs Margaret Blackberry, Dr Trina Jorre de St Jorre and Mr Fahad Almohsen (School of Animal Biology, University of Western Australia) provided valuable assistance with the collection and preservation of tissue.
This project was funded by Canadian government and UWA School of Animal Biology, and Yongjuan Guan was financially supported by a Scholarship for International Research Fees from the University of Western Australia.
YG and GL performed the molecular analysis, data analysis and draft writing; YG and GBM contributed to animal study; LLG contributed to interpretation of the molecular data. All authors contributed to manuscript writing. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
All the sequencing data were deposited in the publicly available NCBI’s Gene Expression Omnibus Database. The data are accessible through GEO Series accession number GSE68274 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE68274).
Phenotypic traits of sheep testes used for WGCNA analysis. “H” indicates sheep fed with high diet, “L” indicates sheep fed with low diet. P1 indicates testis weight (g), P2 indicates sperm number per testis, P3 indicates diameter of seminiferous tubule, P4 indicates volume of seminiferous epithelium (x 1012 μm3), P5 indicates change of scrotal circumference (cm) and P6 indicates apoptotic germ cells/tubule. (XLSX 10 kb)
Transcription profiles plotted across the sheep genome, showing the distribution of the RNA-seq read density along the length of each chromosome. Each vertical blue line represents log2 of the frequency of reads plotted against the chromosome coordinates. (PDF 880 kb)
Differentially expressed mRNAs in testis from sheep fed a low or high diet (N = 8 for each treatment). Note: Fold change (FC) = CPM of low diet group/CPM of high diet group. CPM (Counts per million) = (mRNAs reads number/total reads number per library) × 1,000,000. The significant DE mRNAs were determined by false discovery rate (FDR) < 0.05. (XLSX 168 kb)
The 10 most related function clusters of differentially expressed mRNAs by DAVID analysis. The lower p value indicates greater relevance. (XLSX 12 kb)
qRT-PCR validation of differentially expressed genes. mRNA expressions from qRT-PCR are shown by line graphs on the top and values are shown on the right Y-axis as relative expression (2-ΔΔCt). mRNA expressions from RNA-Seq are shown by bar graphs on the bottom and values are shown on the left Y-axis as log2 (normalised reads number). a, b - indicate the significant difference in the relative expression of mRNAs detected via qRT-PCR at P<0.05; A, B - indicate significant difference in the expression of mRNAs detected from RNA-seq at FDR <0.05. Data are presented as Mean±Standard deviation. (PDF 610 kb)
Identification of putative miRNA-mRNA pairs on the basis of target prediction and the negative regulatory effect of miRNAs on mRNA expression levels. (XLSX 26 kb)
Regulatory relationships for two pairs of miRNAs and mRNAs that were differentially expressed in sheep testis following nutritional treatment: oar-novel-miR-33 with 68 mRNAs, and oar-novel-miR-31 with 52 mRNAs. (PDF 937 kb)
Identification of alternative splicing events with Tophat2 software. (XLSX 1164 kb)
About this article
Cite this article
Guan, Y., Liang, G., Martin, G.B. et al. Functional changes in mRNA expression and alternative pre-mRNA splicing associated with the effects of nutrition on apoptosis and spermatogenesis in the adult testis. BMC Genomics 18, 64 (2017). https://doi.org/10.1186/s12864-016-3385-8
- Pre-mRNA alternative splicing