This study is part of a broader research initiative to evaluate the long-term molecular effects of maternal nutrition on the development and performance of beef cattle offspring. The experimental animals were managed at the Faculty of Animal Science and Food Engineering (FZEA-USP) campus (Pirassununga, SP, Brazil). The Research Ethics Committee of FZEA-USP approved the study under protocol number 1843241117, in accordance with the ethical standards and guidelines of the National Council for the Control of Animal Experimentation.

The nutritional interventions applied to 126 pregnant Nellore cows during gestation were previously described by Schalch Junior et al. (2022). Briefly, cows were artificially inseminated (AI) using semen from four sires, and pregnancy was confirmed 30 days post-AI. The dams were then stratified into three experimental groups (n = 42 per group), balanced for age, body weight (BW), and body condition score. Each group was maintained in grazing paddocks with Urochloa brizantha cv. Marandu, which were equipped with troughs for water and feed supplementation. The cows were assigned to one of three prenatal nutritional treatments: (1) NP (Not Programmed; control) – mineral supplementation only (0.3 g/kg BW per day) throughout gestation; (2) PP (Partially Programmed) – protein–energy supplementation (3 g/kg BW per day) provided only during the third trimester; and, (3) FP (Fully Programmed) – the same protein–energy supplementation (3 g/kg BW per day) provided continuously from pregnancy confirmation until calving. All groups received mineral supplementation at 0.3 g/kg of BW per day; however, in the PP and FP groups, the minerals were incorporated into the protein–energy supplement formulation (Table 1). The quality of the grazing areas was monitored through periodic bromatological analyses to ensure comparable nutrient availability among treatments. The total digestible nutrient (TDN) content averaged 63.07%, 64.1%, and 61.43% for the NP, PP, and FP groups, respectively, while crude protein (CP) levels were 7.38%, 7.82%, and 7.40%, respectively. Further methodological details and the effects of maternal nutritional treatments on dam performance and metabolism are provided in Schalch Junior et al. (2022).

After calving, protein–energy supplementation for the PP and FP dams was discontinued. Thereafter, all groups, including NP, received 0.3 g/kg BW of mineral supplement daily and were managed together in a single grazing system with Urochloa brizantha cv. Marandu. Standardized health and feeding protocols were followed for cow–calf management until weaning at 240 ± 28 days. Post-weaning, calves were separated by sex and reared under uniform management conditions until reaching 570 ± 28 days of age. The young bulls received two types of supplements depending on the season: an energy supplement during the dry season and a protein supplement during the wet season, while continuing to graze on the same pasture. The finishing phase began at 570 ± 28 days of age and included 63 bulls, lasting until slaughter at 676 ± 28 days of age. During this period, the animals received three sequential diets: a 15-day adaptation diet, a subsequent diet for 35 days, and a final diet for 56 days (Polizel et al., 2022a). At the end of the experimental period, all animals were slaughtered at the FZEA-USP slaughterhouse, following the regulations established by the Brazilian Ministry of Agriculture, Livestock, and Food Supply (MAPA; Brasilia, DF, Brazil). Figure 1 illustrates a schematic overview of the experimental design.

Immediately after slaughter, samples of the longissimus thoracis muscle (located between the 12th and 13th ribs) were collected from all bulls for transcriptomic analyses. Tissue samples were obtained within 15 min postmortem, rapidly snap-frozen in liquid nitrogen, and stored at −80 °C until RNA extraction. From the 63 bulls, a subset of 15 animals (n = 5 per maternal nutrition group), all sired by the same bull (to minimize paternal genetic variation), was randomly selected for RNA sequencing (RNA-Seq) analyses.

Total RNA was extracted from 100 mg of longissimus thoracis muscle using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, United States), following the manufacturer’s instructions. RNA concentration and purity were assessed with a DS-11 spectrophotometer (DeNovix, Wilmington, DE, United States), and RNA integrity was evalauted with the Bioanalyzer 2,100 (Agilent Technologies, Santa Clara, CA, United States). Only samples with a RNA Integrity Number (RIN) > 7.0 were retained for subsequent analyses.

For library construction, 0.1–1 µg of high-quality RNA was used following the TruSeq Stranded mRNA Sample Preparation protocol (Illumina, San Diego, CA, United States). Library quantification was performed via quantitative PCR (qPCR) using the KAPA Library Quantification kit (Roche, Basel, Switzerland), and average fragment size was determined with the Bioanalyzer. Sequencing was performed on a single flow-cell lane using the TruSeq PE Cluster kit v3-cBot-HS and a paired-end sequencing configuration on the Illumina HiSeq 2,500 platform. All sequencing procedures were conducted by NGS Soluções Genômicas (Piracicaba, São Paulo, Brazil).

