The BRV strain G8P[7] (RVA/Cow-tc/CHN/10927/2019/G8P[7]) was previously isolated from Chinese diarrheal calves in MA-104 cells and characterized by mass spectrometry. The complete protein sequences of the 11 segments were deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the dataset identifier PXD017047 (Elkady et al., 2021). The accession numbers in the Uniprot database are provided in Supplementary Table 1. The viral titer calculated using the Reed and Muench equation (Reed and Muench, 1938) was 10^4.5 TCID₅₀/50 μl. The MA-104 cell line was purchased from Procell Life Sciences and Technology Co., Ltd. (Wuhan, China). The cells were sub-cultured in six-well tissue culture plates (cat. no. 002019; Corning, Inc., Corning, NY, United States). Dulbecco’s modified Eagle’s media (DMEM) (cat. no. 11965092; Gibco™, Waltham, MA, United States) complemented with 5%–10% fetal bovine serum (cat. no. 26140-079; Gibco™) and antibiotics (1,000 IU penicillin and 1,000 μg/ml streptomycin) (cat. no. 15140122; Sigma-Aldrich, St. Louis, MO, United States) were used.

Cells were used for BRV infection when grown to 80–90% confluence. They were infected with 100 TCID₅₀ of the BRV strain G8P[7]. Uninfected cells were used as the mock-infected group. The cytopathic effect (CPE) of the infected cells was observed at 0, 6, 12, 24, 36, 48, and 72 h post-infection (hpi) in comparison with uninfected cells. Virus propagation and the associated CPE were confirmed by RT-PCR of the VP6 gene (primers: F, 5′-ACCACCAAATATGACACCAGC-3′; R, 5′-CATGCTTCTAATGGAAGCCAC-3′) (Elkady et al., 2021). The PCR products were sequenced and confirmed using an automated sequencer (ABI3730, Applied Biosystems, Bedford, MA, United States; TSINGKE Biological Technology, Beijing, China).

To further detect and confirm the propagation of BRV in the culture fluid, the culture flask that showed typical CPE at 72 hpi was completely frozen at −80°C and then completely thawed in one cycle at 4°C. Thereafter, total RNA was extracted from an almost 4-ml infected cellular fluid using the TRIzol reagent (15596-026; Ambion^®, Life Technologies, Carlsbad, CA, United States) according to the manufacturer’s instruction. The extracted RNA was loaded into an SDS–PAGE 2 × loading buffer. The loaded RNA was screened by 12% resolving polyacrylamide gel using a vertical electrophoresis apparatus and at 80 V, 3 h, 20 mA (two glass plates of 1 mm thickness) (Bio-Rad, Hercules, CA, United States). Finally, the 11 segments of viral double-stranded RNA (dsRNA) were visualized with a silver stain and photographed (Gray and Desselberger, 2000).

The infected and control cell samples (three biological replicates) were collected from six-well plates for each of the following time points: 0 hpi (just after 90 min adsorption) and at 6, 12, 24, 36, and 48 hpi. Three biological replicates of the virus- and mock-infected cultures were prepared at each time point. Total RNA was purified from the cells using the TRIzol reagent (15596-026; Ambion^®, Life Technology). Then, strict quality control was carried out for the RNA samples. The quality control criteria included the following aspects: (1) agarose gel electrophoresis analyzed the integrity of sample RNA and excluded DNA contamination; (2) NanoDrop was used for preliminary quantification, detection of RNA concentration, and purity (OD_260/280 and OD_260/230); and (3) Agilent 2100 Bioanalyzer was used for accurate detection of RNA integrity.

