key: cord-1000095-f7dhquuo
authors: Tian, Saisai; Zheng, Ningning; Zu, Xianpeng; Wu, Gaosong; Zhong, Jing; Zhang, Jinbo; Sheng, Lili; Liu, Wei; Wang, Chaoran; Ge, Guangbo; Han, Jingyan; Zhao, Jing; Li, Houkai; Zhang, Weidong
title: Integrated hepatic single-cell RNA sequencing and untargeted metabolomics reveals the immune and metabolic modulation of Qing-Fei-Pai-Du decoction in mice with coronavirus-induced pneumonia
date: 2022-01-02
journal: Phytomedicine
DOI: 10.1016/j.phymed.2021.153922
sha: 84fb00218f7c118cf30afa161dcb90389b0805d5
doc_id: 1000095
cord_uid: f7dhquuo

BACKGROUND: Although Qing-Fei-Pai-Du decoction (QFPDD) is extensively used clinically to treat COVID-19 patients, the mechanism by which it modulates the immunological and metabolic functions of liver tissue remains unknown. PURPOSE: The purpose of this study is to investigate the mechanism of action of QFPDD in the treatment of mice with coronavirus-induced pneumonia by combining integrated hepatic single-cell RNA sequencing and untargeted metabolomics. METHODS: We developed a human coronavirus pneumonia model in BALB/c mice by infecting them with human coronavirus HCoV-229E with stimulating them with cold-damp environment. We initially assessed the status of inflammation and immunity in model mice treated with or without QFPDD by detecting peripheral blood lymphocytes and inflammatory cytokines. Then, single-cell RNA sequencing and untargeted metabolomics were performed on mouse liver tissue. RESULTS: HCoV-229E infection in combination with exposure to a cold-damp environment significantly decreased the percentage of peripheral blood lymphocytes (CD4(+) and CD8(+) T cells, B cells) in mice, which was enhanced by QFPDD therapy. Meanwhile, the levels of inflammatory cytokines such as IL-6, TNF-α, and IFN-γ were significantly increased in mouse models but significantly decreased by QFPDD treatment. Single-cell RNA sequencing analysis showed that QFPDD could attenuate disease-associated alterations in gene expression, core transcriptional regulatory networks, and cell-type composition. Computational predictions indicated that QFPDD rectified the observed aberrant patterns of cell-cell communication. Additionally, the metabolic profiles of liver tissue in the Model group were distinct from mice in the Control group, and QFPDD significantly regulated hepatic purine metabolism. CONCLUSION: To the best of our knowledge, this study is the first to integrate hepatic single-cell RNA sequencing and untargeted metabolomics into a TCM formula and these valuable findings indicate that QFPDD can improve immune function and reduce liver injury and inflammation.

Coronavirus disease 2019 caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a newly severe global pandemic that has already claimed millions of lives, underscoring the critical need for better therapeutic strategies. SARS-CoV-2 shares ~79% sequence similarity with SARS-CoV and ~50% similarity with the Middle East respiratory syndrome (MERS) coronavirus, all of which cause severe respiratory symptoms Jothimani et al., 2020; Lu et al., 2020) . Interestingly, in addition to the respiratory system, a significant majority of COVID-19 patients had varying degrees of liver damage (Chai et al., 2020; Jothimani et al., 2020) . Numerous studies have reported that patients with severe symptoms are more likely to exhibit abnormal liver function than those with mild symptoms (Wang et al., 2020b) . Additionally, earlier studies demonstrated that ACE2 receptors are expressed at a significantly higher level in cholangiocytes, indicating that the liver might be a target organ for SARS-CoV-2 (Seow et al., 2021) .

Dysregulation of metabolism, including the alteration in glycerophospholipids, sphingolipids, steroid hormones, amino acids, and fatty acids metabolism in COVID-19 patients, has also been reported in COVID-19 patients (Shen et al., 2020; Song et al., 2020; Thomas et al., 2020) , which is associated with inflammation and immune regulation (Cain and Cidlowski, 2017; Sevastou et al., 2013; Weigert et al., 2009 ). Therefore, disordered metabolism is a typical characteristic of patients.

Clinical studies on prospective treatments are currently ongoing, and no effective treatment has been confirmed Shin et al., 2020) . In China, traditional Chinese medicine (TCM) has played a significant role in the treatment of COVID-19, with the seventh edition of the COVID-19 guideline recommending Qing-Fei-Pai-Du decoction (QFPDD), for therapeutic use in light, moderate, and severe patients.

Clinical data from 701 confirmed that QFPDD has a successful cure rate of more than 90% against COVID-19 . Additionally, when QFPDD was combined with western medicine, the extent of multi-organ damage was reduced and the 7 immune regulation and anti-inflammation effects were improved compared to when western medicine was used alone (Xin et al., 2020) , implying that QFPDD played its therapeutic role via systemic regulation on the whole body. Therefore, understanding the mode of action of QFPDD's treatment of COVID-19 might be extremely beneficial in the war against SARS-CoV-2. Previous network pharmacological studies have demonstrated that QFPDD may target liver tissue and modulate immunological and inflammatory responses to function as a detoxifier Zhao et al., 2020) . However, these previous investigations were mostly predictive and lacked comprehensive analysis of the detoxifying organ liver and experimental validation, thus reducing the trust in the conclusions. To address these issues, we performed single-cell sequencing and metabolomic analysis on liver tissue from model mice treated with QFPDD.

