Analysis of miRNAs profiles in the serum of rats treated with QFPDD
High-throughput sequencing and annotation of small RNAs in SD rats
To evaluate the impact of QFPDD on miRNA expression in SD rats, we performed Illumina sequencing on small RNA libraries obtained from QFPDD treated rats and compared the miRNA profiles with untreated rats. Approximately 9.69 ~ 18.16 million reads were obtained from the 20 samples included in the two groups (Additional file 1). And after removal of adaptors, contaminants, and low-quality reads, approximately 8.20 ~ 15.52 million clean reads ranged from 17 ~ 35 nt were obtained (Additional file 2). Intron, exon, tRNA, rRNA, miRNA, snRNA and snoRNA reads were annotated, respectively, as shown in Additional file 3. The length distribution of clean reads peaked at about 22 nt for all the samples.
The identification of known and novel miRNAs
The known mature miRNA and precursor miRNA of the species were archived from mirBase (v21). We used mirDeep2 (v2.0.0.8) to discover novel miRNAs, predicate hairpin structure of precursor-miRNA, and quantify miRNA. To identify known miRNAs, we aligned the small RNAs to the known mature miRNAs and their precursors in miRBase 21 to obtain the miRNA counts as well as the base bias at the first position. Approximately 133,462 ~ 22,046 unique sequences in the two groups’ library were annotated as miRNA candidates (Additional file 3). A total of 585 known miRNA genes were identified (Additional file 4). The miRNAs ranged in size from 17 to 25 nt, with 71.62% from 21 to 22, and the base bias at the first position of the identified miRNAs showed a strong preference for ‘U’ at the 5’-end, which is coordinated with previous studies [20].
After removing the snRNAs, snoRNAs, rRNAs, tRNAs, and known miRNAs, we aligned the remaining unannotated reads against the Rattus norvegicus genome to predict novel miRNAs. A total of 1622 potential novel miRNA reads were mapped in control and QFPDD treated samples. The miRNAs are shown in Additional file 5, along with their Dicer cleavage site, minimum free energy, frequency of reads and typical secondary structures of the characteristic stem-loop hairpins for the predicted precursors.
Changes in miRNA expression profiles and prediction of targets
We normalized the expression of each miRNA (Additional file 6) as reads count per million (CPM) and calculated the QFPDD treated and untreated rat’s serum miRNA expression ratio (Additional file 6). We further analyzed the expression data by calculating the fold changes and P-values based on the expression ratios and plotted these data as a scatter plot. Changes after oral administration of QFPDD were apparent (Fig. 1a). The upregulated miRNAs were more numerous than the downregulated miRNAs.
Using a |fold-change| ≥ 2 and P value ≤ 0.05 to screen the miRNAs, a total of 23 DEMs were identified (Additional file 7 and 8). Of these miRNAs, 20 were statistically upregulated (miR-10b-3p, miR-1-3p, miR-100-5p, miR-199a-3p, miR-1b, miR-205, miR-206-3p, miR-214-3p, miR-218a-5p, miR-34b-5p, miR-3594-3p, miR-383-3p, miR-466b-5p, novel-1149-mature, novel-227-mature, novel-280-mature, novel-595-mature, novel-667-mature, novel-754-mature and novel-98-star), and 3 were significantly downregulated (novel-1138-mature, novel-614-mature, and novel-89-mature) in QFPDD treated group. Ten of the significantly differentially expressed miRNAs are putative novel miRNA (Fig. 1b).
Multi-action modes of QFPDD directly induced by the known DEMs
Of the 13 known DEMs, 11 had been reported to be involved by the regulation of immunity, inflammation, multi-organ fibrosis, and virus infection (Table 1). In particular, miR-1-3p, miR-10b-3p and miR-199a-3p may possess antivirus effect on SARS-CoV-2 [21, 22], and miR-199a-3p, miR-206-3p and miR-214-3p can attenuate the pulmonary fibrosis [23–25]. However, it is contradictory that miR-34b-5p is significantly upregulated by QFPDD, which may promote inflammation and pulmonary fibrosis [26, 27].
Table 1
The potential role of DEMs during SARS-Cov-2 infection induced by reported evidence.
miRNA name
|
Potential role during SARS-CoV-2 infection induced by reported evidence
|
miR-1-3p
|
Involving regulation of immunity [28], inflammation [28–30], multi-organ fibrosis [31–33], and virus infection, e.g., SARS-CoV-2, EBOV, IAV [21, 34–36]
|
miR-1b
|
Involving regulation of inflammation [37]
|
miR-10b-3p
|
Involving regulation of immunity [38–40], HBV-associated liver fibrosis [41] and virus infection, e.g., SPPV, SARS-CoV-2, EBV, KSHV, HIV and many other retroviruses [22, 42–44]
|
miR-34b-5p
|
Involving regulation of immunity [45, 46], pulmonary inflammation [26, 47–49] and fibrosis [27] and virus infection, e.g., ALV-J, RSV [50, 51]
|
miR-100-5p
|
Involving regulation of immunity [52, 53], inflammation [54–56], virus infection, e.g., H5N1, HBV [57, 58]
|
miR-199a-3p
|
Involving regulation of immunity [59–61], inflammation [23, 62–65], multi-organ fibrosis [23, 66–70], virus infection [71], e.g., HBV, PRRSV, SARS-CoV-2, HCV, HBV [21, 67, 72–74]
|
miR-205
|
Involving regulation of inflammation [75–79], multi-organ (especially pulmonary) fibrosis [80–82], virus infection, e.g., HCV, HIV, HBV, EBV, HPV [83–87]
|
miR-206-3p
|
Involving regulation of immune [88], multi-organ (especially pulmonary) fibrosis [24], virus infection, e.g., BoDV [89]
|
miR-214-3p
|
Involving regulation of immune [90], inflammation [91–96],multi-organ (especially pulmonary) fibrosis [25, 93], virus infection, e.g., HCMV [97]
|
miR-383-3p
|
Involving regulation of inflammation [98]
|
miR-466b-5p
|
Involving regulation of inflammation [99]
|
Note: SARS-CoV-2, severe acute respiratory syndrome coronavirus-2; EBOV, ebola virus; IAV, influenza A virus; SPPV, sheep pox virus; KSHV, kaposi sarcoma-associated herpesvirus; HIV, human immunodeficiency virus; ALV-J, Avian leukosis virus subgroup J; RSV, respiratory syncytial virus; H5N1, Influenza A virus H5N1; HBV, hepatitis B virus; PRRSV, porcine reproductive and respiratory syndrome virus; HCV, hepatitis C virus; EBV, epstein-barr virus; HPV, human papillomavirus; BoDV, Borna disease virus; HCMV, Human cytomegalovirus. |
Multi-action modes of QFPDD indirectly induced by the predicted targets of DEMs
Action mode induced by the top 100 predicted targets ranked by score
To investigate the regulated miRNAs functions, we predicted the potential targets of the 23 DEMs. A total of 2088 targets, covering 1636 regulated genes were predicted (Additional file 9). The top 100 regulated genes ranked by score were specifically selected out for further analysis, of which 42% (42/100) were reported to be involved in regulation of at least one kind of virus infection (Additional file 10). Especially, seven target genes might facilitate SARS-CoV-2 infection [100–103] or be related to symptoms of COVID-19 [104–112], including fatty acid synthase (Fasn) [100], VPS39 HOPS complex subunit (Vps39) [101, 102], sodium channel epithelial 1 alpha subunit (Scnn1α or ENaC-α) [109, 113], solute carrier family 6 member 19 (Slc6a19 or B0AT1) [104–108], ADP-ribosylarginine hydrolase (Adprh) [103], C-C motif chemokine receptor 9 (Ccr9) [110, 111], WNK lysine deficient protein kinase 1 (Wnk1) [112]. Fasn, VPS39 and Adprh may be required for SARS-CoV-2 replication[100, 103] or infection [101, 102]. ENaC-α, B0AT1 and CCR9 were strategically mimicked or regulated by SARS-CoV-2 infection, leading to respiratory symptoms or respiratory failure [109–111], intestinal inflammation and nutritional status [104–108] in COVID-19 patients. Moreover, Wnk1 could fight against cytokine storm caused by sepsis in COVID-19 patients [112].
Overview of functions induced by GO enrichment analysis
GO enrichment covers three domains: biological processes partaken by target genes, cellular components containing target genes, and molecular functions possessed by target genes. GO categorization showed that the first five biological process of target genes were the cellular process, single-organism process, metabolic process, response to stimulus, multicellular organismal process and cellular component organization or biogenesis. Five molecular function of target genes, namely binding, catalytic activity, transporter activity, molecular function regulator, and signal transducer activity were most enriched. The genes located at cell, organelle, membrane, macromolecular complex, or membrane-enclosed lumen were most significant among the cellular component GO terms (Fig. 2a). More potential roles covered by target genes in every domain were presented in Fig. 2a and Additional file 11. Then, the gene numbers of each GO term, regulated by DEMs were calculated (Additional file 12), and based on a cut-off criterion of P value < 0.05, top ten GO terms with most gene numbers in every domain were showed in Fig. 2b.
Multi-action modes induced by KEGG pathway analysis
KEGG pathway is a collection of manually drawn pathway maps representing our knowledge on the aspects of molecular interaction, reaction and relation networks for metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, human diseases, and drug development. KEGG categorization showed that the first ten pathways involved signal transduction, infectious diseases, endocrine system, cancers, immune system, signaling molecules and interaction, cellular community, nervous system, transport and catabolism, and development (Fig. 3a, Additional file 13). The gene numbers of each KEGG ID were calculated (Additional file 14, Additional file 15a), based on a cut-off criterion of Q-value < 0.05, top twenty KEGG enrichments were selected for individual analysis (Fig. 3b). The KEGG pathway enrichment analysis revealed some important pathways that were significantly enriched in response to QFPDD treatment, including the NF-kappa B signaling pathway, Cytokine-cytokine receptor interaction, MAPK signaling pathway, signaling pathways regulating pluripotency of stem cells, Jak-STAT signaling pathway, Toll-like receptor signaling pathway, and RAS signaling pathway (Table 2). In addition, a putative network between significantly expressed miRNAs and their target genes was established for the top KEGG enrichments (Additional file 15a, b). According to the network, the top five core genes ranked by edges were selected for function analysis (Table 3). According to the reported evidence showed in Tables 2–3, these important KEGG pathways (Table 2) and the five core genes (Table 3) may play a positive role in SARS-CoV-2 infection through regulating immunity, inflammatory cytokines, virus infection and pulmonary fibrosis.
Table 2
Important KEGG pathways enriched in response to QFPDD treatment and their potential roles in SARS-CoV-2 infection induced by reported evidence.
KEGG pathway
|
Potential role during SARS-CoV-2 infection induced by reported evidence
|
NF-kappa B signaling pathway
|
Involving regulation of immunity [114, 115], inflammatory cytokines [114, 116–120], pulmonary fibrosis [121], and virus infection, e.g., SARS-CoV-2, H1N1 [115–117, 119–124]
|
Cytokine-cytokine receptor interaction
|
Involving regulation of immunity and inflammatory cytokines [125, 126], virus infection, e.g., SARS-CoV-2, DHAV, PDCoV, RSV [125, 127–131]
|
MAPK signaling pathway
|
Involving regulation of immunity and inflammatory cytokines [132–137], virus infection, e.g., PDCoV, HIV, H1N1, SARS-CoV-2, FMDV, IAV [133–144]
|
signaling pathways regulating pluripotency of stem cells
|
Involving regulation of virus infection, e.g., SARS-CoV-2, IAV [145–153]
|
Jak-STAT signaling pathway
|
Involving regulation of immunity [154], inflammatory cytokines [154–160], virus infection [161], e.g., SARS-CoV-2, HBV, IV, Flavivirus [154–157, 159, 160, 162–167]
|
Toll-like receptor signaling pathway
|
Involving regulation of inflammatory cytokines [157, 168], virus infection [71], e.g., HCV, EBV, SARS-CoV-2, IV [157, 168–171]
|
RAS signaling pathway
|
Involving regulation of virus infection, e.g., herpesviruses, HCMV, HTLV, HSV, HIV, reovirus, SARS-CoV-2, IV [172–180]
|
Note: SARS-CoV-2, severe acute respiratory syndrome coronavirus-2; DHAV, duck hepatitis A virus; PDCoV, porcine deltacoronavirus; FMDV, Foot-and-Mouth Disease Virus; RSV, respiratory syncytial virus; HIV, human immunodeficiency virus; IAV, influenza A virus; HBV, hepatitis B virus; IV, influenza virus; HCV, hepatitis C virus; EBV, epstein-barr virus; HCMV, Human Cytomegalovirus; HTLV, human T-cell leukemia virus; HSV, herpes simplex virus; H1N1, Influenza A virus H1N1 2009. |
Table 3
Core genes obtained from KEGG enrichment network analysis.
Gene name
|
Regulated
miRNA
|
Edges
|
Potential role during SARS-CoV-2 infection induced by reported evidence
|
SOS Ras/Rac guanine nucleotide exchange factor 1 (Sos1)
|
novel-754-mature
|
7
|
Predicted to be targeted by SARS-CoV-2 encoded miRNA [164], involving regulation of inflammatory cytokines [181], virus infection, e.g., SARS-CoV-2, enterovirus 71, HCV [164, 182, 183], development of fibrosis [184, 185]
|
epidermal growth factor (Egf)
|
novel-614-mature
|
7
|
Involving regulation of ACE2 [186], immunity [187–189], inflammatory cytokines [190–192], virus infection [193–203], especially SARS-CoV-2 [196, 198–201, 204–206], multi-organ (especially pulmonary) fibrosis [207–210]
|
phospholipase C, gamma 1 (Plcg1)
|
novel-89- mature
|
6
|
Involving regulation of immunity [211, 212], virus infection, e.g., HTLV, DHBV [213–215], liver fibrosis [216]
|
SHC adaptor protein 3 (Shc3)
|
novel-614-mature
|
5
|
Involving regulation of immunity [217, 218]
|
A-Raf proto-oncogene, serine/threonine kinase (Araf)
|
miR-3594-3p/novel-280-mature
|
5
|
Involving regulation of immunity [219, 220], inflammatory cytokines [221, 222], virus infection, e.g., H1N1, HBV [223–225]
|
Note: SARS-CoV-2, severe acute respiratory syndrome coronavirus-2; HTLV, human T-cell leukemia virus; HCV, hepatitis C virus; DHBV, duck hepatitis B virus; HBV, hepatitis B virus; H1N1, Influenza A virus H1N1 2009. |
Analysis of metabolomics and identification of SCMs in the serum of rats treated with QFPDD
In this study, 5339 and 4434 peaks (metabolites) were detected in positive ion mode (POS)(Additional file 16a) and negative ion mode (NEG) (Additional file 16b) respectively, and the raw data was processed according to the reported reference [226]. Briefly, the original data was de-noised based on relative standard deviation, and the missing values were filled up by half of the minimum value. Then, total ion current (TIC) was employed to normalize the data. After these processing, The final dataset contained 4181 and 3446 peaks (metabolites) were left for subsequent analysis, in positive ion mode (POS༉(Additional file 17a) and negative ion mode (NEG) (Additional file 17b), respectively.
The final dataset containing the information of peak number, sample name and normalized peak area was imported to SIMCA16.0.2 software package (Sartorius Stedim Data Analytics AB, Umea, Sweden) for multivariate analysis. Data was scaled and logarithmic transformed to minimize the impact of both noise and high variance of the variables. Then, supervised orthogonal projections to latent structures-discriminate analysis (OPLS-DA) was applied to identify SCMs between QFPDD and control groups [227] (Additional files 18a and 18b). In the OPLA-DA analysis, the value of variable importance in the projection (VIP) of the first principal component was obtained. It summarizes the contribution of each variable to the model. The metabolites with VIP > 1, p < 0.05 (calculated by student t-test) were considered as SCMs. Finally, the 10 SCMs were identified, including phosphorylcholine, 13S-hydroxyoctadecadienoic acid, PC(16:1 (9Z)/16:1(9Z)), PC(P18:1(11Z)/22:5(4Z,7Z,10Z,13Z,16Z)), PC(22:5(4Z,7Z,10Z,13Z,16Z)/P-18:0), beta-thujaplicin, octadecylamine, carnosol, glycyrrhetinic acid, 3-(3-Hydroxyphenyl) propanoic acid, of which, six SCMs belongs to lipids or lipid-like molecules (Table 4).
Table 4
The significantly changed metabolites (SCMs) in the serum of rats treated with QFPDD.
Compound name
|
SuperClass
|
VIP
|
Scan mode
|
Trend
(T/C)
|
P
value
|
Fold change
|
Log2
foldchange
|
Phosphorylcholine
|
Organonitrogen compounds
|
2.14
|
ESI+
|
↑
|
0.04
|
1.30
|
0.38
|
13S-hydroxyoctadecadienoic acid
|
Lipids or lipid-like molecules
|
2.00
|
ESI+
|
↑
|
0.04
|
1.35
|
0.43
|
PC(16:1(9Z)/16:1(9Z))
|
Lipids or lipid-like molecules
|
1.99
|
ESI+
|
↓
|
0.03
|
0.81
|
-0.30
|
PC(P-18:1(11Z)/22:5(4Z,7Z,10Z,13Z,16Z))
|
Lipids or lipid-like molecules
|
2.35
|
ESI+
|
↓
|
0.04
|
0.81
|
-0.30
|
PC(22:5(4Z,7Z,10Z,13Z,16Z)/P-18:0)
|
Lipids or lipid-like molecules
|
2.24
|
ESI+
|
↓
|
0.02
|
0.64
|
-0.64
|
beta-Thujaplicin
|
Hydrocarbon derivatives
|
2.16
|
ESI+
|
↓
|
0.02
|
0.83
|
-0.27
|
Octadecylamine
|
Organonitrogen compounds
|
1.95
|
ESI+
|
↑
|
0.04
|
1.23
|
0.30
|
Carnosol
|
Lipids or lipid-like molecules
|
1.99
|
ESI+
|
↑
|
0.04
|
1.24
|
0.30
|
Glycyrrhetinic acid
|
Lipids or lipid-like molecules
|
2.96
|
ESI−
|
↑
|
0.004
|
3.05
|
1.61
|
3-(3-Hydroxyphenyl) propanoic acid
|
Phenylpropanoids or polyketides
|
2.65
|
ESI−
|
↑
|
0.04
|
1.41
|
0.49
|
Note: ESI+ and ESI− refer to positive ion mode (POS) and negative ion mode (NEG), respectively; Fold change calculated from the arithmetic mean values of each group; Trend ↑, higher than control levels; ↓, lower than control levels. |
Statistical and theoretical linear correlation between DEMs and SCMs
To verify the regulation of DEMs, we performed an analysis of correlation between the 23 DEMs and 10 SCMs. Especially, the metabolites related to the core enzyme gene, phospholipase C, targeted by novel-89-mature, were paid more attention.
It was interesting to find that three SCMs, PC (16:1(9Z)/16:1(9Z)), PC (P18:1(11Z)/22:5(4Z,7Z,10Z,13Z,16Z)), PC (22:5(4Z,7Z,10Z,13Z,16Z)/P-18:0) belonged to a kind of phosphatidylcholines (PCs) which could be metabolized into phosphorylcholine by phospholipase C [228]. It meant that phosphorylcholine (another SCM) was the product of PCs metabolized by phospholipase C. More importantly, phospholipase C, gamma 1 (Plcg1) was predicted to be targeted by novel-89-mature.
To further confirm their correlation, we conducted an analysis of linear correlation between the 10 SCMs and novel-89-mature, along with the other 22 DEMs. Interestingly, novel-89-mature was significantly positively correlated with PC (P18:1(11Z)/22:5(4Z,7Z,10Z,13Z,16Z)), PC (22:5(4Z,7Z,10Z,13Z,16Z)/P-18:0), with the Pearson product moment coefficient (r) being 0.663 (P = 0.001) and 0.536 (p = 0.015) respectively. Meanwhile, no significant correlation was found between novel-89-mature and the other three SCMs (Fig. 4).
Confusingly, miR-10b-3p and novel-89-mature were found to be significantly positively correlated with 13S-hydroxyoctadecadienoic acid and beta-thujaplicin, with r = 0.761 (P < 0.0001) and r = 0.575 (P < 0.01) respectively. Both miR-34b-5p and miR-466b-5p were statistically positively correlated with octadecylamine, with r = 0.579 (P < 0.01) and r = 0.467 (P < 0.05), respectively. Moreover, novel-89-mature was negatively correlation with 3-(3-Hydroxyphenyl) propanoic acid, with r=-0.56 (P < 0.05). However, we only found statistical correlation between these miRNAs and metabolites, waiting for further theoretical interpretations. Except the linear correlation mentioned above, no other linear correlation was found between the 23 DEMs and 10 SCMs.