The small RNA-seq process mainly consists of two parts: database-building sequencing and bioinformatics analysis. After detection of qualified samples, a Small RNA Sample Pre-Kit was used to construct the sRNA library. The special structure of the 3′ and 5′ ends of sRNA was used (the 5′ end had the complete phosphate group and the 3′ end had the hydroxyl group). Total RNA was used as the starting sample, and both ends of sRNA were directly added to the joint. It was then reversely transcribed into complementary DNA (cDNA). After PCR amplification, the target DNA fragments were separated by SDS–PAGE, and the cDNA library was recovered by gel cutting. After library construction, Qubit 2.0 was used for preliminary quantification, and the library was diluted to 1 ng/μl. Then, the highly sensitive Agilent 2100 was used to detect the insert size of the library. After meeting the insert size expectations, the effective concentration of the library was accurately quantified by quantitative PCR (qPCR) (effective library concentration, > 2 nM) to ensure quality. After qualified library inspection, different libraries were pooled into Illumina SE50 sequencing according to the requirements of effective concentration and target offline data amount. The reaction system was simultaneously added with DNA polymerase, coupling primers, and four dNTPs with base-specific fluorescence markers (as Sanger sequencing method). These dNTP 3′-OH groups were chemically protected so that only one dNTP can be added. After the addition of dNTP to the synthetic chain, all unused free dNTP and DNA polymerase were eluted. Then, the buffer was added for fluorescence excitation, the fluorescence signal was stimulated by laser, and it was recorded by an optical equipment. Finally, the optical signal was converted into the sequencing base by computer analysis. After fluorescence signal recording, chemical reagents were added to quench the fluorescence signal and remove the dNTP 3′-OH protective group, so that the next round of sequencing reaction was carried out. The Illumina technology, adding only one dNTP at a time, solves the issue of accurately measuring the length of the homomer. All sRNA families in the samples were sequenced and expressed quantitatively with a high-throughput sequencing technology, and the miRNA, small interfering RNA (siRNA), Piwi-interacting RNA (piRNA), and other non-coding RNAs, as well as the corresponding target sequences, were analyzed. The parametric analysis process of animal sRNA was adopted [quality control of the original sequencing data, statistics of clean read length distribution, sRNA annotation, known miRNA annotation, known piRNA annotation, ncRNA analysis, prediction of new miRNA (with reference genome), differential gene analysis, miRNA family analysis, cluster analysis of the DE miRNAs, prediction of miRNA target genes, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the target genes] (Figure 1). The raw sequence data obtained by high-throughput sequencing were analyzed by base calling and transformed into sequenced reads, called raw data or raw reads. The results were stored in FASTQ format, which contains the sequence information of the sequencing sequences (reads) and their corresponding sequencing quality information. The clean reads of each sample were screened, sRNAs within a certain length were selected for subsequent analysis, and the length distribution of sRNAs was determined; the length range of animal sRNA is 18–35 nt (Bartel, 2018). The sRNAs were aligned and analyzed in comparison with the reference genome sequence in the Sanger miRBase 21.0 database¹ (Langmead et al., 2009). The expression levels of the known and novel miRNAs in each sample were calculated, and TPM (transcript per million) was used to normalize the expression levels (Zhou et al., 2010). For the samples with biological duplication, DESeq2 (Love et al., 2014) based on negative binomial distribution was used for analysis. For samples without biological duplication, DEGseq (Wang et al., 2010) based on the TMM (trimmed mean of M-values) algorithm was used to standardize the read count data and then conduct difference analysis. The thresholds for differences of miRNA filter condition are: p_adj < 0.05 and | log2(fold change) | > 1. Additionally, the DE miRNAs between the different BRV infection time points were identified. Moreover, animal miRNA target genes were predicted using miRanda² (Enright et al., 2003) and RNAhybrid³ (Krüger and Rehmsmeier, 2006).

After obtaining the DE miRNAs from each infection time point, we conducted GO and KEGG enrichment analyses for the collected DE miRNA target genes in each group according to the corresponding relationship between miRNA and its target genes. GO⁴, which is the international standard classification system for gene function, was used to study the function of the target mRNA. The GO enrichment analysis method used was GOseq (Young et al., 2010), based on Wallenius non-central hypergeometric distribution. In addition, KEGG (Kanehisa et al., 2008) was adopted to analyze the biochemical metabolic pathways and signal transduction pathways (protein--protein interaction) involving the target gene candidates; pathway significance enrichment analysis was conducted using KEGG PATHWAY (KOBAS 2.0 annotation tool)⁵. The hypergeometric test/Fisher’s exact test was used to determine which pathway was significantly enriched in the candidate target genes compared with the whole genome background. The calculation formula for this analysis is as follows:

