Whole Transcriptome Analysis on Sepsis Caused by Pyogenic Liver Abscess

Background Pyogenic liver abscess (PLA) is a common clinical disease throughout the world with remaining high mortality and morbidity. This study sought to investigate key genes and non-coding RNAs (ncRNAs) associated with the pathogenesis and treatment of sepsis caused by PLA (SLA). Methods We constructed the expression proles of RNAs by RNA sequencing in the peripheral blood of healthy subjects, SLA and seven day therapy of SLA patients. Then the sequencing data was analysed through integrated bioinformatics strategy. Results 4549 mRNAs, 9808 lncRNAs, 467 circRNAs and 292 miRNAs were differentially expressed (DE) between SLA and control samples. 3199 mRNAs, 6018 lncRNAs, 161 circRNAs and 132 miRNAs were differentially expressed between SLA after 7 days therapy and SLA control subjects. Moreover, DE mRNAs of the former two comparisons were both related to immunity and T cell-mediated immune response; DE lncRNAs were both associated with B cell receptor signaling pathway and Osteoclast differentiation; DE circRNAs were both related to Chagas disease and B cell receptor signaling pathway; DE miRNAs exerts its specic role in the process of SLA by MAPK signaling pathway and Estrogen signaling pathway compared to DE lncRNAs, circRNAs and mRNAs. Finally, lncRNA/circRNA–miRNA–mRNA network was constructed to explore the internal regulating relationships of these ncRNAs and mRNAs.


Introduction
Pyogenic liver abscess (PLA) is a common infectious disease of the liver, which can be divided into amoebic, bacterial, fungal and hydatid liver abscesses [1]. Among them, bacterial liver abscess is the most common, accounting for 80% of the incidence of liver abscess [2]. The incidence of PLA has been on the rise in recent years. PLA has a rapid onset, severe disease and rapid development [3]. The mortality rate has increased signi cantly after the progression to sepsis [4]. There are few studies on PLA-induced sepsis, and there is a lack of biomarkers for early recognition of liver abscess progression to sepsis.
Non-coding RNA (ncRNA) cannot be translated into protein, but it can regulate protein expression [5].
MicroRNAs (miRNAs) are a family of ncDNAs with 21-25 nucleotides, which play an important role in regulating protein expression by binding to speci c target mRNA for cleavage or translation inhibition [6].
Long noncoding RNA (LncRNA) is a kind of ncRNA that accounts for the majority of mammalian ncRNAs, with a length of more than 200 nucleotides, and regulates gene expression through a variety of mechanisms [7]. Circular RNA (circRNA) is an endogenous ncRNA with a closed circular structure, which can regulate linear RNA transcription, downstream gene expression and protein production [8]. Some researchers have shown that mRNA, miRNA, lncRNA and circRNA can interact as competitive endogenous RNA to form ceRNA regulatory networks [9,10]. LncRNA and circRNA can act as ceRNA molecules and have the function of miRNA sponge [11]. They regulate gene transcription and translation by competitively binding to RNA binding protein sites on miRNA, which plays a role in multiple physiological and pathological processes [12]. However, so far, SLA-related RNA molecules are still rarely studied in SLA, and their speci c role in the pathogenesis and treatment of SLA is not clear. Therefore, comprehensive detection and analysis of ncRNAs in the development of SLA are critical for the prevention and treatment of SLA.
In the present study, we used whole transcriptome sequencing analyses to construct RNA expression pro les in peripheral blood of samples. Two comparisons (SLA samples vs. healthy controls, 7 days therapy of SLA vs. SLA) were applied to screen differentially expressed mRNAs, miRNAs, lncRNAs, and circRNAs, and their functions were predicted by bioinformatics analysis. The analysis ow of this study is shown in Fig. 1. The ndings facilitate our understanding of regulatory mechanisms, discovering new targets and providing assistance for the treatment of sepsis caused by liver abscess.

