The cultured primary myoblasts from the fetal beef longissimus dorsi muscle, the transcriptome sequencing datasets after quality control and alignment, and the differential expression analysis results were achieved from our previous studies (31, 32).

All the animal procedures were carried out according to the protocols approved by the Institutional Animal Care and Use Committee (IACUC) of the College of Animal Science and Technology, Northwest A&F University, China. Ninety-day-old fetal cattle were collected from Tumen slaughterhouse in Xi'an, Shaanxi Province. First, we used 75% alcohol and 1% double-antibody sterile phosphate buffer solution (phosphate buffered saline, PBS) to gently wash the epidermis of the fetal cow to eliminate blood stains and bacterial contamination in an extra-large plate. Second, we cut the muscle tissue pieces and placed them in a 50-ml centrifuge tube, with collagenase I digestion in Dulbecco's modified Eagle's medium (DMEM) at 37°C for 1.5 h. Third, the suspension was filtered and centrifuged, and the supernatant was removed. Fourth, we added four times the volume of 0.25% trypsin in the sediment at 37°C for 30 min. Fifth, the digested sample was filtered through 1-mm stainless steel mesh and 100-μm mesh. Sixth, we added 500 μl of a medium containing 15% fetal bovine serum (FBS) to the pellet to terminate the digestion, and the cells were inoculated in a 6-cm Petri dish. After 20-min culturing, we drew the upper culture medium and continued the culturing process in a new 6-cm Petri dish. When the cell density reached 80–90%, we used them for the subsequent experiments.

RNA for each cell sample was isolated for mRNA sequencing to generate 125 or 150-bp paired-end reads and for miRNA sequencing to generate 50-bp single-end reads on an Illumina Hiseq platform (Illumina, USA). We removed the unqualified reads by quality control to achieve clean reads by in-house perl scripts that were used in our previous studies (31, 32). Then, the miRNA tags were mapped to the reference genome of Bos taurus (UMD_3.1.1/bosTau8) by Bowtie software (version 0.12.9) (33). The mapped miRNA tags were used to seek the known miRNAs using miRBase20.0 as the reference, so the potential miRNA and the secondary structures were obtained by miRDeep2 software (version 2.0.0.5) (34).

The expected number of fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) of each gene was calculated to achieve the average FPKM values for the replicates. The threshold of gene expression was set for mRNA, when the FPKM value is larger than 1. Following the normalization formula (35), miRNA expression levels were also estimated by transcript per million (TPM). A differential expression analysis of two groups (case and control) for both mRNA and miRNA was performed using the R package DESeq (version 1.18.0) (36). The P-values were adjusted using the Benjamini and Hochberg's method for controlling the false discovery rate (FDR). Differentially expressed mRNAs and DE miRNAs were defined when the adjusted P-values were < 0.05. In addition, we calculated fold changes (FCs) between the case and control groups based on the averaged FPKM values and TPM values to define the up-regulated (log₂FC > 0) and down-regulated (log₂FC < 0) mRNAs and miRNAs, respectively.

We used the skeletal muscle cells under the polyphenol RSV treatment as the case group, while the skeletal muscle cells without RSV treatment were considered as the control group. Meanwhile, both the case and control groups had three independent experiments separately for the skeletal muscle cell collections. The correlation coefficients among samples were visualized in the heatmaps using log₁₀(FPKM + 1) for mRNA and log₁₀(TPM + 1) for miRNA in Figure 1.

The R package WGCNA (37) was used to construct the co-expression network. It constructs a similarity matrix by Pearson correlation coefficients to measure the similarity between the gene expression profiles and then transforms the similarity matrix into an adjacency matrix (A) raised to a β exponent (soft threshold) based on the free-scale topology model. In this study, a total of 18,329 genes were filtered from 26,332 genes in mRNA data based on the median absolute deviation (MAD) of each gene bigger than 0.01. The β power parameter (soft threshold) was equal to 12 when the R² of the free-scale topology was equal to 0.8 (Figure 2A). In the miRNA data, 650 miRNAs were filtered from 765 genes, and the β power parameter (soft threshold) was equal to 4 (Figure 2B).