Quality control of the raw RNA-Seq reads was performed using FastQC (v0.12.1), and adapter sequences as well as low-complexity reads were removed with SeqyClean (v1.10.09; Zhbannikov et al., 2017). The cleaned reads were then aligned to the Bos taurus reference genome (ARS-UCD1.3) using the STAR aligner (version 2.7.11b; Dobin et al., 2013) in two-pass mode (Veeneman et al., 2016) based on the Ensembl genome annotation (release 113) database. The resulting coordinate-sorted BAM files (“Aligned.sortedByCoord.out.bam”) were assembled into transcripts using StringTie (v.2.1.7; Pertea et al., 2015) with default parameters. Subsequently, the transcriptomes from all the samples were merged using the StringTie “merge” function to produce a unified transcriptome assembly. By default, StringTie assigns the prefix “MSTRG” to novel genes and transcripts.

To classify known and novel transcripts, the unified GTF file was compared to the Bos taurus reference annotation (ARS-UCD1.3.113) using GffCompare (v.0.12.10; Pertea and Pertea, 2020). Predicted transcripts were categorized according to GffCompare class codes as follows: known isoforms of a reference gene (class code “=”), novel isoforms of a reference gene (class codes “c”, “k”, “j”, “m”, “n”, or “o”), novel loci (class codes “i”, “u”, “y”, or “x”), and, potential artifacts (class codes “e”, “s”, or “p”) (Halstead et al., 2021).

The GffRead tool (v.0.12.7; Pertea and Pertea, 2020) was used to generate a transcriptome FASTA file from the StringTie-assembled GTF file and the Bos taurus reference genome FASTA. This custom transcriptome assembly was then indexed and used for transcript quantification with Salmon (v1.10.1; Patro et al., 2017) operating in quasi-mapping mode to estimate transcript-level abundance.

Transcript-level counts from Salmon were imported into R (v4.4.1) using the tximport package (v.1.32.0; Soneson et al., 2016). A DGEList object was then created using edgeR (v.4.2.2; Robinson et al., 2010). Lowly expressed transcripts were filtered out using the “filterByExpr” function and 39,280 transcripts were kept for subsequent analyses. Data normalization was performed using the Trimmed Mean of M-values (TMM) method based on library sizes. Differential expression analyses were conducted in edgeR using the exact test based on the negative binomial distribution. The analysis included the prenatal nutritional group as the only factor of interest. Statistical significance was assessed using the “exactTest” function, and p-values were adjusted for multiple testing using the false discovery rate (FDR) approach.

The normalized counts were used for differential expression (DE) analysis, whereas log₂-transformed counts per million (CPM) values were used for performing principal component analysis (PCA) using the “prcomp” R function. For gene-level analyses, the same computational workflow was employed, but gene expression matrices were obtained using the “extractGeneExpression” function from the IsoformSwitchAnalyzeR package (v.2.2.0; Vitting-Seerup et al., 2019). After filtering, 15,832 genes were retained for downstream analyses. Transcripts and genes were considered significantly differentially expressed when the FDR-adjusted p-value was lower than 0.05 and the absolute log₂ fold change (|log₂FC|) exceeded 1.

The pipeline of IsoformSwitchAnalyzeR (Vitting-Seerup et al., 2019) was used to investigate differential transcript usage (DTU). Only genes with two or more isoforms (“removeSingleIsoformGenes” = TRUE), an average expression >10 transcripts per million (TPM; “geneExpressionCutoff” = 10) across experimental groups, and at least one isoform expressed >3 TPM (“isoformExpressionCutoff” = 3) in at least one group were retained. After filtering, 7,973 isoforms remained for downstream analyses.

Isoform switching in each contrast (FP × NP, FP × PP, and NP × PP) was assessed using the DEXSeq method (Anders et al., 2012), which implements a generalized linear model (GLM) framework based on the negative binomial distribution to test for differences in relative transcript (or exon) usage between conditions. In this approach, count data from transcript regions are modeled so that changes in relative usage across groups can be detected while accounting for biological variability. DEXSeq fits models with and without condition–region (exon) interactions and performs likelihood ratio tests to identify significant differential usage, with p-values adjusted for multiple testing using the false discovery rate (FDR) method. The exon- or region-level differential usage results are subsequently integrated within the IsoformSwitchAnalyzeR framework to infer transcript-level (isoform) usage patterns. The difference in isoform fraction (dIF) was calculated as IF₂ – IF₁, where IF₁ and IF₂ represent the isoform fractions of the first and second groups in each contrast, respectively (e.g., FP in FP × NP or FP × PP). The dIF value serves as an estimate of effect size, analogous to fold change in differential expression analysis. Isoform switches were considered as significant when |dIF| > 0.1 and FDR <0.05.

