This study was conducted at the Teagasc Animal and Grassland Research and Innovation Centre. All procedures involving animals were approved by the Teagasc Animal Ethics Committee and licensed by the Irish Health Products Regulatory Authority in accordance with the European Union Directive 2010/63/EU. All experiments were performed in accordance with relevant regulations and the ARRIVE (Animal Research: Reporting on In Vivo Experiments) guideline.
Experimental design and animal management
The animal model and general management regimen pertaining to the animals used in this study have previously been extensively described [10]. Briefly, thirty Holstein-Friesian bull calves with a mean (±SD) age and bodyweight of 17.5 (± 2.85) days and 48.8 (± 5.3) kg, respectively, were sourced from four commercial dairy farms and blocked on age, sire, liveweight and farm of origin. After seven days of acclimatisation, calves were assigned to either a high (HI) or moderate (MOD) plane of nutrition. Calves were individually fed milk replacer and pelleted concentrate using an electronic calf feeding system (Forster-Tecknik Vario; Engen, Germany). Calves on the HI received 1500 g of milk replacer in 10 L of water daily, together with concentrate ad-libitum. Calves on MOD were allocated 500 g of milk replacer in 4 L of water plus a maximum of 0.5 kg of concentrates daily. All calves had daily access to fresh water and approximately 0.5 kg of hay. At a mean age (±SD) of 87 (±2.14) days, all calves were euthanized following intravenous administration of an overdose of sodium pentobarbitone (300 mg/ml: 0.25 mL/kg bodyweight). The testes were excised from each calf and the tunica albuginea, epididymides and any excess connective tissue removed. The testes parenchyma tissue was harvested avoiding the rete testis and snap-frozen in liquid nitrogen and subsequently stored at −80°C pending further processing.
RNA isolation
Total RNA was isolated from all testes samples (n=30) using the Qiagen RNeasy Plus Universal Mini kit (Qiagen Ltd, Manchester, UK) according to the manufacturers instructions including steps for the purification of Total RNA containing miRNA . Yield of the resultant RNA was determined by measuring the absorbance at 260 nm with a Nanodrop spectrophotometer (NanoDrop Technologies; Wilmington, DE, USA). The quality of the isolated RNA was evaluated on an Agilent 2100 Bioanalyser (Agilent Technologies Ireland Ltd., Dublin, Ireland). All samples were deemed to be of excellent quality for subsequent RNA sequencing, with an overall average RNA Integrity Number of 9.5.
RNAseq library preparation and sequencing
Both mRNA and miRNAseq library preparation and sequencing were undertaken by a commercial sequencing facility (Macrogen Europe Inc., Amsterdam, The Netherlands). Individual cDNA libraries were prepared from all 30 total RNA samples using the Illumina Truseq stranded mRNA kit for mRNAseq and the Illumina small RNA kit for miRNA, according to the kit instructions. High-throughput sequencing was undertaken on an Illumina NovaSeq, incorporating 150 bp paired end sequencing for mRNA and an Illumina HiSeq 2500 sequencer, incorporating 50 bp single end sequencing for small RNA sequencing.
miRNAseq and mRNAseq bioinformatic analyses
Both mRNA and miRNA raw sequence reads were firstly assessed for sequencing quality using FASTQC (version 0.11.8, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The Illumina sequencing adaptor was clipped off all the raw read sequences using Cutadapt (version 1.18). For the miRNA data, reads of lengths shorter than 15 bp, and longer than 28 bp were subsequently removed as short and long reads, respectively, using Cutadapt. The retained reads were then additionally filtered for other bovine short RNA species including ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs) and small nucleolar RNAs (snoRNAs) downloaded from https://rnacentral.org/. To profile miRNA expression in each sample, the miRDeep2 package (version 2.00.8) modules were used, together with the bovine reference genome (ARS-UCD1.2) and the known bovine mature miRNA sequences and their precursor sequences from the miRBase database (release 22.1). The miRDeep2 mapper module (mapper.pl) was then used with default parameters to collapse reads of the sequences into clusters. Bowtie (version 1.1.1) was then employed to align the collapsed reads to the indexed reference genome. Using default parameters and input files including the reference genome, collapsed reads versus reference genome alignment, known bovine mature miRNAs and their precursors sequences (including the hairpin structures) and Bos taurus (bta) as the species of interest, the miRDeep2 module (miRDeep2.pl) was used to quantify bovine miRNAs. Through this, the miRDeep2 quantifier module was used to quantify all known expressed miRNAs in the sequence data, producing read counts for each individual sample. For the mRNA sequencing data, the Spliced Transcripts Alignment to a Reference (STAR) aligner [11] was employed to align sequencing reads to the ARS-1.2 version of the bovine reference genome. Additionally, within the STAR alignment script, the quantmode function was included in order to quantify the number of reads aligned to each gene. Resultant read counts for each sample across each mRNA and miRNA datasets were then merged into two separate files (miRNA reads and mRNA reads) and subsequently assessed separately for differentially expressed miRNA and mRNA using the R (v2.14.1) Bioconductor package, edgeR (version 3.26.7), as previously described [7]. Target genes for differentially expressed miRNAs were predicted using TargetScan (release 7.2; http://www.targetscan.org/vert_72/).
Global proteomics analysis
Sample preparation
The methodology has been previously described in [12] . Global proteomics analysis was undertaken on the same tissue samples as used for RNAseq. For each sample, proteins were extracted using a tissue homogenizer (TissueLyser II, Qiagen) and digested using a commercial iST Kit (PreOmics, Germany). Briefly, 50 µl of ‘Lyse’ buffer and around 50 ug of glass beads (425-500 µm, Sigma) were added to the thawed tissue. After 2 cycles of protein extraction (2 min each, 30 Hz) and 10 min at 95°C, the tubes were centrifuged at 16,100 × g for 15 min, and the supernatant was pipetted into a fresh Eppendorf tube for proteomics analysis. The solubilization of the extracted proteins was enhanced by subjecting the samples with High Intensity Focused Ultrasound (HIFU) for 1 min setting the ultrasonic amplitude to 85%. The protein concentration was estimated using the Qubit® Protein Assay Kit (Life Technologies, Zurich, Switzerland). For each sample, 50 µg of protein were transferred to the cartridge and digested by adding 50 µl of the ‘Digest’ solution. After 60 min of incubation at 37 °C the digestion was stopped with 100 µl of ‘Stop’ solution. The solutions in the cartridge were removed by centrifugation at 3,800 g, while the peptides were retained by the iST-filter. Finally, the peptides were washed, eluted, dried and re-solubilized in 20 µL of injection buffer (3% acetonitrile, 0.1% formic acid).
Liquid chromatography-mass spectrometry analysis
Mass spectrometry (MS) analysis was performed on a Q Exactive HF-X mass spectrometer (Thermo Scientific) equipped with a Digital PicoView source (New Objective) and coupled to a M-Class UPLC (Waters) [12] . Solvent composition at the two channels was 0.1% formic acid for channel A and 0.1% formic acid, 99.9% acetonitrile for channel B. For each sample, 1 μL of peptides was loaded on a commercial MZ Symmetry C18 Trap Column (100Å, 5 µm, 180 µm x 20 mm, Waters) followed by nanoEase MZ C18 HSS T3 Column (100Å, 1.8 µm, 75 µm x 250 mm, Waters). The peptides were eluted at a flow rate of 300 nL/min by a gradient from 8 to 27% B in 85 min, 35% B in 5 min and 80% B in 1 min. Samples were acquired in a randomized order. The mass spectrometer was operated in data-dependent mode (DDA), acquiring a full-scan MS spectra (350−1’400 m/z) at a resolution of 120’000 at 200 m/z after accumulation to a target value of 3,000,000, followed by HCD (higher-energy collision dissociation) fragmentation on the twenty most intense signals per cycle. HCD spectra were acquired at a resolution of 15’000 using a normalized collision energy of 25 and a maximum injection time of 22 ms. The automatic gain control (AGC) was set to 100’000 ions. Charge state screening was enabled. Singly, unassigned, and charge states higher than seven were rejected. Only precursors with intensity above 250,000 were selected for MS/MS. Precursor masses previously selected for MS/MS measurement were excluded from further selection for 30 s, and the exclusion window was set at 10 ppm. The samples were acquired using internal lock mass calibration on m/z 371.1012 and 445.1200.
Protein identification and label free protein quantification
The acquired raw MS data were processed by MaxQuant (version 1.6.2.3), followed by protein identification using the integrated Andromeda search engine [13]. Spectra were searched against a Uniprot Bos taurus reference proteome (taxonomy 9913, version from 2017-08-17), concatenated to its reversed decoyed fasta database and common protein contaminants. Carbamidomethylation of cysteine was set as fixed modification, while methionine oxidation and N-terminal protein acetylation were set as variable. Enzyme specificity was set to trypsin/P allowing a minimal peptide length of 7 amino acids and a maximum of two missed cleavages. MaxQuant Orbitrap default search settings were used. The maximum false discovery rate (FDR) was set to 0.01 for peptides and 0.05 for proteins. Label free quantification was enabled and a 2-minute window for match between runs was applied. In the MaxQuant experimental design template, each file is kept separate in the experimental design to obtain individual quantitative values. Protein fold changes were computed based on Intensity values reported in the proteinGroups.txt file. A set of functions implemented in the R package SRMService [14] was used to filter for proteins with 2 or more peptides allowing for a maximum of 10 missing values, and to normalize the data with a modified robust z-score transformation and to compute p-values using the t-test with pooled variance. If all measurements of a protein were missing in one of the conditions, a pseudo fold change was computed replacing the missing group average by the mean of the 10% smallest protein intensities in that condition.
Integration of miRNA, mRNA and proteomic results
In order to evaluate the relationship between miRNA, mRNA and proteins, the various ‘omics’ datasets pertaining to the current study were integrated. This was achieved by firstly determining a relationship between the miRNA and mRNA results; specifically, the mRNA targets of the DE miRNA were assessed for the presence of DE mRNA genes. A weighted co-regulatory network analysis was then undertaken on the proteomics data and resultant networks mined for DE genes (which were affected/targeted by DE miRNA) using the WGCNA R software program [15]. Label-free quantitation intensity values were Log2 transformed in R. The WGCNA automatic network construction and module detection method was utilised to generate unsigned co-regulated networks. For each separate analysis, pair-wise weighted Pearson correlations were calculated between all pairs of proteins in the testes dataset. Adjacency matrices were calculated to reach scale-free topology of the network (R2>0.09) by raising the co-regulation matrix to a soft-threshold power of 14. Following this, the topology overlap matrix was calculated, providing information on the similarity of the co-regulation between two proteins with all other proteins in the network. Average linkage hierarchical clustering was then applied to the topology overlap matrix resulting in the grouping of modules of co-regulated proteins. Resulting modules or networks of co-expressed proteins were mined for DE mRNA genes affected by DE miRNAs to determine the interactions of these transcripts with the testicular proteome as affected by prevailing early life nutrition.
Gene and protein pathways analysis
In order to determine biological pathways and functions enriched by the dietary treatments imposed, miRNA target genes, DE mRNA genes and protein network analysis results were individually subjected to pathway analysis using Ingenuity Pathway Analysis (IPA) according to the manufacturer’s instructions (Qiagen Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis [15].