Ethical approval.
This research was approved by the Ethics Committee of the rst hospital of Jilin University (20K019-002).
All subjects, enrolled from the rst hospital of Jilin University, agreed to participate in the study after fully understanding the purpose and procedure of the experiment. Full inclusion and exclusion criteria are listed in Table 1. The subjects were divided into three groups: SLA (M0), SLA after 7 days therapy (T7) and the healthy subjects (C). Subsequently, the fasting venous blood of all subjects were obtained and stored at − 80°C until analysis.

RNA isolation, library preparation & RNA sequencing
Total RNA was extracted using Trizol reagent (Life Technologies, MD, USA), and the quality and quantity were assessed by real-time quantitative PCR (RT-PCR) and Agilent 2100 bioanalyzer (Thermo Fisher Scienti c, MA, USA). Ribo-Zero™ rRNA Removal Kits (Epicentre, WI, USA) were applied for rRNA depletion. A total 1 µg RNA of each sample was used for constructing cDNA library. This library was sequenced by Illumina HiSeq4000 sequencer. The raw reads were saved as FASTQ format and processed for quality control. After removing low-quality reads and adaptor sequences, clean data were obtained and used for all subsequent analyses. The fragments per kilobase per million values were used to estimate the level of gene expression. Differential expression of circRNAs, mRNAs, miRNAs and lncRNAs with statistical signi cance (DEGs) was identi ed by DESeq package, based on cut-off criteria of fold change (FC) > 2 and corrected p-value < 0.05.

Gene ontology and pathway enrichment analyses.
To understand the biological signi cance of these DE transcriptome, we performed functional enrichment (KEGG pathway and Gene ontology enrichment) analyses through the clusterPro ler package [13]. And those with P < 0.05 were signi cantly enriched. R language ggplot2 and GOplot package were used to draw the bubble plot of the top 5 enriched GO terms in each category (biological process, molecular function, and cellular component).

