MicroRNA Profile of MA-104 Cell Line Associated With the Pathogenesis of Bovine Rotavirus Strain Circulated in Chinese Calves

Bovine rotavirus (BRV) causes massive economic losses in the livestock industry worldwide. Elucidating the pathogenesis of BRV would help in the development of more effective measures to control BRV infection. The MA-104 cell line is sensitive to BRV and is thereby a convenient tool for determining BRV–host interactions. Thus far, the role of the microRNAs (miRNAs) of MA-104 cells during BRV infection is still ambiguous. We performed Illumina RNA sequencing analysis of the miRNA libraries of BRV-infected and mock-infected MA-104 cells at different time points: at 0 h post-infection (hpi) (just after 90 min of adsorption) and at 6, 12, 24, 36, and 48 hpi. The total clean reads obtained from BRV-infected and uninfected cells were 74,701,041 and 74,184,124, respectively. Based on these, 579 were categorized as known miRNAs and 144 as novel miRNAs. One hundred and sixty differentially expressed (DE) miRNAs in BRV-infected cells in comparison with uninfected MA-104 cells were successfully investigated, 95 of which were upregulated and 65 were downregulated. The target messenger RNAs (mRNAs) of the DE miRNAs were examined by bioinformatics analysis. Functional annotation of the target genes with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) suggested that these genes mainly contributed to biological pathways, endocytosis, apoptotic process, trans-Golgi membrane, and lysosome. Pathways such as the mammalian target of rapamycin (mTOR) (mml-miR-486-3p and mml-miR-197-3p), nuclear factor kappa B (NF-κB) (mml-miR-204-3p and novel_366), Rap1 (mml-miR-127-3p), cAMP (mml-miR-106b-3p), mitogen-activated protein kinase (MAPK) (mml-miR-342-5p), T-cell receptor signaling (mml-miR-369-5p), RIG-I-like receptor signaling (mml-miR-504-5p), AMP-activated protein kinase (AMPK) (mml-miR-365-1-5p), and phosphatidylinositol-3-kinase/protein kinase B (PI3K/Akt) signaling (mml-miR-299-3p) were enriched. Moreover, real-time quantitative PCR (qPCR) verified the expression profiles of 23 selected DE miRNAs, which were consistent with the results of deep sequencing, and the 28 corresponding target mRNAs were mainly of regulatory pathways of the cellular machinery and immune importance, according to the bioinformatics analysis. Our study is the first to report a novel approach that uncovers the impact of BRV infection on the miRNA expressions of MA-104 cells, and it offers clues for identifying potential candidates for antiviral or vaccine strategies.

Bovine rotavirus (BRV) causes massive economic losses in the livestock industry worldwide. Elucidating the pathogenesis of BRV would help in the development of more effective measures to control BRV infection. The MA-104 cell line is sensitive to BRV and is thereby a convenient tool for determining BRV-host interactions. Thus far, the role of the microRNAs (miRNAs) of MA-104 cells during BRV infection is still ambiguous. We performed Illumina RNA sequencing analysis of the miRNA libraries of BRV-infected and mock-infected MA-104 cells at different time points: at 0 h postinfection (hpi) (just after 90 min of adsorption) and at 6, 12, 24, 36, and 48 hpi. The total clean reads obtained from BRV-infected and uninfected cells were 74,701,041 and 74,184,124, respectively. Based on these, 579 were categorized as known miRNAs and 144 as novel miRNAs. One hundred and sixty differentially expressed (DE) miRNAs in BRV-infected cells in comparison with uninfected MA-104 cells were successfully investigated, 95 of which were upregulated and 65 were downregulated. The target messenger RNAs (mRNAs) of the DE miRNAs were examined by bioinformatics analysis. Functional annotation of the target genes with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) suggested that these genes mainly contributed to biological pathways, endocytosis, apoptotic process, trans-Golgi membrane, and lysosome. Pathways such as the mammalian target of rapamycin (mTOR) (mml-miR-486-3p and mml-miR-197-3p), nuclear factor kappa B (NF-κB) (mml-miR-204-3p and novel_366), Rap1 (mml-miR-127-3p), cAMP (mml-miR-106b-3p), mitogen-activated protein kinase (MAPK) (mml-miR-342-5p), T-cell receptor signaling (mml-miR-369-5p), RIG-I-like receptor signaling (mml-miR-504-5p), AMP-activated protein kinase (AMPK) (mml-miR-365-1-5p), and phosphatidylinositol-3-kinase/protein kinase B (PI3K/Akt) signaling (mml-miR-299-3p) were enriched. Moreover, real-time quantitative PCR (qPCR) verified the expression profiles of 23 selected DE miRNAs,