where N is the number of genes with pathway annotation in all genes, n is the number of candidate target genes in N, M is the number of genes annotated as a specific pathway in all genes, and m is the number of candidate target genes annotated as a specific pathway. For the false discovery rate (FDR) correction method, Benjamini–Hochberg (BH) was used for p-Value correction; the smaller the corrected p-Value, the more significant it is. Pathways with values less than 0.05 were defined as significantly enriched in candidate target genes. KEGG enrichment was measured using Rich factor, Q-Value, and the number of genes enriched in a pathway. Rich factor refers to the ratio of the number of genes in the DE genes to the total number of genes in the pathway item in all the annotated genes. The larger the Rich factor, the greater the degree of enrichment. The Q-value is the p-Value after multiple hypothesis testing and correction. The range of Q-values is [0, 1]; the closer the value is to zero, the more significant the enrichment.

We selected some DE miRNAs (n = 23) and their target genes (n = 28) for validation by real-time relative qPCR based on the following criteria: (1) they were previously reported to modulate the cellular machinery (apoptosis, autophagy, cancer, and immunity) with their corresponding target mRNAs of immune importance; (2) they were in the top 20 significantly enriched cellular signaling pathways; (3) some novel DE miRNAs with the highest differential folds; and (4) some of the corresponding target mRNAs associated with pathways responsible for host antimicrobial response. Total RNAs at different time points—0 hpi (just after 90 min adsorption) and 6, 12, 24, 36, and 48 hpi—from the virus-infected and non-infected groups were extracted and purified. They were then reverse transcribed to cDNA using miRNA 1st-Strand cDNA Synthesis Kit (by stem–loop) (#MR101; Vazyme, Nanjing, China), according to the manufacturer’s instruction, using the RT-specific primer. The gene expression and the differential expression of the miRNAs were determined with the miRNA Universal SYBR^® qPCR Master Mix (#MQ101; Vazyme), following the manufacturer’s instruction, using specific reverse and forward primers. U6 was used as an internal reference gene. In addition, for the detection of the DE target mRNAs, RNA was reverse transcribed into single-strand cDNA using the Transcriptor First Strand cDNA Synthesis Kit (cat. no. 04897030001; Roche, Basel, Switzerland) according to the manufacturer’s instructions. The gene expression and the differential expression of the target mRNAs were determined using the AceQ^® qPCR SYBR Green Master Mix (#Q111; Vazyme). β-actin was used as an internal reference gene. The primers used in the reverse transcription and real-time qPCR reactions are listed in Supplementary Table 2. All qPCR reactions were performed using Bio-Rad CFX-96. The relative normalized expressions of the DE miRNAs and their target mRNAs were calculated with the 2^–ΔΔCT method using Bio-Rad software data analysis from three independent experiments. The gene regulation threshold > 2 and an FDR-adjusted p-Value < 0.05 were considered as the parameters for statistical significance of the DE genes.

Unpaired, two-tailed Student’s t-test for grouped data was used with GraphPad Prism 8 software, and comparative analysis of the fold change (FC) values (2^–ΔΔCT) between the non-infected and virus-infected groups was performed. All qPCR data are presented as the mean ± SD and p < 0.05 was considered statistically significant. In addition, sequencing data above the thresholds were screened using default statistical analysis, which was explained in the methodological description.

Compared with the mock-infected (control) cells, no CPE was observed in the BRV-infected MA-104 cells from 0 to 36 hpi (Figure 2A), whereas the typical CPE of BRV appeared within 48 and 72 hpi (Figures 2B,C). To confirm viral replication in MA-104 cells following BRV infection, the cell-associated viruses were harvested before and after the occurrence of CPE and were confirmed via RT-PCR to amplify the viral (VP6) gene (accession. no MW018712) (Figure 3A). Sequencing analysis of the PCR products further validated BRV infection. Additionally, it was confirmed by the typical RNA pattern of this group A BRV strain, G8P[7], in comparison with the negative control using SDS–PAGE and silver stain (Figures 3B,C).

RNAs of grade A quality were used for sequencing and differential analysis of the virus-treated and non-infected MA-104 cells in triplicate. Raw reads from the libraries derived from the virus-treated and control cell groups at different time points after infection (0–48 hpi) were generated. Then, the clean reads were obtained after filtration and sRNAs within a certain length were selected for subsequent analysis, as shown in Supplementary Tables 3a,b. The obtained length of the clean reads was 22–24 nt, mainly with the length of 23 nt at 24 hpi (Figure 4A) and at 0, 6, 12, 36, and 48 hpi (Supplementary Figures 1a–e). The percentages of non-coding nucleic acids [e.g., known miRNAs, novel miRNAs, ribosomal RNA (rRNA), transfer RNA (tRNA), snRNA, small nucleolar RNA (snoRNA), repeats, exons, introns, etc.] in the virus-infected and control groups were shown at 24 hpi (Figure 4B) and at 0, 6, 12, 36, and 48 hpi (Supplementary Figures 2a–e). The mapped known and novel miRNAs of Macaca mulatta species are listed in Supplementary Table 4. After BLASTN analysis by miRBase 20.0 with M. mulatta miRNAs, 579 known and 144 novel mature miRNAs were obtained. The DE miRNAs with true significance (FDR < 0.05 and log2 FC > 1) at the different time points were identified. Basically, the DE miRNAs, in both number and genes, varied with time post-infection, including 20 upregulated and 10 downregulated at 0 hpi, 11 upregulated and 2 downregulated at 6 hpi, 8 upregulated and 14 downregulated at 12 hpi, 8 upregulated and 5 downregulated at 24 hpi, 27 upregulated and 23 downregulated at 36 hpi, and 21 upregulated and 11 downregulated at 48 hpi (Supplementary Tables 5a–f). However, some common DE miRNAs exist between the different time points of BRV infection, such as mml-miR-99a-5p between 0 and 6 hpi, mml-miR-342-5p between 12 and 24 hpi, novel_458 between 12 and 36 hpi, mml-miR-411-3p between 12 and 48 hpi, novel_356 and novel_66 between 24 and 36 hpi, and mml-miR-411-5p between 36 and 48 hpi (shown in Figure 4C for 24 hpi and in Supplementary Figures 3a–e for 0, 6, 12, 36, and 48 hpi).

All target genes and their specific mRNA transcripts in the known and novel miRNAs were bioinformatically expected, and the corresponding relationship between the DE miRNA and the target genes was characterized (Supplementary Table 6a). All annotated target mRNAs are listed in Supplementary Table 6b. GO and KEGG analyses were conducted to characterize the functions of the predicted target mRNAs of the DE miRNAs. The difference in probability was estimated based on the bias of gene length so that the probability of the GO term being enriched by the candidate target genes can be more accurately calculated. The GO enrichment list of the candidate target genes associated with the different infection times is shown in Supplementary Tables 7a–f. Mostly, the most significantly enriched GO terms of the target mRNAs between the infected and control cell groups were involved in biological process (BP), cellular component (CC), and molecular function (MF) at 0, 6, 12, 24, 36, and 48 hpi (Figures 5A–F). At 0 hpi, the targets included in functions such as membrane-bound organelle and plasma membrane could help the virus to attach to the cellular receptor and internalize the cells; the target mRNAs were mainly associated with the apoptotic process at 6 hpi and contributed to the trans-Golgi membrane at 12 hpi. Evidently, the target genes produced at 24 hpi played a role in polymerase activity and genome replication, such as RNA polymerase I regulatory region, RNA polymerase I core element sequence, and binding proteins. Moreover, the targets were involved in the localization and regulation of developmental process at 36 hpi and in DNA packaging, trans-Golgi network membrane, and polymerase III regulatory regions at 48 hpi. On the other hand, all the enriched pathways characterized through KEGG mapping were shown (Supplementary Tables 8a–f), highlighting the top 20 significantly enriched pathways represented by the predicted target mRNAs of each comparison group (Figures 6A–F and Supplementary Tables 9a–f). KEGG enrichment elucidated the pathways that played a part in the pathogenesis of the virus, such as the B-cell receptor signaling pathway, Fc gamma R-mediated phagocytosis, and chemokine signaling pathway at 0 hpi; pathways in cancer, Rap1 signaling pathway, mTOR signaling pathway, Wnt signaling pathway, lysosome, leukocyte trans-endothelial migration, lysine degradation, and Fc gamma R-mediated phagocytosis at 6 hpi; and pathways in cancer, MAPK signaling, primary immunodeficiency, acute myeloid leukemia, cell adhesion molecules (CAMs), long-term potentiation, and Fc gamma R-mediated phagocytosis at 12 hpi. In addition, at 24 hpi, the cAMP, NF-κB, and TNF signaling pathways, Fc gamma R-mediated phagocytosis, B-cell receptor signaling pathway, and chronic myeloid leukemia were involved. On the other hand, the MAPK and Ras signaling pathways, Fc gamma R-mediated phagocytosis, long-term potentiation, and B-cell receptor signaling pathway were included at 36 hpi. Last but not least, pathways in cancer, MAPK signaling, the NF-κB and B-cell receptor signaling pathway, endocytosis, AMPK signaling pathway tight junction, arachidonic acid metabolism, and the chemokine signaling pathway were implicated at 48 hpi.

We validated the DE miRNAs (n = 23) and their corresponding target mRNAs (n = 28) at 0–48 hpi. The expression change values of all validated miRNAs and genes were in agreement with those in the deep sequencing and the bioinformatics analyses. Some of the DE miRNAs showed significant upregulation, such as mml-let-7f-5p (0 hpi; FC = ∼2.7, p = ∼0.00004) and its target (SERINC4; FC = ∼45, p = ∼0.001871), mml-miR-486-3p (6 hpi; FC = ∼11.5, p = ∼0.000011) and its target (PDPK1; FC = ∼42, p = ∼0.02), and mml-miR-486-5p (at 6 hpi; FC = ∼18.8, p = ∼0.001). Interestingly, mml-miR-127-3p (6 hpi; FC = ∼242.5, p = ∼0.000002) and its target (PARD6G; FC ∼ −2.3, p = ∼0.02) showed very high upregulation. In addition, mml-miR-132-3p (12 hpi; FC = ∼4.8, p = ∼0.03666) and its target (AZGP1) were dramatically upregulated (FC = ∼19.34, p = ∼0.0068). Moreover, mml-miR-369-5p (12 hpi) was upregulated (FC = ∼4.2, p = ∼0.003), and its target (CD4) was sharply downregulated with (FC ∼ −15.3, p = ∼0.000719). Additionally, mml-miR-451 (24 hpi; FC = ∼2.06, p = ∼0.040658) and its target ERAL1 (FC = ∼3.9, p = ∼0.01) were upregulated, and novel_352 (24 hpi) showed FC = ∼2 (p = ∼0.05). Meanwhile, mml-miR-21-3p (36 hpi) exhibited FC = ∼2.56 (p = ∼0.000476), and its target (COG6) was upregulated (FC = ∼6.29, p = ∼0.001113). Also, mml-miR-106b-3p (36 hpi; FC = ∼2.7, p = ∼0.0008) and its target (BDNF; FC ∼ −13, p = ∼0.007) and mml-miR-299-3p (48 hpi; FC = ∼14, p = ∼0.0003) and its target (LAMB2; Fc = ∼65.7, p = ∼0.000097) were dramatically upregulated, while novel_366 (48 hpi) and its target (BCL2A1) showed values of FC = ∼8.3 (p = ∼0.02) and FC = ∼27.5 (p = ∼0.003), respectively. mml-miR-197-3p (24 hpi; FC ∼ −2.01, p = ∼0.005337) and its two targets, TSC2 (FC ∼ −2.05, p = ∼0.04) and TRBV16 (FC ∼ −2.39, p = ∼0.05), exhibited downregulation; however, its other target (ORMDL3) was dramatically upregulated (FC = ∼5.7, p = ∼0.003546). Moreover, mml-miR-204-3p (24 hpi) showed FC ∼ −2.47 (p = ∼0.012433), and most of its targets were upregulated: IFITM5 (FC = ∼ = 6.9, p = ∼0.027985, TRIM17 (FC = ∼11.7, p = ∼0.000946), mTORC1 (FC = ∼15.5, p = ∼0.002), OASL (FC = ∼8.9, p = ∼0.000035), MAP3K12 (FC = ∼5.5, p = ∼0.000283), ISG20 (FC = ∼9.8, p = ∼0.003473), PHPT1 (FC = ∼31, p = ∼0.000043), IKBKB (FC = ∼2.12, p = ∼0.000007), and AKT1 (FC = ∼7.2, p = ∼0.000028). However, its target BCL6 was downregulated (FC ∼ −2.5, p = ∼0.02). mml-miR-504-5p (24 hpi) was also downregulated (FC ∼ −6.2, p = ∼0.000007), and its two targets were upregulated: STAT2 (FC = ∼32.9, p = ∼0.000005) and ISG15 (FC = ∼4.8, p = ∼0.006). Additionally, at 36 hpi, mml-miR-365-1-5p showed a downregulation (FC ∼ −2.9, p = ∼0.004), along with its targets: HSPB6; FC ∼ −2.3, p = ∼0.0002) and ACACA (FC ∼ −2.08, p = ∼0.000263). Tables 1a,b list the DE miRNAs whose targets were implicated in the signaling pathways, while the DE miRNAs without targets or without significant targets and not implicated in the signaling pathways, confirmed by qPCR [marked with (−)] are shown in Figures 7A–F, 8A–F and in Supplementary Table 10. Additionally, the expression levels of the common DE miRNAs between the different time points of infection were analyzed and confirmed. The miRNA mml-miR-99a-5p was upregulated at 0 and 6 hpi (FC = ∼6.11, p = ∼0.011099 and FC = ∼7.79, p = ∼0.001, respectively). At 36 hpi, mml-miR-411-5p had FC = ∼3.5, p = ∼0.003, and it showed an almost the same expression level at 48 hpi (FC = ∼3.3, p = ∼0.000735). Similarly, the expression pattern of mml-miR-411-3p at 12 hpi (FC = ∼2.03, p = ∼0.04) was close to that at 48 hpi (FC = ∼2.4, p = ∼0.002). mml-miR-342-5p had FC ∼ −1.4 (p = ∼0.2) at 12 hpi; however, at 24 hpi, it showed a surge in fold change, with FC ∼ −16.2 (p = ∼0.000002), with its target MAPK15 (12 hpi: FC = ∼61.5, p = ∼0.009; 24 hpi: FC = ∼14.7, p = ∼0.000222). Besides, novel_66 experienced very high regulation at 24 hpi (FC = ∼13.4, p = ∼0.0009), which dramatically dropped at 36 hpi (no fold change and p-Value). Moreover, novel_458 (12 hpi: FC = ∼23.7, p = ∼0.000111) showed a decrease to FC = ∼3.1 (p = ∼0.001) at 36 hpi, as well as novel_356 (24 hpi: FC = ∼3.5) decreasing to FC = ∼1.4 at 36 hpi, both without p-Value (Figure 9, Table 2, and Supplementary Table 11).