This work aims to systematically investigate the mechanism of action for QFPDD on the mice pneumonia model. We firstly evaluated the effect of QFPDD on the mice pneumonia model by combining human coronavirus HCoV-229E with cold-damp environment exposure, which has been shown to replicate clinical pneumonia characteristics (Bao et al., 2020a; Bao et al., 2020c) . Based on the single-cell RNA sequencing (scRNA-seq) data and untargeted metabolomics data, we further performed extensive data mining, including variation analysis of cell type composition, differential expression analysis, cell-type-differential expression analysis, pathway analysis, transcriptional regulatory networks analysis, cell-cell communication analysis, cytokine and inflammatory score analysis and metabolic pathway analysis to explore the mechanism of action of QFPDD in multiple dimensions. This study also was the first report with direct experimental evidence supporting the clinical benefits of QFPDD, to the best of our knowledge, for improving the immune function of COVID-19 patients in the context of the pneumonia model of combined cold-damp environment exposure with HcoV-229E virus infection on mice. 8

Mice TNF-α， IFN-γ， IL-6 ELISA kits were purchased from Bio-Techne Company, USA. Human coronavirus 229E (HCoV-229E) was provided by the Institute of Medicinal Biotechnology Chinese Academy of Medical Sciences. PerCP-Cyanine5.5