INTRODUCTION
Rotaviruses (RVs) account for the severe diarrhea in children and young animals globally. RVs belong to the genus Rotavirus, which is one of the 15 genera of the Reoviridae family (Matthijnssens et al., 2012). They are non-enveloped, double-stranded, and segmented RNA viruses. The 11 segments encode 11 or 12 viral proteins, including six structural proteins (VPs) and five or six non-structural proteins (NSPs) (Estes and Greenberg, 2013). Based on VP6, genetically and antigenically, RVs are categorized into 10 groups (species A-J) (Bányai et al., 2017). Besides being a disease cause, species A rotaviruses (RVAs) have a large economic impact on human and animal species (Matthijnssens et al., 2010). Based on the antigenic features of VP4 and VP7, RVs are further classified into the G and P genotypes and serotypes, respectively (Estes and Greenberg, 2013). MicroRNAs (miRNAs) are a unique class of small RNAs (sRNAs), which are endogenous, single-stranded, non-coding RNAs (ncRNAs) of 19-24 nucleotides (nt) in length, expressed by all multicellular eukaryotes. This class has demonstrated significant evolutionary conservation (Friedman et al., 2009), which can inhibit the expressions of complementary target messenger RNAs (mRNAs) via transcriptional destabilization or translational inhibition (Bartel, 2009). Many studies have demonstrated that miRNAs are involved in a variety of biological processes, such as cell differentiation, development, homeostasis, proliferation, apoptosis, and signal transduction, in animals and plants (Axtell et al., 2011). RNA virus infection can mediate changes in the expression of the host miRNA machinery, leading to downstream changes in the host transcriptome that can be favorable to viral replication and pathogenesis. Notwithstanding, these changes in miRNA expression can also lead to increases in antiviral effector activities, resulting in decreased viral replication (Trobaugh and Klimstra, 2017). Mature miRNA guides RNA-induced silencing complexes (RISCs) to mRNAs bearing complementary target sites (Hammond et al., 2000). RISCs, which are minimally composed of a mature miRNA and an Argonaute protein, usually, but not invariably, bind to the 3untranslated region (UTR) of the targeted transcripts. Functional targets are generally fully complementary to nucleotides 2-7, preferably 2-8, at the 5 end of the miRNA, referred to as the miRNA seed (Bartel, 2009). However, there are a few examples of base pairing between the target transcript and the 3 half of the miRNA that can compensate for mismatches in the seed. Perhaps the best example of this phenomenon was provided by the miRNA lin-4 and the well-characterized target mRNA encoding lin-14, where robust inhibition of the expression of lin-14 was conferred by target sites lacking full seed homology (Ha et al., 1996). A single miRNA can target multiple mRNAs, and a single mRNA 3 -UTR can contain binding sites for several miRNAs (Bartel, 2009). The genes for miRNAs are distributed in different genomic locations, such as the intergenic regions, introns, and exons of protein-coding genes. Most miRNA genes are transcribed individually and scattered across the entire genome, whereas a portion of these genes is clustered and transcribed as one long primary miRNA (pri-miRNA) transcript and processed into individual precursor miRNAs (pre-miRNAs) (Altuvia et al., 2005). According to miRBase, humans encode >500 different miRNAs (Griffiths-Jones et al., 2008), which tend to be expressed in a developmental or tissue-specific manner (Landgraf et al., 2007). It has been demonstrated that individual miRNAs are capable of directly downregulating the expressions of hundreds of different mRNAs, and >30% of all mammalian mRNAs are thought to be regulated by miRNAs (Lewis et al., 2005). Some studies have characterized cellular miRNAs after infection with human immunodeficiency virus (HIV) (Chang et al., 2013), hepatitis C virus (HCV) (Norman and Sarnow, 2010), influenza virus (IV) (Skovgaard et al., 2013), West Nile virus (WNV) (Slonchak et al., 2014), Japanese encephalitis virus (JEV) (Zhang et al., 2015), dengue fever virus (DENV) (Liu et al., 2015), bluetongue virus (BTV) (Xing et al., 2016), and vesicular stomatitis virus (VSV) (Otsuka et al., 2007). Notably, one study reported on a comprehensive miRNA analysis of HT29 and MA104 cells infected with the human RV strain ZTR-68 (P[8] G1) (Zhou et al., 2016) and of Caco-2 cells infected with the human RV strain Wa (Tian et al., 2017). Here, our research aimed at analyzing the differentially expressed (DE) miRNAs of host cells before and after bovine rotavirus (BRV) infection. Data on the differential miRNAs induced by the interaction between BRV and host cells could fill the gap in viral pathogenesis, host innate antiviral defense, and viral replication and thus can be used for future development of novel antiviral drugs and vaccines against this virus.