We chose the soft threshold power (β = 12 for mRNA and β = 4 for miRNA) based on the criterion of approximate scale-free topology to construct a weighted gene network and detect the consensus modules with the topological overlap matrix (TOM). The minimum module size was set at 30 for both miRNA data and mRNA data, and the maximum module sizes were set at 18,329 and 650 for mRNA data and miRNA data, respectively. Based on the dissimilarity between module eigengenes (MEs), the modules can be merged, where the first principal component of each module represents the gene expression profiles within the modules (38). Here, we set the cut height for module merging at 0.25, so the modules whose eigengenes are correlated above 0.75 will be merged.

A module association analysis was conducted between the ME and the RSV treatment status (i.e., 0, 0, 0, 1, 1, 1 for three control and three case groups, respectively) to calculate the correlations for the relevant module identifications. We calculated the module significance (MS) [i.e., the average absolute gene significance (GS) of all the genes involved in the module, where GS is measured as log P-value in the linear regression between gene expression and RSV treatment status] to evaluate the correlation strengths. Normally, the module with the highest MS score is the key module (37). Module significance genes in the association analysis (P-value < 0.1) were assigned for functional enrichment analysis. The hub genes were defined as the TOM values up to 0.8.

R package clusterProfiler (version 3.6) (39) was used to test the Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichments. Significant enrichment was defined with the adjusted P-value < 0.1 for both GO terms and pathways.