PPI network construction and sub-network identi cation.
To discovery the co-differentially genes in the two comparisons (SLA samples vs. healthy controls, 7 days therapy of SLA vs. SLA), we conducted venn analysis of the differential expressed mRNAs (http://bioinformatics.psb.ugent.be/webtools/Venn/). Then, the Search Tool for the Retrieval of Interacting Genes (STRING) database [14] was used to construct the protein-protein interaction (PPI) network of the common and opposite differential expressed mRNAs between the two comparisons. The interaction pairs with the PPI combined score > 0.4 were selected as signi cant and the resulting network was performed with Cytoscape(version 3.2.1) [15]. The MCODE can effectively nd densely connected regions of a molecular interaction network [16]. MCODE was used to discovery the PPI network function module, which used a cut-off value for the connectivity degree of nodes (proteins in the network) greater than 3 and score > 4.

ceRNA analysis.
Cytoscape was used to construct the lncRNA/circRNA-miRNA-mRNA ceRNA regulatory network based on the DE transcriptome in the blood between SLA vs. healthy samples and 7 days therapy of SLA vs. SLA patients. The screening criteria of ceRNA were as follows: 1. differentially expressed mRNAs, lncRNAs, circRNAs and miRNAs; 2. a targeted relationship and negative correlation between miRNA and lncRNA/circRNA (lncRNA-miRNA-mRNA: pvalue < 0.05 and cor ≤ -0.7, circRNA-miRNA-mRNA: pvalue < 0.05 and cor ≤ -0.6); 3. a targeted relationship and negative correlation between miRNA and mRNA (pvalue < 0.05 and cor ≤ -0.7); 4. a signi cant correlation between lncRNA and mRNA targeted to the same miRNA (lncRNA-miRNA-mRNA: pvalue < 0.05 and cor ≥ -0.7, circRNA-miRNA-mRNA: pvalue < 0.05 and cor ≥ 0.6). Then, to discovery the co-differentially lncRNA/circRNA-miRNA-mRNA in the two comparisons (SLA samples vs. healthy controls, 7 days therapy of SLA vs. SLA), we conducted venn analysis of the differential expressed ceRNA regulatory network. Finally, ceRNA was further constructed using the two intersections (up-regulated ceRNA in SLA patients and down-regulated after therapy, downregulated ceRNA in SLA patients and up-regulated after therapy) 2.6 Statistical Analysis.
This study used SPSS software (version 20.0) and GraphPad Prism 8 (GraphPad Software, CA) for the statistical evaluations. The values were expressed as mean ± standard deviation. Statistical differences were performed with two-tailed Student's t-tests (two groups) and One-way ANOVA (three groups). Differences with p ≤ 0.05 were indicated signi cant.

Clinical characteristics of the participants.
In this study, eight healthy subjects, eight SLA patients and eight SLA patients after 7 days therapy were enrolled. All patients ful lled the diagnostic criteria for SLA. There was no signi cant difference in age and gender between the SLA and healthy control groups. The characteristics of all subjects are shown in Table 1. Compared with the SLA patients, the SLA patients after 7 days therapy exhibited signi cant increases in platelets level and decreases in procalcitonin and CPR levels.

Quality assessments and mapping results
To construct expression pro les of miRNAs, lncRNAs, circRNAs, and mRNAs of the three groups, transcriptome data sets were generated by whole-transcriptome sequencing. Subsequently, quality control and mapping analysis were performed for the sequencing output (Supplementary Table 1). For miRNA clean data, the base percentage of Q20 was > 97%, Q30 was > 92.3%. For total RNA clean data, the base percentage of Q20 was > 97.2%, Q30 was > 92.3%, and the average mapping rate of the 24 samples was 97.73%. These results could indicate the high quality of transcriptome sequencing data with suitable mapping. of hierarchical clustering of DE mRNAs, lncRNAs, circRNAs or miRNAs was generated to visualize the overall pattern of gene expression (Fig. 2). Moreover, the PCA loading plot showed that DE mRNAs, lncRNAs, circRNAs or miRNAs can be independently used to distinguish SLA after 7 days therapy samples, SLA and healthy control samples (Fig. 3).

Functional and pathway enrichment analyses of DEGs.
For the DE transcriptome in SLA vs. healthy samples, the GO terms (Fig. 4) showed that DE mRNAs were mainly related to neutrophil activation, neutrophil mediated immunity, and secretory granule membrane et al. (Supplementary Table 2); DE lncRNAs were signi cantly related to neutrophil activation, secretory granule membrane, and transcription factor complex et al. (Supplementary Table 3); DE circRNAs were signi cantly related to ubiquitin-like protein transferase activity, secretory granule membrane, and protein polyubiquitination et al. (Supplementary Table 4); DE miRNAs were signi cantly related to vacuolar membrane, NAD-dependent histone deacetylase activity, and acetylcholine receptor binding et al. Table 5). On the other hand, the top 5 enriched KEGG pathways (Fig. 5) Table 9).