Virus and Cell Culture
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)

Detection of Rotavirus dsRNAs by SDS-PAGE
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 R , 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 doublestranded RNA (dsRNA) were visualized with a silver stain and photographed (Gray and Desselberger, 2000).

RNA Purification and Quality Control
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 R , 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.

Small RNA Sequencing
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 basespecific 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 1 (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 2 (Enright et al., 2003) and RNAhybrid 3 (Krüger and Rehmsmeier, 2006).

Enrichment Analysis of the Differential MicroRNAs Target Genes
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 4 , 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) 5 . 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.

Real-Time Relative Quantitative PCR
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 R 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 R 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.

Statistical Analysis
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.

Confirmation of Bovine Rotavirus Replication in MA-104 Cells
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).      , microRNAs (miRNAs) by real-time quantitative PCR (qPCR) from 0 to 48 hpi. The differences in expression were determined using the 2 -CT method and with unpaired, two-tailed Student's t-test for grouped data. All experiments were performed in triplicate (n = 3). U6 was used as an internal reference. Data represent the mean ± SD (error bar). Significant difference thresholds: log2(fold change) > 2 and FDR-adjusted p < 0.05.

MicroRNA Expression and Differential Expression
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) Tables 5af). 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).

Interactions Between MicroRNAs and Messenger RNAs
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 FIGURE 8 | (A-F) Validation of target messenger RNAs (mRNAs) by real-time quantitative PCR (qPCR) from 0 to 48 hpi. The differences in expression were determined using the 2 -CT method and with unpaired, two-tailed Student's t-test for grouped data. All experiments were performed in triplicate (n = 3). β-actin was used as an internal reference. Data represent the mean ± SD (error bar). Significant difference thresholds: log2(fold change) > 2 and of FDR-adjusted p < 0.05. 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)

