Non-invasive metabolomics biomarkers of production efficiency and beef carcass quality traits

Study design and samples details

Research protocols were approved and monitored by the USDA, ARS, U.S. Meat Animal Research Center Institutional and Animal Care Committee in accordance with the Guide for the Care and Use of Agricultural Animals in Agricultural Research and Teaching (2010). All authors complied with the ARRIVE guidelines. Calves (n = 75) used in this experiment were the steer progeny of Red Angus, Charolais, and Simmental sires bred by artificial insemination to four dam breeds: MARC II Composite (¼ Hereford, ¼ Angus, ¼ Gelbvieh, ¼ Simmental), MARC III Composite (¼ Hereford, ¼ Angus, ¼ Pinzgauer, ¼ Red Poll), ½ MARC II X ½ MARC III, and ½ MARC II X (¼Red Angus X ¼ Angus).

Stage-1 (forage based-diet)

After weaning, steers received an implant containing 80 mg of trenbolone acetate and 16 mg estradiol (Revalor-IS; Merck), and were housed in a facility with Calan Broadbent electronic headgates (American Calan, Inc.,Northwood, NH) to measure individual feed intake. Steers were trained to use Calan Headgates during an adaptation period (approximately 21 days).

At the beginning of the study (0 day), study population mean age and weight of cattle (± standard deviation) was 271.4 ± 5.0 days and 301.7 ± 31.9 kg, respectively. The diet consisted of 30.0% chopped alfalfa, 69.8% corn silage, and 0.2% salt on a DM basis.

Feed intake was measured for 84 days with cattle weighed on day 0, 1, 21, 42, 63, 83 and 84 of the experiment. Daily feed offered was recorded and orts determined weekly. A urine sample was collected at day 0, 42, and 83 during the scheduled weigh dates.

Stage-2 (concentrate based-diet)

Following the completion of stage 1, steers were transitioned to a concentrate-based ration over a 21-day period and remained on that diet for the remainder of stage 2. On a DM basis, the diet consisted of 8.00% chopped alfalfa, 67.75% rolled corn, 20.00% wet distillers grain with solubles, and 4.25% vitamin/mineral supplement with Rumensin (Elanco, Greenfild, IN). Forty-two days after stage 1 steers received an implant containing 120 mg of trenbolone acetate and 24 mg of estradiol (Revalor-S, Merck Animal Health, Madison, NJ) and feed intake and body weight gain data were collected for 84 days. Body weight measurements were determined at 0, 1, 21, 42, 63, 83 and 84 days on study. Daily feed offered were recorded and orts determined weekly. At the beginning of stage-2 (0 day), mean age and weight of cattle (± standard deviation) was 397.3 ± 5.0 days and 448.6 ± 41.8 kg, respectively. A urine sample was collected at day 0, 42, and 83 during the scheduled weigh dates.

Urine collection

On days 0, 42, and 83 for each stage-1/2 urine samples (~ 50 mL) were collected by preputial stimulation and immediately centrifuge 10,000g for 10 min to remove particulates (4 °C) before storage at − 80 °C. No urine sample could be collected by stimulation 23% (n = 19) of the selected steers in stage-1 and 11% (n = 9) of the selected steers in stage-2.

Steer selection by divergent growth for each stage

Steers were selected according to differences in ADG at each stage with those with the greatest and least ADG, whose dry matter intake was within 0.32 SD of the mean intake (Fig. 7). After the end of intake study (stage-2), selected steers in both stages received the same ration ad libitum and remained in the same pen until slaughter (5–8 days) at an abattoir. Carcass data was collected using the VBG2000 beef carcass grading camera (Vision For You, L.L.C., Dakota Dunes, SD).

Sample preparation for untargeted metabolomics analysis

All urine sample extractions and LC–MS analysis were randomized across time sampling to avoid batch effects. Duplicate samples of each urine were extracted for metabolomics analysis within batches. Both extraction and chromatographic solvents were UPLC-MS Optima Grade (Fisher Chemical Ltd., Waltham, MA). Sample preparation was adapted from previous reports⁴⁶. Briefly, 50 µL urine aliquots were diluted with 150 µL of water, vortexed and centrifuged at 16,000×g for 10 min at 4 °C. The supernatant was filtered at 0.22 µm by centrifugation in Spin-X LC (Fisher Chemical Ltd., Waltham, MA) membrane-tube filters (5000×g for 5 min at 4 °C) before transfer to LC-vials.