Further analyses of the switched isoforms were performed to evaluate their potential consequences on gene function, including intron retention (IR), changes in the amino acid sequences of open reading frames (ORFs), susceptibility to nonsense-mediated decay (NMD), modifications in protein domain structures, and variations in coding potential. A comprehensive examination of notable isoform switches was conducted using multiple bioinformatics tools: IUPred2A web server (Mészáros et al., 2018) for predicting intrinsically disordered regions; SignalP (v.5.0; Teufel et al., 2022); for signal peptide prediction; CPC2 (v.1.0.1; Kang et al., 2017) for coding potential evaluation; Pfam (v.37.0; Finn et al., 2014); for identifying conserved protein domains; DeepTMHMM (v.1.0.44; Hallgren et al., 2022); for protein topology prediction; and DeepLoc2 (v.2.0; Thumuluri et al., 2022); for subcellular location inference. In addition, AS events associated with isoform switching were investigated, including alternative 3′ and 5′ splice sites (A3, A5), alternative transcription start and termination sites (ATSS, ATTS), exon skipping (ES), and intron retention (IR). These events can reveal potential modifications in gene function mediated by altered exon junctions and transcript structural variation.

As the primary objective of this study was to evaluate the effects of maternal nutrition at the transcript level, rather than through conventional gene-level analyses, the over-representation analysis (ORA) was performed using transcript-based data. Genes corresponding to significant transcripts identified in differential transcript expression (DTE) or DTU analyses were used as input for enrichment testing. ORA was conducted separately for each contrast (FP × PP, PP × NP, and FP × NP) and for both analysis type (DTE and DTU). Functional enrichment was carried out using Gene Ontology (GO) analysis including the three main GO categories: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). Furthermore, metabolic pathways analyses were performed using the Reactome (REAC) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases via the g:Profiler platform (Kolberg et al., 2023). Enriched processes were considered significant when FDR was lower than 0.05. The overlap of significant transcripts between DTE and DTU analyses, as well as the intersection of enriched biological processes, was visualized using the R packages “eulerr” (v7.0.2) and “ggvenn” (v0.1.10). A schematic overview of the main bioinformatics workflow is presented in Figure 2.

A total of 27,412 genes were identified in the assembled muscle transcriptome, including 26,080 known genes and 1,332 novel genes (Figure 3A). At the transcript level, 111,185 transcripts were detected, of which 75,284 transcripts were classified as known and 35,901 as novel (Figure 3B). Among the novel transcripts, 32,705 originated from reference genes, 3,196 from novel genes, and 177 were categorized as potential artifacts (Figure 3C). Comprehensive details for each gene and transcript are provided in the Supplementary Material: Additional File 1 (known transcripts), Additional File 2 (novel loci), Additional File 3 (novel transcripts from reference genes), and Additional File 4 (potential artifacts).

The PCA of both gene-level (Figure 4A; PC1 = 14% and PC2 = 12.4%) and transcript-level (Figure 4B; PC1 = 10.7% and PC2 = 9.3%) expression data showed similar distribution patterns, with considerable overlap among the experimental groups. This overlap indicates the lack of a clear global distinction between the maternal nutrition treatments at either gene or transcript level.

Differential expression analysis enabled the identification of 13 DEGs between the PP and NP groups (FDR <0.05 and |log₂FC| > 1), two genes between FP × NP, and one gene between FP × PP (Figure 4C; Additional File 5). At the transcript level (Figure 4D; Additional File 6), 14 transcripts were DE in the FP × NP contrast, 12 in PP × NP, and 10 in FP × PP. Notably, among the 13 DEGs, five corresponded to novel genes (MSTRG.19105, MSTRG.2378, MSTRG.10953, MSTRG.16026, and MSTRG.13630), while at the transcript level, 10 unique novel transcripts were identified (Figure 4).

A total of 120 transcripts exhibited DTU across the three maternal nutrition contrasts: NP × PP, FP × NP, and FP × PP (Figure 5A; Additional File 7). Among these, 30 transcripts were detected in NP × PP comparison, 87 in FP × NP, and three in FP × PP. We further assessed the functional impact of isoform switching for each significant isoform detected in the DTU analysis (Figure 6A; Additional File 8). Although some isoform switching events occurred in the FP × PP comparison, none resulted in detectable functional consequences (Figure 6A). In contrast, 21 switch consequences were detected in the NP × PP comparison, involving different domains, coding potential, ORF similarity, intron retention, and other features across 11 genes. For the FP × NP contrast, 93 switch consequences were detected, affecting 38 unique genes. Representative examples of switch consequences in key genes are illustrated in Figure 5B (SORBS3 gene), Figure 5C (SLC7A8 gene), and Figure 5D (SLC25A30 gene). These genes were selected as illustrative cases because they showed robust and biologically meaningful isoform switching events (high dIF values and significant FDR), and because they are functionally linked to metabolic, transport, and cellular adaptation processes previously reported to be sensitive to maternal nutritional status. As such, they provide representative examples of how prenatal nutrition may influence transcript usage with potential functional consequences in offspring muscle. Figure 6A provides an overview of the number of genes exhibiting isoform switching with associated functional consequences across all pairwise comparisons.

