Enhanced pre-pubertal nutrition upregulates mitochondrial function in testes and sperm of post-pubertal Holstein bulls

Supplemental energy and protein during calf-hood (2–30 wk) in dairy bulls hastened puberty (~1 mo), upregulated steroid biosynthesis, concentrations of reproductive hormones and Sertoli cell maturation, with larger testes and greater sperm production (~25%) in mature bulls. The objective was to evaluate effects of feeding high (20.0% crude protein [CP], 67.9% total digestible nutrients [TDN]), control/medium (17.0% CP, 66.0% TDN) and low (12.2% CP, 62.9% TDN) diets from 2 to 30 wk on post-pubertal testes of Holstein bulls. Based on RNA sequencing, 497 and 2961 genes were differentially expressed (P < 0.1) in high- vs low- and high- vs medium-diet groups, respectively. According to KEGG analysis, oxidative phosphorylation and ribosome pathways were upregulated in high- vs medium- and low-diet groups, with majority of upregulated genes encoding for essential subunits of complex I, III, IV and V of OXYPHOS pathway. In addition, mitochondrial translation, mitotic nuclear division and cell division were enriched in high- vs medium-diet groups. Consistent with these results, a greater percentage of sperm from high-diet bulls were progressively motile and had normal mitochondrial function compared to medium-diet sperm (P < 0.1). Thus, enhanced early life nutrition upregulated mitochondrial function in testes and sperm of post-pubertal Holstein bulls.

Differentially expressed genes. In general, of 18288 genes expressed, 2960 were differentially expressed (P < 0.1) between HD and MD groups compared to 496 in HD versus LD groups (Supplementary Datasets 2 and 3). A scatterplot denoting transcriptomic differences across diet groups is provided (Supplementary File 1). To visualise extent of differential expression and relationship between mean gene expression levels and log2fold changes, MA plots were generated (Fig. 1). Of the 496 DEGs in the HD versus LD analysis, 446 genes were upregulated in the HD group (genes with higher log2fold change were: COX5B, COX6B1, ATP5J, ATP5G1, MRPL11 and HSD17B1). When comparing HD and MD groups, of the 2960 DEGs, 1828 genes were upregulated in the HD group (genes with higher log2fold change: COX2, COX3, COX5B, HSD17 and CCNA1). The log2fold changes varied from ± 0.21 to ± 1.17 among all DEGs. No genes were differentially expressed between MD and LD groups. Therefore, all further analyses were done using two datasets, HD vs LD and HD vs MD groups.
Gene Ontology term enrichment. To systematically examine enriched genes associated with differential pre-pubertal diets on post-pubertal testes, gene ontology (GO) analysis was done using the functional annotation and clustering option in DAVID. "Mitochondrial translation initiation" and "elongation" were the major biological functions enriched in HD vs LD and HD vs MD groups. In the HD vs MD group comparison, "mitotic nuclear division", "cytoplasmic translation" and "ATP synthesis coupled protein transport" were also enriched. Complete lists of GO terms enriched in HD vs LD and HD vs MD datasets are provided in Tables 1 and 2, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/ Over-represented KEGG pathways. Eight pathways were associated with HD vs MD groups and six pathways in HD vs LD groups when considering genes (P < 0.1, Tables 3 and 4). "Oxidative phosphorylation" and "ribosome pathways" were major pathways enriched in both datasets.
In addition to above, oxidative phosphorylation and ribosome pathways were enriched in the upregulated genes from both datasets (HD vs MD and HD vs LD, Supplementary Files 2 and 3). However, nine downregulated pathways were identified in the HD vs MD dataset (Supplementary File 4). Major pathways identified were focal adhesion, regulation of actin cytoskeleton and axon guidance.      www.nature.com/scientificreports www.nature.com/scientificreports/ PI-(P < 0.1) and Rhl23-/PI + sperm (P < 0.05; Table 6). Percentage of Rhl23−/PI− sperm (live sperm loosing mitochondrial function) was higher in MD vs HD groups (15.1 ± 2.6 vs 8.6 ± 0.9). However, there was no significant difference among groups for percentage of Rhl23 + /PI-sperm (live sperm with normal mitochondrial function). The percentage of Rhl23-/PI + sperm (dead or necrotic sperm) in the MD was 66.2 ± 3.5 and was lower (P < 0.05) than both LD and HD groups (78.6 ± 2.2 and 78.7 ± 2.2, respectively). A representative image of flow cytometry explaining sperm populations is provided (Fig. 4).