UPLC-qTOF spectral acquisition of urine samples

The UPLC/MS analysis was performed using a Waters ACQUITY ultra-performance liquid-chromatography (UPLC) system (Waters Corporation, Milford, MA) equipped with an auto sampler and coupled with a hybrid triple quadrupole-time of-flight mass spectrometry (XEVO-G2-S-qTOF; Waters). Instrument calibration was performed before running the samples using 0.5 nM sodium formate. To obtain information regarding system suitability and stability, quality control (QC) samples were injected at regular intervals (every ten samples) throughout the analytical run. Quality control samples were prepared by mixing ~ 100 urine aliquots (10 µL), producing separate QC samples.

For analyses, samples were held at 4 °C, and analytes in 10 µL of urine extracts were separated on a 2.1 × 50 mm × 1.7 µm Acquity UPLC BEH C18 column (Waters Corp.) held at 40 ºC using a 0.4 mL/min flow of 0.1% formic acid in water (mobile phase A) and 0.1% formic acid in acetonitrile (mobile phase B) with the following gradient: 99% A from 0 to 2 min; 85% A at 4 min; 50% A at 8 min; 5% A at 12 min; 99% A at for 16 min followed by 2 min of re-equilibration. All gradient ramps were linear.

Mass spectrometry was performed in both positive and negative modes. The capillary voltage was 3.2 kV and 2.4 kV for positive and negative mode, respectively. The system parameters were set as follows: source of temperature 120 °C, desolvation temperature 350 °C, and cone gas flow (nitrogen) 25 L/h and desolvation gas flow (nitrogen) 900 L/h. Data were collected in a centroid mode using the lockspray to ensure accuracy and reproducibility. Leucine enkephalin was used as lock-mass in a 2 ng/mL concentration solution. The lockspray frequency was set at 15 s, and the lock mass data were average over 15 scans for correction. The scan mass range was from 50 to 1200 m/z using an extended dynamic range. The MS/MS analysis was carried out by ramping the collision energy from 10 to 50 V using argon as a collision gas.

Untargeted-metabolomics data processing