(Supplementary
For the DE transcriptome in 7 days therapy of SLA vs. SLA patients, the GO terms showed that DE mRNAs were mainly related to neutrophil activation, neutrophil mediated immunity, and secretory granule membrane et al.; DE lncRNAs were signi cantly related to methyltransferase activity, membrane lipid metabolic process, and transcription factor complex et al.; DE circRNAs were signi cantly related to amyloid precursor protein metabolic process, secretory granule membrane, and cadherin binding et al.; DE miRNAs were signi cantly related to vacuolar membrane, NAD-dependent histone deacetylase activity, and regulation of metal ion transport et al. (Fig. 4). On the other hand, the top 5 enriched KEGG pathways terms of DE mRNAs were primary immunode ciency, osteoclast differentiation, hematopoietic cell lineage, NF-kappa B signaling pathway and PD-L1 expression and PD-1 checkpoint pathway in cancer; DE lncRNAs were signi cantly related to Glycosphingolipid biosynthesis-lacto and neolacto series, B cell receptor signaling pathway, central carbon metabolism in cancer, Glycosaminoglycan biosynthesiskeratan sulfate, and various types of N-glycan biosynthesis; DE circRNAs were signi cantly related to transcriptional misregulation in cancer, proteoglycans in cancer, PD-L1 expression and PD-1 checkpoint pathway in cancer, Legionellosis, and axon guidance; DE miRNAs were signi cantly related to 2-Oxocarboxylic acid metabolism, oxytocin signaling pathway, glycosaminoglycan biosynthesischondroitin sulfate, human immunode ciency virus 1 infection, and GnRH signaling pathway (Fig. 5).

Construction and analysis of PPI networks.
PPI networks were constructed through a STRING database to explore the functional relationship between the SLE therapy and SLE differential genes. Cytoscape was used to analyze PPI networks of the 364 dysregulated genes, including 243 U2D and 121 D2U genes. The results showed that there were 269 nodes and 696 edges (Fig. 6). MCODE was used to discovery the PPI network function module and four signi cant modules (nodes > 3 and score > 4) presented (Fig. 7). Module 1 (score: 12.5) consisted of 13 nodes and 75 edges. The seed gene was CYSTM1; Module 2 (score: 7.931) consisted of 30 nodes and 115 edges, seed gene was IRAK2; module 3 (score: 4.5) consisted of 5 nodes and 9 edges, and the seed gene was S100A8; Module 4 (score: 4.5), consisted of 5 nodes and 9 edges, and the seed gene was ACER3. These four seed genes were all up-regulated in SLE compared to healthy controls, and downregulated after therapy.
As the method described, we constructed the ceRNA regulatory network of M0 vs. C and T7 vs. M0, respectively. The lncRNA-miRNA-mRNA network contained 68275 pairs in M0 vs. C and 1632 pairs in T7 vs. M0. A total of 1031 common and opposite pairs were screened (Fig. 8). Speci cally, 459 ceRNAs were upregulated in M0 vs. C (lncRNA_down-miRNA_up-mRNA_down) and down-regulated in T7 vs. M0

Discussion
In this study, we used whole transcriptome sequencing to investigate the signi cantly differentially expressed mRNAs, circRNAs, lncRNAs, and miRNAs in the PBMC of patients with SLA. Then, these differentially expressed RNAs were to predict their potential biological functions through GO and KEGG pathway analysis. In addition, we also performed PPI network and circRNA/lncRNA-miRNA regulatory network. Our results showed that ncRNAs may participate in the pathogenesis of SLE and provide some potential targets for the treatment of SLA.
In this study, 243 genes were upregulated in SLA patients and downregulated after seven days therapy, and 121 genes were downregulated in SLA patients and upregulated after seven days therapy. And the common GO terms of these two comparisons included neutrophil activation, neutrophil mediated immunity, T cell activation, et.al; the common KEGG terms contained NF-kappa B signaling pathway, Th1 and Th2 cell differentiation, T cell receptor signaling pathway. Consistent with our results, previous study has found that circulating Th1 and Th2 subset accumulation kinetics was related to immunosuppression and prognosis in sepsis [17]. It is known that sepsis is the result of host immune response to infection, and the original site of infection might be the main driver of secondary organ injury or even death [18]. Moreover, another study shows that circulating Th1 and Th2 subset accumulations differ in septic patients with different infection sites [19]. Therefore, these results indicated that the pathogenesis and treatment of SLE may be related to immunity and T cell-mediated immune response.
LncRNA, a mRNA-like transcript, functions in various biological process [20]. Through its colocated and coexpressed protein-coding genes, we predicted the biological function of lncRNA. The common GO terms of the two comparisons (T7 vs. M0 and M0 vs. C) included transcription factor complex and nuclear transcription factor complex; the common KEGG terms contained Chagas disease, B cell receptor signaling pathway, Osteoclast differentiation, RNA degradation and Central carbon metabolism in cancer. Previous study has demonstrated that B-cell receptor-associated protein 31 (BAP31) played a critical role in mitochondrial function and ER homeostasis in sepsis-related myocardial injury [21]. A transcriptomic meta-analysis study reveals that osteoclast differentiation was upregulated in human septic shock [22]. Therefore, DE lncRNAs appear to function in the pathogenesis and therapy of SLA by regulating B cell receptor signaling pathway and Osteoclast differentiation. In our future work, the speci c role of lncRNAs of these two functions needs to be further studied.
CircRNA expression was signi cantly altered in PBMCs of sepsis patients [23]. To identify a benchmark map of circRNA expression dynamics in the therapy of SLA patients, we identi ed DE circRNAs in two comparisons (T7 vs. M0 and M0 vs. C). The common GO terms were secretory granule membrane and inclusion body; the common KEGG terms contained Chagas disease, B cell receptor signaling pathway, Axon guidance, Transcriptional misregulation in cancer, PD-L1 expression and PD-1 checkpoint pathway in cancer, and Proteoglycans in cancer. Interestingly, Chagas disease and B cell receptor signaling pathway were signi cantly dysregulated between lncRNA and circRNA functional enrichment, suggesting that these two pathways may function in SLA through lncRNA and circRNA.
DE miRNAs in the therapy of SLA patients were performed by functional enrichment analyses of predicted target genes. The common GO terms of the two comparisons (T7 vs. M0 and M0 vs. C) included vacuolar membrane, NAD-dependent histone deacetylase activity and NAD-dependent protein deacetylase activity; the common KEGG terms contained MAPK signaling pathway, Rap1 signaling pathway, Oxytocin signaling pathway, Estrogen signaling pathway, GnRH signaling pathway et.al.. MAPK, a serine-threonine kinase, plays a critical role in regulating stress response, in ammation, cell proliferation and apoptosis [24]. A recent study showed that inhibition of MAPK pathway could ameliorate sepsis-induced cardiac injury in AT1R-knockdown rats [25]. It is known that estrogen counteracts in ammation and clinical symptoms in various diseases including chronic liver injury [26], neurovascular/neuroimmune disease [27], and aging [28]. A previous study revealed that estrogen could protect sepsis-related disorder through decreasing in ammation [29]. Moreover, a recent study proved that estrogen signi cantly ameliorated sepsis-induced liver injury through blocking oxidative stress mediated activation of pyroptosis signaling pathway [30]. Interestingly, these common functional terms were different from those of lncRNA, circRNA and mRNA. Taken together, the results revealed that DE miRNAs exerts its speci c role in the process of SLA by MAPK signaling pathway and Estrogen signaling pathway compared to DE lncRNAs, circRNAs and mRNAs.
In recent years, the role of ceRNA in sepsis has become increasingly prominent. Fang et al. proved that lncRNA-H19 increase the expression of AQP1 through competitive binding miRNA-874 in sepsis, restoring dysregulated in ammatory responses [31]. Circular RNA VMA21 ameliorated sepsis-associated acute kidney injury via regulating miR-9-39/SMG1/in ammation axis and oxidative stress [32]. In this study, lncRNA/circRNA-miRNA-mRNA network was also constructed to discovery the regulatory relations of these ncRNAs and mRNAs. In the lncRNA-miRNA-mRNA network, miR-150-5p was downregulated in SLA patients compared to healthy control, and recover after 7 days therapy. Consisting with our results, Zhu et al. [33]found that miR-150-5p protected sepsis-induced myocardial depression by alleviating apoptosis. Interestingly, miR-150-5p show the same trend in circRNA-miRNA-mRNA network, suggesting that miR-150-5p may have a key role in the therapy of SLA through lncRNA/circRNA-miRNA-mRNA network. In addition, many researches revealed several miR-370-3p related ceRNA in sepsis progression, including lnc-MALAT1/miR-370-3p/HMGB1 [34], lnc-NEAT1/miR-370-3p/ Irak2 [35], lnc-NEAT1/miR-370-3p/TSP-1[36], which is consistent with our results of circRNA-miRNA-mRNA network. However, the mechanism of ncRNAs is complicated and we will further explore these predicted ceRNA mechanisms. Finally, it will help to discovery the therapeutic targets and internal mechanism of SLA.

Conclusion
In conclusion, we used whole transcriptome sequencing analysis to explore the expression pro le of RNAs in the peripheral blood of seven days therapy of SLA, SLA and healthy subjects. DE mRNAs, lncRNAs, circRNAs and miRNAs were screened, and their predicted biological functions were performed via bioinformatics analysis. These results may provide new perspectives on the pathogenesis of SLA and accelerate the development of new therapeutic approaches on SLA.  Analysis ow chart of this study Figure 2 The expression pro ling changes of RNAs in peripheral blood sample of SLA patients. The heatmap plot of differentially expressed mRNAs (A), lncRNAs(B), circRNAs (C) and miRNAs (D).