The bovine genomic sequence (release-99) and the gene annotation file were downloaded from the Ensembl FTP site (http://www.ensembl.org/index.html). We used TBtools to obtain the 3′ UTR sequence of the bovine genomic transcripts (40). Then, the transcript stable IDs were converted to the Ensembl stable IDs using the BioMart website (http://www.ensembl.org/biomart). Accordingly, we found the 3′ UTR sequence of the genes in the turquoise modules. The binding capability of miRNAs and their target genes in the turquoise modules was assessed by RNAhybrid (version 2.1.2) (41), with the minimal free energy hybridization under −20 and the helix constraint from 8 to 12.

From our previous study (31, 32), a total of 3,856 DE mRNAs were identified from 18,329 mRNAs based on the threshold of adjusted P-value < 0.05; meanwhile, 93 DE miRNAs were also identified from 650 miRNAs based on the same thresholds (Table 1). The details of log₂FC, P-value, adjusted P-value, and DE mRNAs with FPKM and miRNAs with TPM of each sample are listed in the Supplemental File 1.

Using 18,329 mRNA and 650 miRNA data for the sample clustering, we found that the samples with RSV treatment were clustered together in mRNA analysis (Figure 3A), while RSV-treated samples were not clustered together by miRNA data (Figure 3B). Generally, 18,329 mRNAs were grouped into 32 modules that had similar co-expressions using the average linkage hierarchical clustering algorithm (Figures 3C,E), where 8,311 mRNAs were grouped into turquoise module as the key module, followed by 2,210 mRNAs into blue module, etc. (Figure 3E; Table 2). However, miRNAs were only grouped into eight modules (Figures 3D,F; Table 2), where 285 miRNAs were grouped into turquoise modules as the key module, followed by 100 miRNAs into blue modules. The mRNAs and miRNAs that were not assigned to any modules were grouped into gray modules (Figures 3C–F).

The eigengene adjacency heatmap indicated that these modules of mRNAs and miRNAs could be clustered further together into groups (Figure 4). After incorporating the RSV treatment trait, we found that the treatment status was clustered with blue and turquoise modules for mRNAs (Figure 4A) and with turquoise modules in a single cluster for miRNAs (Figure 4B).

The module–trait relationship results revealed that RSV treatment had high correlations with turquoise module (0.91, P-value = 0.01), blue module (0.93, P-value < 0.01), and tan module (−0.81, P-value = 0.05) using mRNA data (Figure 5A), whereas only turquoise module (0.96, P-value < 0.01) was highly correlated with treatment status using miRNA data (Figure 5B). Therefore, the turquoise module showed a strongly positive relationship with RSV treatment no matter whether mRNA data or miRNA data is used.

We also used the 3,856 DE mRNAs and 93 DE miRNAs for the weighted gene network construction based on the same soft threshold power (β = 12 for mRNA and β = 4 for miRNA) and then visualized them in the TOM clusters (Figure 6). Here, the minimum module size and maximum block size were set at 30 and 3,856, respectively, for mRNA and were set at 30 and 93, respectively, for miRNA. Only two cluster modules were displayed for both mRNAs and miRNAs, i.e., turquoise and blue modules. In the turquoise module, 2,579 DE mRNAs and 59 DE miRNAs were found, while 1,277 DE mRNAs and 34 DE miRNAs were found in the blue module. Furthermore, the network heatmap of DE mRNAs and DE miRNAs showed a high level of overlap densities among the two clusters (Figure 6).

Based on the 2,579 DE mRNAs in the turquoise module, we performed an enrichment analysis to reveal significant GO terms and KEGG pathways. We found that the three most significant GO terms were actin filament-based process (GO:0030029, adjusted P-value = 1.86 × 10⁻⁶) with 37 DE genes that were enriched in, followed by actin cytoskeleton organization (GO:0030036, adjusted P-value = 4.33 × 10⁻⁶) with 34 DE genes, and actin filament organization (GO:0007015, adjusted P-value = 6.86 × 10⁻⁵) with 24 DE genes in the down-regulated category (Figure 7A). Similarly, the three most significant pathways were regulation of actin cytoskeleton (bta04810, adjusted P-value = 1.40 × 10⁻¹¹) with 47 genes, tight junction (bta04530, adjusted P-value = 2.32 × 10⁻⁶) with 33 genes, and axon guidance (bta04360, adjusted P-value = 2.32 × 10⁻⁶) with 33 genes (Figure 7B).

We selected and constructed the network of four genes within the turquoise module including two significantly up-regulated genes, i.e., sushi domain containing 4 (SUSD4) gene and diacylglycerol O-acyltransferase 1 (DGAT1) gene, and two significantly down-regulated genes, i.e., fibroblast growth factor 18 (FGF18) gene and actin gamma 1 (ACTG1) gene (Figure 8). SUSD4 has associations with 186 genes including 8 up-regulated genes and 176 down-regulated genes (Figure 8A). DGAT1 was closely connected with 70 genes where 32 genes were up-regulated (e.g., proliferating cell nuclear antigen, PCNA) and 38 genes were down-regulated (e.g., phosphatase and tensin homolog, PTEN) (Figure 8B). FGF18 and ACTG1 were negatively connected with 36 and 51 up-regulated genes, respectively (Figures 8C,D). The selected hub gene ACTG1 was negatively related with skeletal muscle myosin heavy chain 3 (MYH3) gene, fibroblast growth factor 2 (FGF2) gene, and uncoupling protein-2 (UCP2) (Figure 8D).

Based on the top significant RSV-induced up- and down-regulated DE miRNAs in the turquoise module, we predicted their target genes. The top four up-regulated miRNAs were bta-miR-34c, bta-miR-432, bta-miR-2344, and bta-miR-154c that targeted 21, 78, 22, and 49 down-regulated genes in the turquoise module, respectively (Figure 9A). Likewise, the top four down-regulated miRNAs, i.e., bta-miR-2310, bta-miR-452, bta-miR-1814c, and bta-miR-199b targeted 59, 62, 58, and 15 up-regulated genes in the turquoise module, respectively, where bta-miR-2310 and bta-miR-1814c targeted the same genes (n = 57) (Figure 9B). Three up-regulated and three down-regulated miRNAs targeted the top up-regulated CDKN1A [adjusted P-value = 2.50 × 10⁻¹⁰⁴ and log₂(FC) = 1.97], while five up-regulated and one down-regulated miRNAs targeted the top down-regulated KCNK12 [adjusted P-value = 5.92 × 10⁻⁶³ and log₂(FC) = −2.32] (Figure 9C).