Raw data obtained from UPLC-qTOF analysis was analyzed using Progenesis QI v1.0 software (Waters Corp.). The data was aligned, deconvoluted, and normalized using total ion intensity. Solvent blanks were run between samples and each mass was checked against the blank run to exclude possible sources of contaminations. The coefficient of variation (CV) of the mass abundance for all biological replicates, was calculated across pooled repeated samples for each feature and those with a CV above a 20% cutoff were removed. Metabolites were identified by comparison with online Bovine Metabolome Database (http://www.cowmetdb.ca/) using exact m/z values and retention times. The identities of selected metabolites were confirmed by MS/MS fragment ion analysis using Mass-Fragment application manager software (Water MassLynk v4.1, Waters Corp.). The MS/MS fragmentation of the molecules was compared, with ChemSpider database (www.chemspider. com), by way of chemically intelligent peak-matching algorithms.

Urine bile acid and steroid quantification

Endogenous urinary steroids and bile acids were quantified using isotopic enrichment and authentic calibration standards detected by LC–MS/MS with electrospray ionization and multiple reaction monitoring (MRM) on an API 6500QTRAP (Sciex, Framingham, MA, USA) as previously reported for plasma⁴⁷. Urine extraction was performed by solid-phase extraction (SPE) using Oasis-HLB SPE 96-wellplates 30 mg (Milford, MA, USA). Briefly, 100 µL urine aliquots were spiked with 10 µL of methanol containing a suite of deuterated bile acids and steroids at 12.5 nM, vortexed, and loaded onto SPE cartridges pre-conditioned with 1 mL MeOH, followed by 1 mL H₂O. Loaded cartridges were then washed with 1 mL H₂O and eluted with 2 mL MeOH. The eluate was evaporated under vacuum at room temperature and reconstituted in 100 µL of 50:50 MeOH:H₂O (v/v). Extracts were separated on a 2.1 × 100 mm, 1.7 μm Acquity C18 BEH column using a Shimadzu Nexera X2 UPLC (Shimadzu, Kyoto, Japan). Analytes were quantified using internal standard methods with surrogate/analyte associations for response ratio calculations. Data was processed with AB Sciex MultiQuant v 3.0.1.

Statistical and machine learning data analysis

Animal traits

The feed conversion ratio (FCR) was calculated as the quotient of DMI divided by ADG. The residual feed intake (RFI) was calculated as the residual obtained from the linear regression of DMI on ADG and calculated metabolic body weight at the midpoint of each stage. Feed intake, body weight traits, FCR, RFI, age, and carcass traits were analyzed within stage using the mixed procedure of SAS (Cary, NC). The model included ADG classification as a fixed effect and animal nested within classification as a random effect with the denominator degrees of freedom set to the Kenward-Roger method.

Untargeted urine metabolomics analysis

The resulting data matrix was imported into MetaboAnalyst version 5.0⁴⁸ for subsequent univariate/multivariate analysis. To account for repeated measures structure of the data, the focus of our methodology for analysis was to explain the metabolic changes across steers over time and stage. ANOVA-simultaneous component analysis (ASCA) was conducted to identify the effect of time or stage (diet) on the metabolome urine profile of an ADG-classification (least-ADG vs. greater-ADG steers) using a significance level of P < 0.05. Permutation testing (1000 times) was performed to minimize the possibility that observed effects were by chance. Leverage was used to evaluate the importance of the metabolite to the model, and square prediction error (SPE) was used as a test of the fitness of the model for the particular metabolite. Variables with low SPE and higher leverage had a significant contribution to the model and were selected as influentially affected compounds. Two-way ANOVA univariate analysis was conducted to evaluate the effect of breed type on metabolomic profiles of selected ADG-groups. Post hoc analysis with false discovery rate (FDR) correction was performed on both main effects level and interaction of the metabolites identified. Hierarchical cluster analyses were used to group metabolite by ADG-classification time-dependent changes (FDR < 0.05). The interactive principal components analysis (iPCA) was conducted in both two-way-ANOVO and ASCA to identify and visualize stage and breed type effects of least-ADG vs. greater-ADG steers. Pathway analysis was performed using a Bos taurus pathway library, which integrate global pathway enrichment analysis and relative between centrality pathway topology analysis from MetaboAnalyst 3.0 software. The identification and visualization of the top altered pathways were based on KEGG (http://www.genome.jp/kegg/) database sources⁴⁹. Receiver operating characteristic curve (ROC) analysis was performed to evaluate the minimum number of untargeted metabolites that have the sensitivity and specificity to distinguish between ADG-groups using a linear support vector machine (SVM) algorithm. These models were assessed by area under the curve (AUC) and prediction accuracies by cross-validation. The discriminating ability of AUC were considered as follow: 0.9–1.0 = excellent; 0.8–0.9 = good; 0.7–0.8 = fair. A ROC curve plot the true positive rate (Sensitivity) as a function of the false positive rate (100-Specificity) for cut-off points of selected metabolites calculated with a 95% confidence interval. The significance of the model was further validated using permutation test.

Targeted urine bile acids and steroids

BA and steroid concentrations were log-transformed to achieve normality. The generalized linear mixed model (GLIMMIX) procedure of SAS (SAS Institute, Cary, North Carolina) were used to examine the bile acid and steroid responses to experimental variables. The model included time, stage, breed type, ADG-classification and interactions as fixed effect, residual of steers nested by stage as a random effect, with time as the repeated measure, and spatial power covariance as structure of sampling interval.

Pearson’s correlations were performed to identify the associations between urine bile acids and steroid concentrations with performance and carcass traits. Correlation strength was considered strong (|r|≥ 0.7), moderate (0.7 >|r|≥ 0.4) or weak (0.4 >|r|≥ 0.1) (15). A predictive equation for ADG-classification using bile acids and steroids log transform concentrations was performed by stepwise nominal logistic regression analysis with a JMP Pro version 15.0 (SAS Institute, Cary, North Carolina). Basically, for each ADG-category, multiple regression is run to predict the odds ratio into that category. The odds ratios are then converted into probabilities, and the most likely class is chosen. For ADG-prediction model building, data were divided into a 75/25 training:validation set splits to avoid over-fitting and minimize model reliance on sample selection. Validation set was stratified by ADG-classification to balance training set representation and used to test model fit probabilities. The model sensitivity (true positive), specificity (true negative) and accuracy (proportion of true results, either true positive or true negative) was evaluated from confusion matrix based on the final model training set randomly chosen.

An overview of the experimental design and data analysis is presented in the Fig. 8.