Integrated Analysis of microRNAs and Metabolomics in Rat’s Serum Reveals Multi-action Modes of Qingfei Paidu Decoction for COVID-19 Treatment

Background During the ght against coronavirus disease 2019 (COVID-19) in China, Qingfei Paidu decoction (QFPDD) has been widely applied to treat COVID-19 patients. Retrospective studies showed that QFPDD could improve clinical outcomes of COVID-19. Thus, it is necessary and interesting to explore the action mode of QFPDD for further application and development.

water. The rats were orally administrated by warm QFPDD for 12 days at the dosage of 10 μL/g/d, which was equivalent to six times human dosage used in the treatment of COVID, and the control rats were parallelly treated with warm distilled water.

Rats' serum preparation and Constructed groups
On the rst day post the 12-days treatment for the rats, their blood was collected from the abdominal aorta in the absence of anticoagulants. The blood samples were kept static at 4°C for 2 h. Then the serum was separated from these blood samples by centrifugation at 4°C at 1600 × g for 10 minutes and stored at -80℃ individually until detection and analysis. In such a way, two groups were constructed: 12-day treated serum after oral administration of QFPDD (QFPDD group) and 12-day parallel serum (control group), which included 9 and 11 samples, respectively.
Sequencing of miRNAs pro le RNA isolation Total RNA was extracted from the blood serum using the Total RNA Extractor Trizol kit (B511311, Sangon, China) according to the manufacturer's protocol, and treated with RNase-free DNase to remove genomic DNA contamination.
RNA integrity was evaluated with a 1.0% agarose gel. Thereafter, the quality and quantity of RNA were assessed using a Nano Photometer ® spectrophotometer (IMPLEN, CA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). The high-quality RNA samples were subsequently submitted to the Sangon Biotech (Shanghai) Co., Ltd. for library preparation and sequencing. Concentration of total RNAs extracted from the serum of QFPDD and control group ranged from 0.534~0.720 μg/μL in 30 μL total volume, and integrity showed that the samples meet the requirements for sample preparation.
Small RNA library preparation and sequencing A total of 2μg total RNA per sample was used as input material for small RNA library preparations. Sequencing libraries were generated using NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB, USA.) following manufacturer's recommendations. Total RNA was used as the starting sample, the small RNA ends are directly connected with the adapter, followed by reverse transcription synthesis into cDNA. DNA fragments of 140-150 bp were separated by PAGE gel electrophoresis, and the cDNA library was recovered. Finally, library quality was assessed on the Agilent Bioanalyzer 2100 system, and libraries were sequenced on an Illumina HiSeq X-ten platform.

Raw Reads Quality Control
The rst step in our data analysis is removing adapters with software cutadapter (1.14), then low quality ( 20) bases from both 3' end and 5' end were trimmed with Trimmomatic (v0.36). Reads length ranging from 17bp and 35bp were reserved.
We used FastQC and in house scripts to generate QC report, which could help us to spot problems originating from either library preparation or sequencing process. We reviewed sequence quality metrics including per base quality scores (Q30), GC content, N content, Kmer content and duplicate levels, base bias at each position, copy number, length distribution, bases proportion, etc. Fastq to fasta conversion and removal of duplication were conducted for further data analysis.

Reads Filter
As other small RNA and fragments from mRNA degradation may get introduced in library preparation, reads mapped to other small RNA (rRNA, sRNA, snRNA, snoRNA) and mRNA sequence using blastn (2.6.0) and bowtie (1.1.1) respectively, were removed. The remaining reads were mapped to species genome sequence, and unmapped reads were ltered.

MiRNA analysis
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. The alignment result (mrd format) was used to detect miRNA modi cation, which originated from SNP and RNA-editing. As miRNA are evolutionary conserved, we summarized numerical counts and reads counts of known miRNA family.

Differentially Expressed MicroRNAs (DEMs) Detection
We used an R package edgeR (v3.18.1) to compare expression levels between sample pairs. We used fold change > 2 and p value < 0.05 as thresholds to de ne signi cantly differentially expressed miRNAs. To give an overview of the difference status of gene expression, we used volcano plots, MA plot and scatter plot to present the genes with differential expression between the paired samples.