DISCUSSION
The expression of the host miRNA and the differential expression induced by viral infection were used to investigate the virushost interaction (Ghosh et al., 2009;Skalsky and Cullen, 2010). Some cellular miRNAs can mitigate viral replication either indirectly by regulating the host immune response or by directly targeting the virus (Slonchak et al., 2014). In contrast, other miRNAs can trigger viral replication by modulating the cellular environment (Jopling et al., 2005;Norman and Sarnow, 2010). Notwithstanding, there is no previous report on the relationship between BRV and cellular miRNAs. Therefore, in this study, the cellular DE miRNAs in response to BRV infection of MA-104 cells were identified using small miRNA sequencing of BRV-infected cells. We obtained 579 known miRNAs and 144 novel miRNAs and discovered 160 DE miRNAs in BRV-infected cells in comparison with the uninfected group, of which 95 were upregulated and 65 were downregulated. Bioinformatics analysis through miRanda and RNAhybrid revealed the target genes of the obtained DE miRNAs, and GO and KEGG analyses were performed to study the target mRNA functions and the biochemical metabolic pathways and signal transduction pathways (protein-protein interaction) of the target mRNAs involved in BP, CC, and MF. Moreover, real-time qPCR was used to validate the 23 DE miRNAs and 28 target mRNAs, which were harmonized with the deep sequencing and bioinformatics analyses. On the one hand, some of the validated miRNAs were upregulated [mml-let-7f-5p (0 hpi), mml-miR-99a-5p (0 and 6 hpi), mml-miR-486-3p (6 hpi), mml-miR-486-5p (6 hpi), mml-miR-127-3p (6 hpi), mml-miR-132-3p (12 hpi), mml-miR-369-5p (12 hpi -miR-197-3p (24 hpi), mml-miR-204-3p (24 hpi), mml-miR-504-5p (24 hpi), mml-miR-365-1-5p (36 hpi), and mml-miR-342-5p (12 and 24 hpi)]. In other studies, miR-21-3p was found to be downregulated in the MA-104 cell line in response to both human RV strain ZTR-68 (P[8] G1) (Zhou et al., 2016) and strain Wa (Tian et al., 2017); surprisingly, this miRNA was upregulated in our study, which could be an indication of the different types of pathogenesis between the human and BRV strains. This deserves further investigation. It was demonstrated that EBV induces miR-21 during latency, linking this miRNA with viral persistence (Cameron et al., 2008). Some studies have reported the collaboration between cellular miRNAs and the biological processes and cellular signal transduction pathways. It was elucidated that miRNAs are closely linked to biological pathways such as stress response, proliferation, differentiation, apoptosis, and autophagy (Frankel and Lund, 2012). Previous research works have demonstrated that, despite autophagy being able to act as an antiviral mechanism, some viruses use this machinery in favor of viral replication (Shintani and Klionsky, 2004). Also, there were some speculations over the suppression of apoptosis being good for RV replication. For example, the inhibition of premature apoptosis by the NSP1dependent PI3K3/Akt3/NF-κB pathway was suggested be favor viral growth of the simian RV strain SA11 (2-8 hpi) (Bagchi et al., 2010). Similarly, supportive proof has been reported for the PI3K/Akt/mTOR signaling pathway exerting an anti-autophagic role in SA11 rotavirus pathogenesis at 48 hpi, which in turn improves viral replication (Yin et al., 2018). Notwithstanding, another recent study has indicated that the decreased level of mTOR during early virus infection (4-6 hpi) could stimulate the autophagic pathway and help in the replication of RV (Mukhopadhyay et al., 2019). In addition, it was mentioned that the upregulation of let-7f-5p repressed several pro-apoptotic proteins (Tie et al., 2018). It was also confirmed that miR-99a-5p could target mTOR, while miR-486-3p and miR-486-5p were implicated in oncological and non-oncological conditions (ElKhouly et al., 2020). The downregulation of hsa-miR-127-3p could be the leading factor in the overexpression of type I interferon (IFN) . Besides, miR-132-3p and miR-369 suppressed the level of apoptosis, and miR-132-3p reduced the response of type I IFN to pave the way for influenza virus replication (Zhang F. et al., 2019;Qu et al., 2021;Wang J. et al., 2021). Moreover, miR-451-3p could mitigate apoptosis and autophagy (Tsai et al., 2018;Lv et al., 2021). Interestingly, in our study, mml-miR-127-3p, mml-miR-486-3p, and mml-miR-486-5p exhibited a surge in regulation at 6 hpi. mml-miR-99a-5p, mml-miR-132-3p, miR-369-5p, and mml-miR-451 likely play a role in cell survival and, thus, RV growth at 6, 12, and 24 hpi, which deserves further detailed study. Evidently, at 24 hpi, we observed that the DE miRNAs and their mRNAs were mainly of immune importance for defense against viral replication. On the contrary, the miR-197-3p and miR-504-5p in our study showed downregulation, with a previous report demonstrating that the inhibition of miR-197-3p represses cell apoptosis (Li et al., 2019). Additionally, the stimulation of miR-504 can enhance the apoptotic machinery (Hu L. et al., 2021). Together, it implies that the downregulation of miR-197-3p and miR-504-5p induced by BRV infection would suppress cell apoptosis to improve viral replication. Notably, the autophagic machinery was possibly operated via the collaboration between the downregulated miR-204 and the upregulated LC3-II protein (Xiao et al., 2011); therefore, the downregulation of miR-204-3p in our study could be implicated in the signaling transduction pathways for the virus strain replication. Because miR-21 has been reported as a pro-apoptotic miRNA that acts as a binding target of Bcl-2 , there was speculation around its participation as a pro-autophagic factor (Yang and Liang, 2015). As a result, the exact mechanism of the effect of both miR-204-3p and miR-21-3p on BRV replication needs further detailed investigation. At 36 hpi, our results showed the downregulation of mml-miR-365-1-5p and the upregulation of mml-miR-106b-3p. It has been mentioned that the upregulation of miR-365 promotes apoptosis via inactivation of the PI3K/AKT/mTOR pathway (Wang et al., 2020). Besides, miR-106b-3p has been used in the diagnosis of pancreatic cancer (Cao et al., 2016). Last but not least, mml-miR-299-3p and novel_366 exhibited a notable dramatic upregulation at 48 h, and it has been mentioned that miR-299-3p suppresses cell progression and induces apoptosis (Wang Z. et al., 2021). Evidently, mml-miR-342-5p was slightly downregulated at 12 hpi; however, it showed a marked downregulation at 24 hpi. Similarly, some novel DE miRNAs, such as novel_356 and novel_66, were upregulated at 24 hpi in comparison with 36 hpi, and novel_458 was also highly upregulated at 12 hpi compared to 36 hpi. Meanwhile, mml-miR99a-5p, mml-miR-411-3p, and mml-miR-411-5p exhibited the same upregulation levels at different time points (at 0 and 6 hpi, at 12 and 48 hpi, and at 36 and 48 hpi, respectively). Other reports have mentioned that miR-342-5p and miR-411-3p/5p were implicated in the formation of cancer (Lindholm et al., 2019;. In our study, real-time qPCR validated the bioinformatics results of the target mRNAs of the DE miRNAs that were implicated in the cellular biological processes and signaling pathways. Some of the upregulated targets are as follows: SERINC4 (0 hpi); PDPK1 (surge upregulation at 6 hpi); AZGP1 (dramatic upregulation at 12 hpi); ERAL1, ORMDL3, IFITM5, TRIM17, and mTORC1 (surge upregulation), OASL, MAP3K12, ISG20, and PHPT1 (Janus superfamily, surge upregulation), IKBKB, AKT1, and STAT2 (surge upregulation), and ISG15 (all at 24 hpi); and COG6 (at 36 hpi). MAPK15 was dramatically upregulated at 12 hpi in comparison with 24 hpi. Additionally, LAMB2 and BCL2A1 both exhibited a surge in upregulation level at 48 hpi. The other targets were downregulated, as follows: PARD6G (at 6 hpi); CD4 (at 12 hpi); TRBV16, TSC2, and BCL6 (at 24 hpi); and HSPB6, ACACA, and BDNF (dramatic decrease at 36 hpi). Interestingly, PDPK1 has been implicated in autophagosome biogenesis (Hu B. et al., 2021). It can modulate the autophagy pathway of our RV strain, which will need further detailed research, whereas PARD6B was shared in the cancer pathways (Cunliffe et al., 2012). COG6 and SERINC4 have been suggested to be inhibitors of HIV-1 replication (Liu et al., 2014;Qiu et al., 2020), while the expression of AZGP1, which has been reported to be associated with HPV status, was significantly high (Poropatich et al., 2019). ERAL1 suppressed RNA virus infection by facilitating RIG-I-like receptor signaling , while ORMDL3 was reported to be regulated by IRF-3 as a result of RSV infection by direct binding to the promoter of ORMDL3 . Evidently, STAT2, as a transcription regulatory factor, could increase the levels of the interferon-stimulated gene (ISG) and IFN transcripts in response to RV, establishing an antiviral state (Angel et al., 2012). The IFITM family (interferoninduced genes) was found to respond to type I and II IFNs (Yánez et al., 2020), and OASL, a family member of ISGs, was documented as one of the key contributors in the control of antiviral innate immunity, specially ISG20 and ISG15, which were reported to be induced or regulated by both type I and II IFNs (Zhu et al., 2015;Zheng et al., 2017). Moreover, TRIM17 can not only promote the degradation of the antiapoptotic protein MCL1 but also prevent the ubiquitination of other proteins and stabilize them by binding to other TRIM proteins and inhibiting their E3 ubiquitin ligase activity. This dual action confers several pivotal roles to TRIM17 in crucial cellular processes such as apoptosis and autophagy (Basu-Shrivastava et al., 2021). On the contrary, the MAPK family members have a close relationship with other viral signaling pathways, such as NF-κB, PI3K/Akt, and Janus kinase/signal transducer and activator of transcription (JAK-STAT) (Chander et al., 2021). Notwithstanding, the transcription factor and proto-oncogene B-cell CLL/lymphoma 6 (BCL6) could mediate the transcriptional repression of LITAF, which may inhibit autophagy in B cells during the germinal center reaction (Bertolo et al., 2013). HSPB6 was reported to suppress the MAPK family, including JNK and the PI3K/Akt pathway (Matsushima-Nishiwaki et al., 2011, and the downregulated ACACA curbed prostate cancer growth . Last but not least, the BCL2 protein was reported as an anti-apoptotic factor via its regulation of NF-κB target gene that exerts important pro-survival functions (Vogler, 2012).

CONCLUSION
In conclusion, our study elucidated a novel approach to reveal the miRNAs of the MA-104 cell line induced by BRV infection and the crosstalk between the DE miRNAs, target mRNAs, and the signaling pathways during BRV infection. At 6 hpi, BRV upregulated mml-miR-486-3p, which may, in turn, elevate the expression of target PDPK1, implicated in autophagosome formation (mTOR signaling pathway) in favor of viral replication. At 24 hpi, BRV downregulated mml-miR-204-3p, which may, in turn, upregulate its target AKT1. At 48 hpi, BRV upregulated both mml-miR-299-3p and novel_366, resulting in the upregulation of the targets LAMB2 and BCL2A1, respectively. These results (at 24 and 48 hpi) exerted anti-apoptotic and pro-survival functions via the PI3K/Akt and NF-κB signaling pathways, improving viral replication. On the other hand, at 24 hpi, the host immune system was activated against viral replication through the downregulation of both mml-miR-204-3p, which activated IKBKB (NF-κB signaling), and mml-miR-504-5p (RIG-I-like receptor signaling pathway), which triggered STAT2 (JAK-STAT signaling pathway), providing antiviral status (Figure 10). Thus, it could be of value to reduce the gap in BRV pathogenesis and to determine potential targets for the development of novel antiviral drugs and vaccines.

DATA AVAILABILITY STATEMENT
The analyzed data that supports the findings of this study are available as the Supplementary Material of this article. The data that support the findings of this study were deposited in GEO data repository under accession GSE196536, using the following contact links: (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi? acc=GSE196536) and (https://www.ncbi.nlm.nih.gov/geo/info/ linking.html). And our data are available from the corresponding author upon reasonable request.

AUTHOR CONTRIBUTIONS
GE designed and performed the practical experiments and drafted the manuscript. YC, CH, JC, and XC coordinated the project. AG contributed to the supervision and revision of the