Anti-Mouse CD4 (RM4-5), APC Anti-Mouse CD8a (53-6.7), PE Anti-Mouse CD19 (1D3), and RBC Lysis Buffer (10X) were purchased from TONBO Biosciences (USA). Methoxyamine HCl, fatty acid methyl ester (C7-C30, FAMEs) standards, pyridine, anhydrous sodium sulfate were obtained from Sigma-Aldrich (St. Louis,

trimethylchlorosilane (MSTFA, with 1% TMCS), hexane, dichloromethane, chloroform, and acetone were purchased from Thermo-Fisher Scientific (FairLawn, NJ, USA).

The chemical constituents of QFPDD were identified using a UHPLC-Q-Orbitrap High-Resolution Mass Spectrometry method. A total of 405 compounds in the aqueous extract of QFPDD were identified, including 40 different alkaloids, 162 flavonoids, 44 organic acids, 71 triterpene saponins, and 88 other different compounds. More detailed methods and results were shown in the published article (Liu et al., 2021) .

Animal experiments were supervised by Professor Xiaolan Cui and were conducted in an ABSL-2 laboratory at China Academy of Chinese Medical Sciences.

All procedures were carried out in accordance with the manual for the Care and Use Briefly, SPF grade BALB/c mice (10~12g) were purchased from the Beijing Huafukang Bioscience company. After one week of acclimatization, mice were randomly assigned into 5 groups (n = 8): Control (Con) group, Model group, low dosage QFPDD treatment (QFPDD_L) group, middle dosage QFPDD treatment (QFPDD_M) group, and high dosage QFPDD treatment (QFPDD_H) group. QFPDD (Dry powder) was provided by the JOINTOWN Pharmaceutical Group (Wuhan, China). Except for mice in the Con group, mice were kept in an artificial climate box with 90% ± 3% of relative humidity, (4 ± 2) ℃ for 4 h each day for 7 days. On the 5th and 6th days, mice were anesthetized with ether and infected via nose-dripping and were anesthetized with ether and infected with 50 μL of HcoV-229E (100TCID50) via nose-dripping method, whereas mice in the Con group were concurrently infected with culture solution in the absence of HcoV-229E. QFPDD was administered orally to mice in the QFPDD_L, QFPDD_M, and QFPDD_H groups at a dose of 0.88 g, 1.76 g, 3.52 g dry powder/kg/d once a day, respectively, which was equivalent to 0.5, 1, 2 times the clinical dosage when infected with HcoV-229E. The Con and Model groups received a similar volume of saline. On the 4 th day after infection, all mice were weighed and blood was obtained via the orbital vein into saline to determine the lymphocyte type ratio. Lung and liver tissues were collected for inflammatory cytokine analysis, scRNA-seq and untargeted metabolomics, respectively.

Three drops of blood were added into 10 mL of PBS, and the sample was centrifuged at 1600 r·min -1 for 5 min at room temperature. The supernatant was obtained, and 1 mL of red cell lysate was added to resuspend the red precipitate. After 5-10 min of lysis at room temperature, 10 mL of PBS was added to stop the lysis. The sample was centrifuged for 5 min at 2000 r·min -1 at 4 ℃. The supernatant was discarded and the precipitate was re-suspended in 10 mL of PBS. The sample was centrifuged for 5 min at 2000 r·min -1 and 4 ℃. The supernatant was discarded and 200 μL of blocking solution (PBS containing 5% of FBS) was added to suspend the residue. The sample was transferred into a 1.5 ml tube, blocked for 30 min at 4 ℃, and centrifuged for 5 min at 2000 r·min -1 and 4 ℃. The supernatant was discarded, and 50 μL of antibody dilutant containing 0.3 μL of anti-mouse CD19-PE antibody, 0.3 μL of anti-mouse CD4-PerCP-Cy5.5 antibody, and 0.3 μL of anti-mouse CD8a-APC antibody was added to the residue. After 30 min, the sample was resuspended in 1 mL of PBS and centrifuged at 2000 r·min -1 for 5 min at 4℃. The supernatant was discarded, and the cells were re-suspended in 200 μL PBS (containing 2% FBS) and transferred to tubes for flow cytometry analysis.

50 mg of lung tissue was added to 500 μL of saline and homogenized using an ultrasonic cell fragmentation device. The samples were centrifuged at 1000g·min -1 for 10 min at 4℃. The supernatant was collected and TNF-α, IL-6, IFN-γ levels were determined using the ELISA technique according to the manufacturers' instructions.

Liver tissue (each n=1) was randomly selected from each of the following 3 groups: model group, QFPDD group (At high dosage QFPDD treatment), and Con group mice. The tissues were dissociated into a single-cell suspension and the 

ScRNA-seq data were aligned to the mouse genome (version mm10), and gene counts were quantified using the default parameters of Cell Ranger (version 2.0.1, 10X Genomics). Following the Cell Ranger analysis, the gene-barcode matrices were processed using the Seurat R package (Butler et al., 2018) . To prepare the matrix for further analysis, we first excluded cells with less than 500 detected genes and cells with a mitochondrial gene ratio of more than 10%. Second, we used the-doubletFinder_v3‖ function in the DoubletFinder package (version 2.0.3) to exclude doublets from each sample, and the number of artificial doublets (pN) was set to 0.25 (McGinnis et al., 2019) . Third, the -NormalizeData‖ function was used to normalize the data using the -LogNormalize‖ method. Then, using -FindVariableFeatures‖ the top 2000 highly variable genes were identified. Finally, to mitigate the batch impact, we developed an integrated dataset using the -FindIntegrationAnchors‖ and -IntegrateData‖ functions. Following that, we used the -RunPCA‖ function to reduce the dimension of the combined after scaling and centering features. Cells were clustered at an appropriate resolution and visualized using a 2-dimensional t-distributed stochastic neighbor embedding (t-SNE) algorithm.

In this study, we used the abundance of known marker genes to annotate the cell types.

The number of cells of each kind in each group (Model, Con, and QFPDD) was tallied, as well as the percentage of each cell type. The cell ratio was then computed to determine the cell types that differed across groups. In this study, ratio C/M denoted the Control group divided by the Model group. If the ratio C/M is greater than 1, it indicates a reduction in cell numbers in the Model group and vice versa. Similarly, the ratio Q/M was also computed. As a result, -rescue cell types‖ were defined as the cell types that were changed in the Model group but were rescued in the QFPDD group.

The Wilcoxon rank-sum test was used to compare the expression of different groups (Model/Con and QFPDD/Model) using the ‗‗FindMarkers'' function of the Seurat package (version 3.2.3). Firstly, we identified differentially expressed genes (DEGs) across the Model and Con groups to generate Model-related DEGs (Model DEGs). We identified DEGs between the QFPDD and Model groups to generate QFPDD-related DEGs (QFPDD DEGs). ‗‗Rescue DEGs'' were defined based on the above results as the downregulated or upregulated genes among the Model DEGs that were upregulated or downregulated respectively, upon QFPDD. Similarly, for each cell type between different groups, differential expression analysis was also performed to generate cell-type-DEGs, and ‗‗cell-type-rescue DEGs''.

Metascape (Zhou et al., 2019) (version 3.5) (http://metascape.org) was used to perform GO biological process enrichment and pathway analyses with the ggplot2 R package used to visualize the results (p < 0.01). The enriched GO terms were shown for several groups/cell types.

SCENIC analysis was performed in Python (pySCENIC) using the GENIE3 and RcisTarget databases with raw count matrix as input (Aibar et al., 2017) . GENIE3 was used in a convenient way to construct gene regulatory networks from gene expression matrices. RcisTarget was used to identify potential target genes (regulons) and enrich transcription factor-binding motifs. The AUCell program was used to classify the activity of each regulon group within each cell. Additionally, the phyloseq package was used to determine regulon specificity scores for each cell type.

We downloaded a gene set named -HALLMARK INFLAMMATORY RESPONSE‖ from MSigDB and extracted cytokine genes based on these references to establish inflammatory and cytokine scores (Liberzon et al., 2011 ) (see Table   S1 -S2). Then, using Seurat's AddModuleScore function, we determined the inflammatory and cytokine scores. In this study, we identified the cell types as the most promising hyper-inflammatory cell type if the median value of its inflammatory score and cytokine score were both greater than zero.

CellChatDB was used to predict cell-cell interactions based on the expression of well-characterized ligand-receptor pairs in a variety of cell types (Jin et al., 2021) . To conclude, we loaded the counts' data into CellChat using the established procedure and utilized standard settings for the preprocessing functions -identifyOverExpressedInteractions‖, and -identifyOverExpressedGenes‖. We used the pre-compiled mouse Protein-Protein-Interactions as a prior network knowledge and selected the Secreted Signaling pathways for the database. The primary analyses were conducted using the core functions -computeCommunProbPathway‖, -computeCommunProb‖, and -aggregateNet‖ were used for the main analyses.

Finally, discrepancies in cell-cell interactions were determined using the functions -compareInteractions‖ and -netVisual_diffInteraction‖.

The untargeted metabolomics profiling of liver samples by gas chromatography-time-of-flight mass spectrometry (GC-TOF/MS) was carried out using the XploreMET platform (Metabo-Profile, Shanghai, China) based on previously reported methodologies with slight modifications (Qiu et al., 2014) . The metabolites were identified by comparison to known structures in the JiaLib metabolite database.

Data were shown as means ± SEM. Multiple comparisons were performed using one-way ANOVA and Dunnett's multiple comparisons test. The student's t test was used to compare between two groups. p < 0.05 was considered statistically significant.

We investigated the proportion of peripheral blood lymphocytes and lung inflammatory cytokines in mice pneumonia model to determine the effect of QFPDD on the immune system and inflammation. We observed a significant decrease in the proportion of peripheral lymphocytes (CD4 + T, CD8 + T, B cells) as well as a significant increase in lung TNF-α, IL-6, and IFN-γ in the Model group, indicating that the immune system was dysregulated in this pneumonia model. Importantly, we found that QFPDD significantly increased the percentages of CD4 + and CD8 + T lymphocytes in peripheral blood in mice in a dose-dependent manner, and increased the percentage of B lymphocytes in the QFPDD_M group compared to the Model group ( Figure 1A) , suggesting that QFPDD has an immune-enhancing impact.

Additionally, lung TNF-α significantly decreased in the QFPDD_L and QFPDD_H groups. IFN-γ also decreased in the QFPDD_M and QFPDD_H groups, while IL-6 only decreased in the QFPDD_L group compared to the Model group ( Figure 1B ),

indicating that QFPDD had an anti-inflammatory impact.

To obtain a better understanding of the effect of QFPDD on the makeup of cell types and transcriptional patterns in liver tissue, we created single-cell atlases for the Neutrophils cells (identified by S100a9, S100a8, G0s2 and Retnlg) and T cells (identified by Cd3e, Cd3g, and Cd3d) (Figure 2A-C) . In summary, we developed a tissue characterization of the cellular diversity of the liver, established a comparative framework, and developed a cellular strategy for future QFPDD research.

We independently analyzed the proportions of each cell type in each of the three groups to characterize the dynamics of cell-type composition. Figure S1A depicts the Single-Cell Atlases for each group. The cell ratio of various cell types in each group was analyzed and we identified the proportional variations in multiple cell types between the Model group and the Con group, and QFPDD was able to rescue many of these differences (Figures 2D-E) . For example, the Model group had a larger percentage of hepatocytes than the Con group, which was decreased by QFPDD.

Additionally, the percentage of HSC cells increased in the Model group compared to the Con group but decreased in the QFPDD group. QFPDD was also able to rescue the increased macrophage cells. T cells and B cells percentages decreased in the Model compared to Con but increased in the QFPDD group. These findings suggested that QFPDD is capable of reversing the populations of different types in mouse liver tissue and reconstructing the mice's cellular ecosystem.

To elucidate the molecular mechanism associated with QFPDD, hundreds of 16 DEGs were found between model and control mice and between QFPDD and model mice dubbed as ‗‗Model DEGs'' and ‗‗QFPDD DEGs,'' respectively. In this, ‗‗rescue DEGs'' were identified, which were defined as genes among model DEGs that were partially reversed by QFPDD. As shown in Figure 3A , the majority of DEGs have been rescued by QFPDD, and Figure 3B -C depicts the Venn diagram. Additionally, we examined two sets of genes: upregulated model DEGs that were elevated by QFPDD and downregulated model DEGs that were lowered by QFPDD. These DEGs accounted for around 1.2 percent of the model DEGs, showing that QFPDD had minor adverse effects in mice (Figures S1B-S1C) . Following that, we explored the biological implications of QFPDD using metascape (http://metascape.org/). The Following that, we assigned model, QFPDD, and rescue DEGs to each cell type to separate the effects of QFPDD. Surprisingly, as indicated by the rose diagrams, hepatocyte, macrophage, and neutrophil were the most substantially affected cell types in the liver tissues by Model and QFPDD, as well as rescue DEGs ( Figure   3F-H) , indicating that these cell types may be the primary target cells. Furthermore, enrichment analysis revealed that the DEGs of QFPDD were associated with metabolic pathways in hepatocyte cells, including cofactor metabolic process, amino acids, derivatives of metabolism, respiratory electron transport, and citric acid (TCA) cycle ( Figure 3I) . In macrophage cells, pathways associated with inflammation and innate immune response included angiogenesis, inflammatory response, and IL−17 signaling pathway ( Figure 3J ) and in neutrophil cells, pathways enriched with SRP-dependent cotranslational protein targeting the membrane, ribosome assembly, and regulation of cell−cell adhesion ( Figure 3K) . Notably, Alb, Trf, and Car3 were DEGs in all cell types affected by QFPDD, indicating that these proteins are critical in the therapeutic mechanism of QFPDD ( Figure 3L ).

Taken together, these findings, offer a broad characterization of QFPDD in mouse liver tissue at the single-cell level and demonstrate the value of large-scale, parallel gene expression studies for investigating complex biological events. These investigations gave insights into the cell-specific effects of QFPDD and may help paint a more comprehensive picture to gain a deeper understanding of QFPDD.

Transcription factors and their downstream-regulated genes constitute a complex and intertwined network of gene regulation that establishes and maintains cell identity.

We used single-cell regulatory network inference and clustering (SCENIC) to predict the core transcription factors in the liver tissues to further infer the activity of regulon (a TF together with its target genes comprise a regulon) for the most strongly affected hepatocyte, macrophage, and neutrophil. In this study, we created a regulon specificity score (RSS) based on Jensen-Shannon divergence and selected the top 15 regulons with higher RSS values as key transcription factors for each cell type. Rara, Jdp2, and Nfil3 were identified as the most specific regulons associated with neutrophils by network analysis (Figure 4A) . The TSNE plot adds to the evidence that the activities of these regulons are extremely specific to neutrophils ( Figure 4B) . We determined the ratio C/M of core transcription factors, which represents the control divided by the Model group. If the value of ratio C/M is greater than 1, it reflects a reduction in the activity of these transcription factors in the Model group and vice versa. Similarly, the ratio Q/M was calculated, and the findings revealed that the majority of the key transcription factors dysregulated during the Model group could be rescued by QFPDD in neutrophils ( Figure 4C) . The macrophage is another well-known cell type and the most specific regulons revealed by our network investigation were Gcm1 and Rxra ( Figure 4D-E) , and the regulatory activity was likewise rescued by QFPDD ( Figures 4F). The other, hepatocyte cell, having identified Nr1i2 and Foxa3 as the most specific regulons (Figures 4G-H) , most core-specific regulons could still be rescued by QFPDD (Figures 4I) . These findings reveal that QFPDD may reprogram 18 disease-compromised transcriptional regulatory networks and that QFPDD can restore the activity of transcription factors that have been deactivated by disease.

The proportions of the immune cells were then estimated, and the immune cells account for about half of the total, indicating a strong immune-inflammatory storm in the liver tissue. Based on the cytokine and inflammatory response genes we collected, we then looked into the likely causes of inflammatory storms. We first defined a cytokine and an inflammatory score for each cell type, then utilized these two interrelated scores as markers to assess each cell type's potential contribution to the inflammatory storm. We discovered that macrophages and neutrophils had seemingly increased expression of inflammatory and cytokine genes, and their scores were all greater than zero based on our scRNA-seq data ( Figures 5A-D) , suggesting that these cells were the likely primary sources of the inflammatory storm. According to previous research, neutrophils, which have been linked to the inflammatory response in COVID-19 patients, may influence platelet functioning at the disease stage (Ren et al., 2021) .

We then investigated the effects of QFPDD on the inflammatory storm and discovered that QFPDD significantly reduced the scores of hyper-inflammatory cell types ( Figure 5E-F) , indicating that QFPDD might reduce the inflammatory storm effect. Taken together, neutrophils and macrophages may be important mediators of the inflammatory storm in liver tissue, and QFPDD may considerably lessen the inflammatory storm phenomena of liver tissue.

The R package -CellChat‖ was used to thoroughly examine cell-cell communications between different cell types in different groups. We determined the intensity of ligand-receptor interactions in each pair of cell types in the Model, Con, and QFPDD groups. Computational analysis revealed an increase in the interaction intensity of cell-cell communication in the Model group compared with the Con group ( Figure 5G) . QFPDD, on the other hand, reversed several forms of cell-cell communication between certain cell types ( Figure 5H ). In particular, the interactions 19 between macrophage and neutrophils with other cell types in the Model group were significantly elevated and may be reversed by QFPDD, implying that QFPDD might modulate the interaction intensity of hyper-inflammatory cell types. Figure 5I -J shows the core pairings. Some pairs, such as Il34-Csf1r, Csf-Csf1r, Cxcl2-Ackr1, and

Osm-(Lifr+Il6st) were significantly reversed in neutrophils by QFPDD, whereas

Pdgfb-Pdgfrb, Lgals9-Cd45, Il34-Csf1r, Igf1-Igf1r, Hbegf-Egfr, and Csf1-Csf1r were significantly reversed in macrophages by QFPDD. These findings suggested that QFPDD reduced aberrant patterns of cell-cell communication, particularly in hyper-inflammatory cell types, demonstrating its therapeutic potential.

Because the liver is a key organ in maintaining body metabolic homeostasis, we used the GC-TOF/MS technique to examine the regulatory effect of QFPDD on hepatic metabolism in the Con, Model, and QFPDD groups. The metabolic profiles of the three groups showed good separation in the PLS-DA model (Figure 6A) ,

indicating that QFPDD influenced liver metabolism in this pneumonia model. There were 55 liver metabolites found between the Model and Con groups, and 13 metabolites identified between the Model and QFPDD groups (Figure 6B -C, Table   S3 , Table S4 ). Amino acid-related pathways (aminoacyl-tRNA biosynthesis pathway, arginine biosynthesis pathway, alanine, aspartate, and glutamate metabolism pathway)

were altered in the Model group ( Figure 6D) . Furthermore, the liver glutathione metabolism pathway was also dysregulated in the Model group (Figure 6D) , which plays a significant role in the antioxidation and detoxification processes of the liver.

The different metabolites between QFPDD and Model metabolites, on the other hand, were mainly involved in nucleotide metabolism (purine metabolism pathway, pyrimidine metabolism pathway) ( Figure 6E ). The relative intensity data revealed that numerous amino acids were downregulated in the liver of the Model group ( Figure 7A) . Furthermore, arginine biosynthesis metabolites (oxoglutaric acid, L-glutamine, L-citrulline, aspartic acid, ornithine), glutathione metabolism (ornithine, putrescine, spermidine), alanine, aspartate, and glutamate metabolism (L-glutamine, pyruvic acid, aspartic acid, oxoglutaric acid, succinic acid), arginine and proline 20 metabolism (ornithine, putrescine, spermidine, L-proline, pyruvic acid), D-Glutamine and D-glutamate metabolism (L-glutamine, oxoglutaric acid), TCA cycle (oxoglutaric acid, pyruvic acid, succinic acid), beta-Alanine metabolism (1,3-diaminopropane, aspartic acid, spermidine) were significantly altered in the Model group (Figure 7B) ,

suggesting that extensive metabolic disorders occurred in the liver of the Model group.

Furthermore, we found no significant restorative effect of QFPDD on these metabolites. However, purine metabolism metabolites in the Model group and QFPDD group were higher than in the Con group ( Figure 7C ).

HcoV-229E and SARS-CoV-2 belong to the same coronavirus family, and

HcoV-229E, which causes mild respiratory disease and is a major cause of common colds (10~30% of the cases), was first reported in the 1960s (Hamre and Procknow, 1966; McIntosh et al., 1967; van der Hoek, 2007) . Meanwhile, COVID-19 was first reported in a humid seafood market during winter, suggesting that a cold and humid environment may be favorable for COVID-19 development. HcoV-229E, combined with cold-damp environments, successfully induced pneumonia in BALB/c mice models (Bao et al., 2020a; Bao et al., 2020b) . These two factors were also shown to induce a typical increase in lung pro-inflammatory cytokine levels and a decrease in the abundance of immune cells in peripheral blood, clinical features that are comparable to those of COVID-19 patients Qin et al., 2020; Song et al., 2020; Wang et al., 2020a) . Clinically, QFPDD has beneficial immune modulatory effects (Xin et al., 2020) , such as improving the levels of total lymphocyte and white blood cell counts and mitigating the extent of multi-organ impairments in COVID-19 patients.

In previous study, we found that most of the research tends to stay in clinical retrospective study, network pharmacology and molecular docking, and there are not many experimental studies. In other words, the previous investigations were just mostly predictive and lacked experimental validation, especially the comprehensive analysis of the liver tissue. This study was the first report with direct experimental 21 evidence supporting the clinical benefits of QFPDD for improving the immune function of COVID-19 patients in the context of the pneumonia model of combined cold-damp environment exposure with HcoV-229E virus infection on mice. This study is also the first to integrate hepatic single-cell RNA sequencing and untargeted metabolomics into a TCM formula.

Firstly, QFPDD significantly increased the abundance of immune cells in the peripheral blood of mice with pneumonia, and decreased the levels of lung pro-inflammatory cytokines. These results prove the immune enhancement and anti-inflammatory effects of QFPDD. We further explored the underlying mechanisms of QFPDD using scRNA-seq analysis at a single-cell level. In this study, we used a total of 42,942 cells from three groups of liver tissues to investigate QFPDD-associated therapeutic mechanisms and biological changes. We constructed a single-cell atlas of liver tissues and identified 7 major cell types and gene expression signatures for each group. QFPDD was shown to reconstitute the cellular ecosystem, reverse gene expression changes as well as changes in cell-type-specific gene expressions. In addition, the activity of transcript factors was also recovered by QFPDD treatment, which could be the mechanism through which QFPDD leads to gene expression recovery in liver tissues. Moreover, QFPDD treatment reversed the aberrant cell-cell communication patterns. Figure 3L shows that Alb, Trf, and Car3

were the DEGs in all QFPDD-regulated cell types, implying that these proteins play a key role in therapeutic mechanisms of QFPDD. The most abundant protein in the blood is encoded by Alb (Albumin). This protein is involved in blood plasma colloid osmotic pressure regulation and serves as a carrier protein for various endogenous compounds, including metabolites, hormones, and fatty acids, as well as exogenous drugs. Alb is highly correlated with acute lung injury and, in this study, its expression levels were suppressed by QFPDD treatment Xia et al., 2020) . By interfering with antithrombin/SERPINC1-mediated inhibition of coagulation proteases, Trf (transferrin) promotes coagulation through a mechanism that is independent of its activity as an iron transporter. An elevated transferrin/antithrombin ratio has been associated with COVID-19-related coagulopathy and worsening disease 22 outcomes in older individuals, particularly in males (McLaughlin et al., 2020) .

Carbonic anhydrase isozymes are encoded by carbonic anhydrase III (Car3) gene, a multigene family member. These carbonic anhydrases are differentially expressed in many cell types and are a class of metalloenzymes that catalyze carbon dioxide reversible hydration. Evidence has shown that towards the renin angiotensin system (RAS), carbonic anhydrases have a major contribution and subsequently participate in SARS-CoV-2 pathogenesis in major organs such as renal, blood circulation and respiratory systems (Deniz et al., 2021; Reyes et al., 2020) .

The upregulated genes of QFPDD were enriched in SRP-dependent co-translational proteins, which targets the membrane, positive regulation of lymphocyte differentiation, lymphocyte activation, B cell differentiation, and CD4-positive, alpha-beta T cell differentiation. The downregulated genes of QFPDD were enriched in many metabolic pathways, including drug metabolic processes, citric acid (TCA) cycle, and respiratory electron transport. These results indicate that metabolic disorders and the immune microenvironment can be regulated by QFPDD.

The enriched primary functional pathways were strongly associated with the natural functions of the cells. Hepatocyte cells were associated with metabolic-related pathways, including cofactor metabolic processes, amino acids, and derivative metabolism, respiratory electron transport, and the citric acid (TCA) cycle; macrophage cells were associated with inflammation and innate immune responses, such as angiogenesis, inflammatory responses, and IL-17 signaling pathways;

neutrophil cells were related with SRP-dependent cotranslational protein targeting to membrane, ribosome assembly and cell-cell adhesion regulation. We also investigated cellular origins of a potential inflammatory storm in liver tissues. It was established that macrophages and neutrophils might be key sources of the inflammatory storm, and that QFPDD can significantly reduce the inflammatory storm effect.

On the other hand, amino acid metabolism was found to be altered in mice with pneumonia. The demand for glucose and amino acids has been shown to be increased in activated T cells due to proliferative expansion, and induction of a range of effector functions (Kedia-Mehta and Finlay, 2019). Therefore, the decrease in liver amino 23 acids in the Model group indicated an increase in amino acid expenditure, which might have been as a result of immune activation in the pneumonia model. In addition, the increase in liver pyruvic acid and oxoglutaric acid suggested enhancement of the TCA cycle and energy metabolism in this pneumonia model. In this study, we observed a significant increase in liver putrescine and spermidine levels in Model group, which have been respectively proven to inhibit and favor pro-inflammatory M1 macrophage polarization of (Latour et al., 2020) . Therefore, elevated liver putrescine and spermidine levels suggest a relationship between liver metabolic alterations and immune responses in the pneumonia model. Importantly, QFPDD mainly regulated liver purine metabolism related metabolites. Purine nucleosides, such as adenosine and its primary metabolite, inosine, belong to purine nucleosides and act as monomeric precursors for DNA and RNA. Adenosine is a powerful signaling molecule that primarily acts on its receptors, thereby playing a role in inflammatory regulation. Moreover, adenosine inhibits the synthesis of pro-inflammatory IL-12 and TNF-α cytokines (Haskó et al., 2000a) . Activation of adenosine receptors by agonists inhibits LPS-induced inflammation in vivo, and animals deficient in A2α adenosine receptor were established to be sensitive to sub-threshold doses of inflammatory stimuli (Ohta and Sitkovsky, 2001) . In summary, adenosine and its receptors are important in inflammatory attenuation. Inosine is transformed from adenosine via deamination, a process that mainly occurs during high intracellular concentrations of adenosine, which have been associated with ischemia, hypoxia, and other forms of cell stress (Deussen, 2000) . Inosine suppressed lymphocyte and neutrophil macrophage activation in vitro and counteracted the pro-inflammatory effects of LPS in vivo (Haskó et al., 2000b) . Moreover, it prevented systemic endotoxin-induced lung and gut injury as well as liver and vascular failure in mice (Soriano et al., 2001) . In this study, elevated liver adenosine and inosine levels in the Model group might be a stress adaptation characteristic to overcome inflammation in the liver, and QFPDD might protect the liver from injury by increasing adenosine and inosine levels.

Adenosine monophosphate (AMP) is an important metabolite that is also involved in purine metabolism, and can be generated through the dephosphorylation of adenosine 24 triphosphate (ATP) and adenosine diphosphate (ADP). Elevated cellular AMP/ATP ratio results in AMP-activated protein kinase (AMPK) activation, thereby regulating glucose and lipid metabolism in the liver. Thus, elevated AMP levels in the QFPDD group indicated the regulatory effects of QFPDD on liver metabolism through AMPK.

Liver malonic acid levels were also shown to be downregulated in the Model group, compared to the Con group, and QFPDD treatment significantly increased malonic acid levels, relative to the Model group. Malonic acid can be transformed into malonyl-coenzyme A (malonyl-CoA), which can be decarboxylated within the TCA cycle to acetyl-CoA and therefore fully oxidized. Besides, in mammalian metabolism, malonyl-CoA is positioned at a central regulatory node to coordinate fatty acid oxidation and synthesis. Thus, QFPDD might influence the TCA cycle and fatty acid metabolism pathways in the liver by elevating malonic acid levels. These findings were confirmed by scRNA-seq analysis data. Expressions of the key gene associated with metabolic pathways were also found to be regulated by QFPDD, which might result in enhanced immune functions.

In summary, in this paper, we for the first time integrated multidisciplinary technologies to explore the pharmacological mechanism of QFPDD. The novelty of our research rooted in providing an experimental and computational basis for further discovery of TCM formula and botanical drug for the treatment of COVID-19. These novel results may also guide future research on TCM in infectious diseases and hyperimmune-related diseases. Simultaneously, these findings provided key information and guidance for further investigation on the pharmacologically active substances and clinical applications of QFPDD.

This study is associated with some limitations. First, the relatively small number of single cell-sequences might affect the stability of the results. In future studies, we will expand the sample size. Second, more experimental studies should be performed to confirm the functional roles of the key genes. Despite the limitations, our findings elucidate on the effects of QFPDD on coronavirus-induced pneumonia in mice, and will form the basis for investigating the mechanisms of QFPDD. To the best of our knowledge, this is the first study to report the therapeutic mechanisms of QFPDD at 25 single-cell level.

In purine metabolism related metabolites (C). * p < 0.05, ** p < 0.01 compared with Con group. # p < 0.05, ## p<0.01 compared with Model group, based on Student's t test. 38 Blue arrows imply significant downregulation while red arrows imply significant upregulation in the Model group, compared to the Con group. Results are presented as means ± SEM.

SCENIC: single-cell regulatory network inference and clustering

Integrating single-cell transcriptomic data across different conditions, technologies, and species

Immune regulation by glucocorticoids

Specific ACE2 expression in cholangiocytes may cause liver damage after 2019-nCoV infection

Is carbonic anhydrase inhibition useful as a complementary therapy of Covid-19 infection?

Metabolic flux rates of adenosine in the heart

A new virus isolated from the human respiratory tract

Adenosine inhibits IL-12 and TNF-α production via adenosine A2a receptordependent and independent mechanisms

Inosine inhibits inflammatory cytokine production by a posttranscriptional mechanism and protects against endotoxin-induced shock

Characteristics of SARS-CoV-2 and COVID-19

Inference and analysis of cell-cell communication using CellChat

COVID-19 and the liver

Competition for nutrients and its role in controlling immune responses

The role of polyamines in the regulation of

Molecular signatures database (MSigDB) 3.0

Comprehensive profiling and characterization of the absorbed components and metabolites in mice serum and tissues following oral administration of

Clinical and biochemical indexes from 2019-nCoV infected patients linked to viral loads and lung injury

Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding

DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors

Recovery in tracheal organ cultures of novel viruses from patients with respiratory disease

COVID-19-related coagulopathy-is transferrin a missing link? 10

Role of G-protein-coupled adenosine receptors in downregulation of inflammation and protection from tissue damage

Dysregulation of Immune Response in Patients With Coronavirus

A distinct metabolic signature of human colorectal cancer with prognostic potential

Traditional Chinese medicine for COVID-19 treatment

COVID-19 immune features revealed by a large-scale single-cell transcriptome atlas

Contribution of hypoxia inducible factor-1 during viral infections

Single-Cell RNA-seq Reveals Angiotensin-Converting Enzyme 2 and Transmembrane Serine Protease 2 Expression in TROP2+ Liver Progenitor Cells

Lysoglycerophospholipids in chronic inflammatory disorders: the PLA(2)/LPC and ATX/LPA axes

Proteomic and Metabolomic 30 Characterization of COVID-19

Association between early treatment with Qingfei Paidu decoction and favorable clinical outcomes in patients with COVID-19: A retrospective multicenter cohort study

Omics-Driven Systems Interrogation of Metabolic Dysregulation in COVID-19 Pathogenesis

Inosine improves gut permeability and vascular reactivity in endotoxic shock

COVID-19 infection alters kynurenine and fatty acid metabolism, correlating with IL-6 levels and renal status. JCI Insight 5. van der Hoek, L

Novel Coronavirus-Infected Pneumonia in Wuhan, China

Regulation of macrophage function by sphingosine-1-phosphate

Preliminary exploration of the mechanism of Qingfei Paidu decoction against novel coronavirus pneumonia based on network pharmacology and molecular docking technology

An increased neutrophil/lymphocyte ratio is an early warning signal of severe COVID-19

Clinical retrospective study on the efficacy of Qingfei Paidu decoction combined with Western medicine for COVID-19 treatment

Investigating mechanism of Qing-Fei-Pai-Du-Tang for treatment of COVID-19 by network pharmacology

Metascape provides a biologist-oriented resource for the analysis of systems-level datasets

This study was supported by the National Key Research and Development Program of China (2020YFC0845400), the National Natural Science Foundation of China (82004215) and Shanghai Municipal Health Commission Project (20204Y0326).We express great gratitude to Professor Xiaolan Cui and her team at the China Academy of Chinese Medical Sciences (ABSL-2 laboratory) for conducting animal and inflammatory cytokine assays in lung tissues as well as for detecting peripheral blood lymphocytes via flow cytometry. In addition, we also thank OE Biotech Co.Ltd (Shanghai, China) for providing single cell RNA-seq and the Home for Researchers editorial team (www.home-for-researchers.com) for the English check of this manuscript.

Li, Gaosong Wu, Jingyan Han revised the manuscript. Jing Zhong and Lili Sheng helped in improving the language. The order of these co-authors is based on the length of time spent on the project. All data were generated in-house, and no paper mill was used. All authors agree to be accountable for all aspects of the study.