Functional enrichment analysis
First, all genes were mapped to each term of the GO database (http://www. geneontology.org/) and the gene numbers of each GO term were calculated to get a gene list. Then the hypergeometric test was used to nd signi cantly enriched GO terms in differentially expressed genes (DEGs) compared with the genome background. After Bonferroni correction, GO terms with corrected P≤0.05 were de ned as signi cantly enriched GO terms. Similarly, by applying a cut-off criterion of Q-value< 0.05, pathway enrichment analysis was performed by comparing the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations with the whole gene background. To better understand the function of the signi cantly differentially expressed miRNA, we used R package clusterPro ler (3.4.4) to do functional enrichment of target genes. Bar plot and scatter plot were created to demonstrate top signi cantly enriched term (top 10 gene ontology and top 20 pathways). Top term gene-term and term-term interactive network plot were implemented in R package igraph (1.1.2).

Detection of metabolomics
Metabolites Extraction 100 μL of serum was transferred to an EP tube. After the addition of 400 μL of extract solution (acetonitrile: methanol = 1:1, containing isotopically labelled internal standard mixture), the samples were vortexed for 30 s, sonicated for 10 min in ice-water bath, and incubated for 1 h at -40 ℃ to precipitate proteins. Then the sample was centrifuged at 12000 rpm (RCF=13800(×g), R= 8.6cm) for 15 min at 4 ℃. The resulting supernatant was transferred to a fresh glass vial for analysis.
The quality control (QC) sample was prepared by mixing an equal aliquot of the supernatants from all the samples. The auto-sampler temperature was 4 ℃, and the injection volume was 3 μL.
The QE HFX mass spectrometer was used for its ability to acquire MS/MS spectra on information-dependent acquisition (IDA) mode in the control of the acquisition software (Xcalibur, Thermo). In this mode, the acquisition software continuously evaluates the full scan MS spectrum. The ESI source conditions were set as following: sheath gas ow rate as 30 Arb, Aux gas ow rate as 25 Arb, capillary temperature 350 ℃, full MS resolution as 60000, MS/MS resolution as 7500, collision energy as 10/30/60 in NCE mode, spray Voltage as 3.6 kV (positive) or -3.2 kV (negative), respectively.

Data preprocessing and annotation
The raw data were converted to the mzXML format using ProteoWizard and processed with an in-house program, which was developed using R and based on XCMS, for peak detection, extraction, alignment, and integration. SIMCA software (v16.0.2) was used for multivariate statistical analysis. Orthogonal partial least squares discriminant analysis (OPLS-DA) and student's t-test were employed to identify signi cantly changed metabolites (SCMs).

Results
Analysis of miRNAs pro les 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 pro les with untreated rats. Approximately 9.69 ~ 18.16 million reads were obtained from the 20 samples included in the two groups (Additional le 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 le 2). Intron, exon, tRNA, rRNA, miRNA, snRNA and snoRNA reads were annotated, respectively, as shown in Additional le 3. The length distribution of clean reads peaked at about 22 nt for all the samples. The identi cation 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 rst position. Approximately 133,462 ~ 22,046 unique sequences in the two groups' library were annotated as miRNA candidates (Additional le 3). A total of 585 known miRNA genes were identi ed (Additional le 4). The miRNAs ranged in size from 17 to 25 nt, with 71.62% from 21 to 22, and the base bias at the rst position of the identi ed 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 le 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 pro les and prediction of targets We normalized the expression of each miRNA (Additional le 6) as reads count per million (CPM) and calculated the QFPDD treated and untreated rat's serum miRNA expression ratio (Additional le 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.
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 le 9). The top 100 regulated genes ranked by score were speci cally 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 le 10). Especially, seven target genes might facilitate SARS-CoV-2 infection [100][101][102][103] or be related to symptoms of COVID-19 [104][105][106][107][108][109][110][111][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 B 0 AT1) [104][105][106][107][108], ADP-ribosylarginine hydrolase (Adprh) [103], C-C motif chemokine receptor 9 (Ccr9) [110,111], WNK lysine de cient protein kinase 1 (Wnk1) [112]. Fasn, VPS39 and Adprh may be required for SARS-CoV-2 replication [100,103] or infection [101,102]. ENaC-α, B 0 AT1 and CCR9 were strategically mimicked or regulated by SARS-CoV-2 infection, leading to respiratory symptoms or respiratory failure [109][110][111], intestinal in ammation and nutritional status [104][105][106][107][108] in COVID-19 patients. Moreover, Wnk1 could ght 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 rst ve 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 signi cant 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 le 11. Then, the gene numbers of each GO term, regulated by DEMs were calculated (Additional le 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 rst 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 le 13). The gene numbers of each KEGG ID were calculated (Additional le 14, Additional le 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 signi cantly 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 signi cantly expressed miRNAs and their target genes was established for the top KEGG enrichments (Additional le 15a, b). According to the network, the top ve 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 ve core genes (Table 3) may play a positive role in SARS-CoV-2 infection through regulating immunity, in ammatory cytokines, virus infection and pulmonary brosis. Table 2 Important KEGG pathways enriched in response to QFPDD treatment and their potential roles in SARS-CoV-2 infection induced by reported evidence.

Analysis of metabolomics and identi cation 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 le 16a) and negative ion mode (NEG) (Additional le 16b) respectively, and the raw data was processed according to the reported reference [226]. Brie y, the original data was de-noised based on relative standard deviation, and the missing values were lled up by half of the minimum value. Then, total ion current (TIC) was employed to normalize the data. After these processing, The nal dataset contained 4181 and 3446 peaks (metabolites) were left for subsequent analysis, in positive ion mode (POS (Additional le 17a) and negative ion mode (NEG) (Additional le 17b), respectively.
The nal dataset containing the information of peak number, sample name and normalized peak area was imported to  (Table 4). 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 nd that three SCMs, PC  (Fig. 4).

Discussion
In this study, 20 SD rats were randomly divided into two groups, namely QFPDD group (n = 9) and control group (n = 11). They were parallelly treated for 12 days with QFPDD and warm distilled water, respectively. The preparation and treatment of QFPDD in our study were in accordance with the clinical administration of QFPDD during COVID-19 patients' treatment [7]. At the endpoint, there is no difference between the two groups in the terms of weight, whole blood cell count, routine biochemistry indicators of the functions of liver and kidney (the data not showed here), suggesting the safety of QFPDD, to a certain extent. Notably, three rats in the QFPDD group (3/9, approximately 33.33%) were symptomized with mild diarrhea during the 12 days, of which one recovered itself and the other two got the symptom at the day of the endpoint. However, this phenomenon just obeyed the treatment principle of QFPDD from the view of Chinese medicine, while mild diarrhea is deemed to be detoxifying COVID-19 patients only if the times of diarrhea were no more than 3 per day and feces do not look like watery stool.
This study was conducted according to a study ow chart mainly divided into three steps: 1) obtaining DEMs from QFPDD treated serum of rats, 2) analyzing the multi-action modes of QFPDD though analysis of potential roles of DEMs or their regulated genes, 3) validating the regulation of DEMs by SCMs, especially the products of key enzyme, Plcg1, targeted by novel-89-mature. A total of 23 DEMs in response to QFPDD were identi ed in the rst step, with 13 known miRNAs and 10 novel miRNAs (Fig. 1). Then, the potential role of the 13 known miRNAs were summarized based on existing evidence (Table 1). Apart from miR-3594-3p and miR-218a-5p, the left 11 miRNAs have been reported to involve regulation of immunity, in ammation, multi-organ (including pulmonary) brosis, or virus infection (Table 1). Especially, miR-1-3p, miR-10b-3p and miR-199a-3p may possess antivirus effect on SARS-CoV-2 [21,22]. These potential roles of DEMs were further supported by the analysis of their predicted target genes as follows.
On the other hand, GO enrichment and KEEG enrichment analysis for the target genes allowed us to take an overview of functions of DEMs (Figs. 2-3). Important KEGG pathways and core genes obtained from KEGG enrichment network analysis were selected for further analysis (Tables 2 and 3). Based on literature review, both the important KEGG pathways and core genes regulated by DEMs involved the regulation of immunity, in ammation, multi-organ (including pulmonary) brosis, or virus infection (Tables 2 and 3). Impressively, each important KEGG pathway has already been reported to involve SARS-CoV-2 infection (Table 2). In addition, the core gene, SOS Ras/Rac guanine nucleotide exchange factor 1 (Sos1), was predicted to be targeted by SARS-CoV-2 encoded miRNA [164], and epidermal growth factor (Egf) involved regulation of angiotensin-converting enzyme 2 (ACE2), a receptor for SARS-CoV-2 infection [229] (Table 3).
Interestingly, one of the core genes, phospholipase C, gamma 1 (Plcg1) was a kind of lipid metabolic enzyme [228], which provided a possibility to validate the regulation of DEMs through the metabolites of enzyme of core genes ( Table 3).
Analysis of metabolomics in the serum of rats treated with QFPDD identi ed 5 SCMs (Table 4) Phospholipase C, gamma 1 (Plcg1) was targeted by novel-89-mature which was downregulated by QFPDD in SD rats.
Generally, the downregulation of novel-89-mature can decrease degradation of mRNA of phospholipase C, enhancing the activity of phospholipase C, which means that the less of novel-89-mature expressed, the higher the activity of phospholipase C [19]. Supportively, the enzymatic reaction from PCs to phosphorylcholine was strengthened by the enhancement of activity of phospholipase C, leading to the reduction of PCs and elevation of phosphorylcholine in the serum of rats. Thus, our ndings of the reduced levels of the three kinds of PCs and the elevated level of phosphorylcholine (Table 4)

Conclusion
Taking all the ndings together, we propose that QFPDD can modify immunity, in ammation, virus infection, pulmonary brosis through the regulation of miRNAs to improve treatment outcomes of COVID-19 (Fig. 5). The regulation of miRNAs was partly validated by change of substrates (PCs) and product (phosphorylcholine) of a core target gene, Phospholipase C, gamma 1 (Plcg1).

Declarations Acknowledgments
The ideas presented here were developed during completion of research projects funded by the Jiangxi Emergency Research Project TCM Application for COVID-19 Prevention and Treatment (NO. 2020 J001). This project was nancially supported by Jaingxi University of Chinese Medicine. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Authors' contributions
Conceptualization and formal analysis: Tie-Long Xu.