Discussion
To provide a more comprehensive assessment of post-pubertal effects of pre-pubertal feeding, RNA-sequencing was done on testicular tissues collected at 72 wk of age. To our understanding, this is the first study to evaluate effects of differential feeding during early pre-pubertal period on the post-pubertal bovine transcriptome using a next-generation sequencing technique. Bulls on the HD (fed supplemental protein and energy) had greater oxidative phosphorylation and mitochondrial protein synthesis within their testes compared to LD and MD bulls. Furthermore, sperm from HD bulls had greater progressive motility and a greater percentage had normal mitochondrial function compared to sperm from MD bulls.
Oxidative phosphorylation (OXYPHOS) occurs in mitochondria, generating nearly 80% of cellular energy. OXYPHOS, the functional unit of mitochondria, connects the electron transport chain (ETC) to cell respiration and ATP synthesis. Enzyme complexes constituting OXYPHOS are comprised of mitochondrial and nuclear encoded polypeptides 9 . Elevated mitochondrial and cytoplasmic translation in HD vs MD were indicative of increased protein availability for OXYPHOS formation. A similar study in Holstein bulls also identified increased OXYPHOS in adipose tissue along with enhanced testes development, cholesterol and androgen biosynthesis in the testes of pre-pubertal bulls fed high vs low planes of nutrition [10][11][12] .
The integral part of OXYPHOS, the ETC, is composed of five multi-protein complexes and two electron carriers (ubiquinone and cytochrome) transporting electrons to oxygen, creating a proton motive force. The resulting membrane potential enables energy production when these protons flow back via ATP synthase 13 . Proton transfer at sites other than ATP synthase reduces mitochondrial efficiency, but is critical for ROS production. Testicular mitochondria consume less oxygen than other tissues to maintain the same electric potential and have higher phosphorylation efficiency with age 14 . In mature testes, spermatocytes and spermatids use OXYPHOS, due to limited glucose availability in seminiferous tubule fluid, whereas Sertoli cells, spermatocytes and mature sperm have high glycolytic activity 15,16 . The majority of upregulated genes in HD (vs MD and LD) were significantly enriched for OXYPHOS.
In sperm, mitochondrial activity has been associated with sperm motility 17 , viability 18 , capacitation 19 and fertilization potential 20 . Within germ cells, mitochondrial ATP synthesis was reported to regulate the apoptotic pathway 21 . Higher MMP has been associated with increased sperm motility and production, as sperm require ATP 17 . In addition, MMP is correlated with activity of the ETC and ATP synthase (components of the OXYPHOS). Increased mitochondrial activity beyond its threshold could induce oxidative stress, which could alter DNA. The majority of genes involved in the ETC had high expression in HD bulls compared to MD and LD bulls. In addition, major mitochondrial functions, including ATP synthesis coupled protein transport, hydrogen ion transmembrane transport, mitochondrial electron transport and oxidation reduction processes were altered in HD vs MD bulls.
A positive association between mitochondrial function and motility has been reported in bovine sperm 22 . In our study, sperm from HD bulls had greater progressive motility (PM) compared to sperm from MD and LD bulls; the association between motility and mitochondrial function, reflected our gene expression study. Furthermore, positive associations between OXYPHOS, ETC and sperm motility have been reported 18,23,24 . We also evaluated MMP, for which Rhodamin123, a mitochondria-specific fluorescent dye in combination with PI, was used 25 . When comparing live cells, a higher percentage of sperm from MD bulls had lowered MMP than sperm from HD bulls. Surprisingly, percentages of dead/necrotic sperm were significantly elevated in both HD and LD groups compared to the MD group (Table 6). An increased percentage of dead sperm in our HD group was reported in our phenotypic study 8 . A possible explanation could be the increased ROS production favoured by the elevated mitochondrial function in high diet bulls promoting cell death 26 . Both progressive motility and MMP results were consistent with our gene expression study, where the majority of genes involved in mitochondrial function were upregulated in HD versus MD groups. Positive associations between MMP and motility have also been reported in sperm 22,27 .
In addition to altered ETC function, majority of genes associated with ribosome pathway were upregulated in HD vs MD groups. We reported that supplemental energy and protein in diets of pre-pubertal bull calves hastened puberty and increased both testes size and sperm production potential, with no significant differences in classical assessments of sperm quality or function 4,8 . The role of mitochondrial function in regulating various aspects of sperm function, including meiosis, spermatogenesis, energy for survival, motility, capacitation, maturation and quality control has been reviewed 28 . Increased mitochondrial activity is essential to support elevated www.nature.com/scientificreports www.nature.com/scientificreports/ sperm production potential previously reported in HD bulls. High ribosomal biogenesis is also intimately associated with cellular activities, including cell growth and division 29 . All genes involved in the ribosome pathway were upregulated in the HD group, implying high ribosomal biogenesis.
In addition to its specific role in sperm function, mitochondria also have an endocrine role in testes 30 . The rate limiting and first step in steroid biosynthesis, converting cholesterol to pregnanolone, takes place within mitochondria and is catalysed by the enzyme cytochrome P450 side cleavage 31 . StAR promotes cholesterol movement in mitochondria (from outer to inner membranes) and after pregnanolone is formed, it is transported to endoplasmic reticulum to complete steroid biosynthesis 32 . Increased steroid production (mainly testosterone) is critical for male reproductive development and sperm production.
In addition to being critical for cell growth and metabolism, mitochondrial function is closely associated with epigenetics 33,34 . Epigenetic changes can be induced by various factors in the environment, e.g. nutrition, that can alter gene expression by incorporating changes into the epigenome. The epigenome is a record of chemical changes on the DNA and histone proteins that may be transferred to the next generation 35,36 . Substantial epigenetic changes, including methylation of DNA methylation and modification of histones, use various metabolic intermediates, including s-adenosyl methionine (SAM) and acetyl coA 34 . Based on the energy status of a cell, these metabolites can vary and are to an extent regulated by mitochondrial function. For instance, energy production within mitochondria converts dietary calories into ATP, acetyl coA, SAM and reduced NAD. With a high-energy diet, there is phosphorylation and acetylation of chromatin by ATP and acetyl-CoA, thereby making nDNA more available for transcription and replication. However, if dietary energy is limited, phosphorylation and acetylation of chromatin do not proceed, thereby supressing gene expression. Methylation of DNA facilitated by SAM can also be regulated by mitochondrial function 37,38 . Thus, alterations in mitochondrial function could cause epigenetic changes with long-term consequences. In addition, recent studies in our lab are indicative of DNA methylation changes in sperm of bulls fed differential pre-pubertal diets (unpublished data). It is possible that the sperm epigenome keeps memory of diet during pre-pubertal period in genes critical for sperm function. However, more and more studies in humans have demonstrated the ability of paternal environment to influence the epigenetic/non-genetic information transmission carried by sperm cells and modify the trajectory of offspring [39][40][41] .
To summarise, enhanced pre-pubertal feeding in Holstein bulls upregulated genes associated with mitochondrial function, including oxidative phosphorylation, ETC and mitochondrial protein synthesis. In addition to enhanced mitochondrial function, elevated mitotic activity and ribosome synthesis could have supported greater sperm production potential of HD bulls, as evidenced in our phenotypic study. In addition to greater progressive motility, percentage of live sperm losing mitochondrial function was much lower in HD vs MD groups. We concluded that enhanced pre-pubertal nutrition followed by a control diet enhanced mitochondrial function and sperm quality in Holstein bulls. Bull phenotype. Body weights, scrotal circumference and paired testes volume from 11-71 wk of age have been reported 4 . Briefly, high-diet bulls had greater body and paired testes weights at slaughter (71 wk). Sperm production potential was significantly greater in HD versus LD bulls, but there were no significant differences among groups regarding age when post-thaw semen first achieved minimum standards (based on progressive motility and morphology) or on sperm viability, protein content or in vitro fertilization potential 8 . total RNA extraction. Testicular tissue (75 mg) was homogenized and total RNA extracted using Trizol reagent, as per manufacturer's instructions (Invitrogen, Canada, Burlington, ON, Canada). Total RNA was suspended in water, quality and quantity measured using a Nanodrop via spectrophotometry (2000/2000c, Thermo Fisher Scientific, Waltham, MA, USA). The RNA integrity number (RIN) was assessed on an Agilent 2100 bioanalyzer using an Agilent RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA). mRNA sequencing of testicular tissue. RNA-sequencing library preparation and sequencing. For each sample (n = 24), total RNA (1 µg) was used to create mRNA libraries using the TruSeq Stranded mRNA Library Prep Kit (Illumina Inc., San Diego, CA, USA) and unique adapter indices (Illumina) at the University of Lethbridge sequencing facility (Lethbridge, AB, Canada). Libraries were sequenced as 75 base pair single end reads using the NextSeq. 500 system (Illumina Inc.). For this purpose, individual libraries were pooled to obtain a single library by diluting individual libraries to a final concentration of 4 nM and adding equal volumes of all samples (n = 24). Base calling and demultiplexing was performed with CASAVA 1.9 (Illumina Inc.), using default settings. Quality of sequencing reads was assessed using FastQC software 42 . Quality of sequencing libraries was considered adequate to proceed with downstream analysis without adapter and quality trimming. Very limited adapter contamination was identified.
Mapping and annotation of reads. Sequencing reads were mapped to the bovine genome (Bos taurus, UMD3.1.87, downloaded from Illumina iGENOME website) using TopHat v2.0.10 43 and with Bowtie v.1.1.2 44 used as an internal aligner. Mapping rate was uniformly good, with >93% of reads successfully mapped to the genome.
Detection of differentially expressed genes. Counts of reads mapping to genes were obtained using featureCounts software from Subread package 45 . Count data were loaded into R, with normalization and variance stabilizing transformation 46 applied using DESeq2 package. Normalized and variance stabilized data were used for quality control and sample clustering analyses. Euclidean sample distances were calculated in R and sample clustering was done using hclust 47 function in R. Sample to sample distances were visualized as heatmaps. Principal components plots were made with PCAplot function in DESeq2 package 48 .
Genes differentially expressed between diets were detected using Wald tests on normalized count data. Experimental conditions were compared with one-way ANOVA; thereafter, results for specific contrasts (HD vs LD, HD vs MD and MD vs LD) were extracted. Correction procedures for multiple comparisons were done with Benjamini-Hochberg correction implemented in R. Genes with adjusted P-values < 0.1 were considered differentially expressed. Functional analysis. The Ensembl identities of DEGs from HD vs LD and HD vs MD analyses were uploaded to Database for Annotation, Visualization and Integrated Discovery (DAVID) to enrich for Gene Ontology (GO) Biological process and cellular components 49 . For most analyses, we used the functional annotation and clustering option and reported all GO terms that were P < 0.05 and molecule number ≥2.
Quantitative real-time polymerase chase reaction (RT-PCR). RT-PCR was performed as described previously 7 . In brief, reverse transcription and subsequently PCR with SYBR Green (Fast SYBR ® Green Master Mix; Applied Biosystems, Foster City, CA, USA) was done to validate mRNA expression of five differentially expressed genes: TJP2, COX5B, ATP5J, ATP5G1 and CRHR1. The National Center for Biotechnology and

Data availability
All data generated or analysed during this study were deposited in the publicly available NCBIs Gene Expression Omnibus Database GSE137690 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE137690.