AS events associated with isoform switching are shown in Figure 6B, with detailed information for each affected gene and isoform provided in Additional File 9. In the FP × NP contrast, the most frequent AS event was ATTS in the less-used isoform. For the FP × PP comparison, only two events AS were detected: ATSS in the more-used isoform and ES in the less-used isoform. Finally, in the NP × PP contrast, the predominant AS event was ATSS in the more-used isoform.

In the FP × NP comparison, ORA of the genes identified through the DTE analysis revealed significant enrichment of pathways related to biotin metabolism (GO:MF:0004077 [FDR = 0.0196], GO:MF:0009374 [FDR = 0.0196], and KEGG:00780 [FDR = 0.0226]). Additional enriched pathways were associated with transport of small molecules (REAC:R-BTA-382551 [FDR = 0.0260]), and with cell adhesion and cadherin regulation (REAC:R-BTA-9759476 [FDR = 0.0260], REAC:R-BTA-9764260 [FDR = 0.0260], REAC:R-BTA-9762292 [FDR = 0.0260], and REAC:R-BTA-9759475 [FDR = 0.0260]. The endosomal/vacuolar metabolic pathway was also significantly enriched (REAC:R-BTA-1236977, FDR = 0.0260).

The DTU analyses further revealed significant enrichment of pathways related to intracellular trafficking and vesicle-mediated transport (REAC:R-BTA-1236977 [FDR = 0.00047] and REAC:R-BTA-1236974 [FDR = 0.0021]), as well as cellular structural organization (GO:CC:0043229 [FDR = 0.0024], GO:CC:0031982 [FDR = 0.0024], GO:CC:0043226 [FDR = 0.0024], GO:CC:0110165 [FDR = 0.0035], GO:CC:0005622 [FDR = 0.0038], and GO:CC:0005575 [FDR = 0.0038]). In addition, immune-related pathways were significantly enriched (REAC:R-BTA-983169 [FDR = 0.0033] and REAC:R-BTA-983170 [FDR = 0.0033]).

In the FP × PP contrast, over-represented processes identified from the DTE analysis were predominantly associated with amino acid metabolism and transmembrane transport (GO:BP:0035524 [FDR = 0.0191], GO:BP:1903801 [FDR = 0.0191], GO:BP:1904273 [FDR = 0.0191], GO:BP:1904557 [FDR = 0.0191], GO:CC:0015827 [FDR = 0.0191], and GO:CC:0015829 [FDR = 0.0191]). Additional enriched processes were associated with inositol metabolism and biosynthesis (GO:MF:0004512 [FDR = 0.0109], GO:BP:0006021 [FDR = 0.0191], GO:BP:0006020 [FDR = 0.0191], and GO:MF:0016872 [FDR = 0.0109]).

In the PP × NP contrast, genes identified through DTE analysis were predominantly enriched in pathways related to biotin metabolism (GO:MF:0004077 [FDR = 0.0050], GO:BP:0009374 [FDR = 0.0050], and KEGG:00780 [FDR = 0.0155]). Additional enriched biological processes were associated with inositol metabolism and biosynthesis (GO:MF:0004512 [FDR = 0.0050], GO:MF:0016872 [FDR = 0.0050], and GO:BP:0006021 [FDR = 0.0331]). Regarding DTU-based ORA, significant enrichment was observed for molecular functions related to adiponectin binding and small nuclear ribonucleoprotein (snRNP) binding (GO:MF:0055100 [FDR = 0.028] and GO:MF:0070990 [FDR = 0.028], respectively). A summary of all over-represented biological processes and molecular functions across the three nutritional contrasts and both analytical approaches (DTE and DTU) is presented in Figure 7. Complete details of the ORA are provided in Additional File 10.

Among the significant isoforms identified in both the DTE and DTU analyses, only MSTRG.1066.2, an isoform of the SLC7A8 gene, was found to overlap between them (Figure 8A). The ORA also revealed processes shared across comparisons and analysis types (Figure 8B), with detailed information provided in Table 2. In the FP × NP comparison, the DTU and DTE analyses shared 15 over-represented processes, indicating convergence in transcript-level and isoform-usage responses to maternal nutrition. In contrast, no overlap was observed between DTU and DTE in the PP × NP comparison. As described in Section 3.4, no over-represented processes were identified for the DTU analysis in the FP × PP contrast; therefore, this comparison was excluded from the overlap assessment.