INTRODUCTION

Acute lung injury (ALI) is one of the severe clinical complications distinguished by pulmonary vascular permeability [1]. Meanwhile, ALI exhibits predisposing and precipitating factors, such as acute pneumonia, acute pancreatitis, and severe trauma, and acute pancreatitis, and ultimately develops into acute respiratory distress syndrome (ARDS) [2]. To make matters worse, ALI and ARDS are linked to many destructive clinical disorders with high morbidity and mortality rates [3]. Lipopolysaccharide (LPS) is the main ingredients of the cell wall in gram-negative bacteria which acts as the pathogenic infection resulting in ALI [4]. Triggered by LPS stimulation, inflammatory cytokines, such as tumor necrosis factor-α (TNF-α) and interleukin-1β (IL-1β), might facilitate ALI development [5, 6] Tremendous core genes have been characterized to date and deepen our understanding of the pathogenesis of ALI [7,8,9]. Nevertheless, the immunological mechanism is still unclear.

Stunning advances in high-throughput technologies and bioinformatic mining have produced huge volumes of data. Conventional gene expression analysis generally focuses on a specific gene, discounting the effect of interactions among individual genes. Thereby, efficient and robust tools are the key to integrate available data and provide insight into complex diseases. Recently, a comprehensive approach is constructed using protein–protein interaction (PPI) data to interpret the interactive patterns across multiple datasets. This network meta-method has been developed as a feasible way to analyze large gene profiles and successfully applied in the health bioinformatics field [10,11,12].

To characterize hub genes and explore the molecular mechanisms of LPS-treated ALI, we built a compendium of genes related to ALI using gene expression profiles from six microarray studies. We then used the above PPI networked comprehensive analysis to detect the hub genes and functional modules. The ALI mice model was performed by intratracheally instilling with LPS. Seven candidate genes were confirmed using an established mouse model of ALI and qRT-PCR validation. The results pave the way for future research into the pathophysiology of ALI, which might influence far-reaching individual anti-ALI treatment. Despite a lot of putative biomarkers were yielded by screening a large amount of data, determining the functional effect of each gene will require in-depth verification in animal models and clinical samples.

MATERIAL AND METHODS

Microarray Data Processing

To mine the related genes of LPS-induced ALI in mice, gene expression microarray studies were derived from the GEO (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo) repository by following strategies: (1) keyword search: “LPS-induced ALI” and “Mus musculus,” “lung,” and “expression profiling by array”; (2) it contained at least three specimens per group. Consequently, six available profiles (GSE102016, GSE2411, GSE16409, GSE104214, GSE17355, and GSE18341) met the inclusion criteria and were selected for the next analyses [13,14,15,16,17,18]. Detailed information (GEO ID, platform information, and sample accession) of these datasets is summarized in Table 1.

Table 1 Characteristics of Six Eligible GEO Studies Composing the Compendium
Full size table

Integrated Multi-microarray Analysis

The integrated microarray analysis was conducted through a visual analytics platform NetworkAnalyst 3.0 [19]. Firstly, all gene probes were converted to Entrez ID. Secondly, we preprocessed the expression profiles by log2 transformation and auto-scaling before quantile normalization. For the meta-analysis based on NetworkAnalyst 3.0, three analytics methods were employed including the Fisher method, fixed-effect model (FEM), and vote-counting (combined p-value < 0.05 or vote counts ≥ 2 were considered to be significant) after adjusting study batch effect. Combined p-value < 0.05 or vote counts ≥ 2 were considered to be significant. In FEM, the estimated effect size (ES) in each group is assumed to come from an underlying true effect size plus measurement error. Finally, the common DEGs defined by above three methods were regarded as the final candidate genes.

Functional Enrichment Analyses of DEGs

According to the hypergeometric distribution algorithm, GO (Gene Ontology, http://www.geneontology.org/) enrichment analyses divided into biological process (BP), molecular function (MF), and cell component (CC) were implemented using “enrichGO” function of “clusterProfiler” R package [20]. To simplify enriched GO results, we then applied “enricher” function to reduce redundant GO terms (The threshold parameter was set by default to “cutoff = 0.7, by = “p.adjust”, select_fun = min”). Furthermore, the pathways of the identified proteins were classified via “enrichKEGG” function of “clusterProfiler” package for KEGG (Kyoto Encyclopedia of Genes and Genomes) annotation.

To test a set of related genes in coordinated way, GSEA (Gene Set Enrichment Analyses) of GO and KEGG were enforced by “gseGO” and “gseKEGG” functions of “clusterProfiler” R package, separately. The Hallmarks gene set was downloaded from MSigDB v7.1 (Molecular Signatures Database, https://www.gsea-msigdb.org/gsea/msigdb/index.jsp). P-value < 0.05 were considered as statistically significance.

Network-Based Meta-analysis and Extracting Co-expressed PPI Modules

A network-based analysis was performed using STRING search function in NetworkAnalyst website (parameters were set to “confidence score cutoff = 900” and “require experimental evidence”). The extended network construction was restricted to include only the original seed proteins. Thus, those with a high level of confidence were retained and regraded as valid interactions. Besides, PPI networks of core seven DEGs were predicted via the STRING v10.5 (Search Tool for the Retrieval of Interacting Genes/Proteins, http://www.string-db.org/). Proteins that linked to each other were detected based on corresponding gene neighborhood, co-occurrence, co-expression, protein homology, experimental determination, curated databases, and text mining.

The Murine Model Establishment of ALI and qRT-PCR Validation

The regents used and the detailed procedures of the murine model establishment and qRT-PCR were performed as before [8]. The H&E immunohistological staining images of the control group and the ALI group are shown in Figure S1. All experiment protocols conformed to the guidelines of the China Council on Animal Care and Use. These animal studies were approved by the Institutional Animal Research Committee of Union Hospital. Mice were randomly allocated into ten control samples and ten LPS samples and then feeding a week of mice to adapt the environment. Control mice were exposed to an aerosol of phosphate buffer saline (PBS) alone. For acute LPS exposure, mice were exposed to PBS containing 0.5 mg/mL LPS for 2 h, in a custom-built cuboidal chamber. The LPS solution was aerosolized with a constant output ultrasonic nebulizer (model: 402B, Yuwell, China) at a flow rate of 30 ml/h. LPS was purchased from Sigma-Aldrich (extracted from Escherichia coli O55: B5, L2880). The chamber was 18 cm long, 12 cm wide, and 12 cm high. The 2−ΔΔCt method was used, and the PCR reactions were normalized to the GAPDH gene. All primers were sourced from PrimerBank (https://pga.mgh.harvard.edu/primerbank/). The reverse transcription step of qRT-PCR was performed by the UEIris RT mix with DNase (All-in-One) (US Everbright® Inc., Suzhou, China). The qRT-PCR reaction (10 µL) was formulated by the 2 × SYBR Green qPCR Master Mix (US Everbright® Inc., Suzhou, China). All qRT-PCRs were carried out on a CFX96 Touch™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The primer sequences are provided in the Supplementary Materials (Table S1).

Statistical Analysis

The Benjamini–Hochberg application of the FDR (false discovery rate) was applied to correct the p-value. Significantly enriched GO terms and KEGG pathways were distinguished using hypergeometric tests. Statistical analysis of significant differences between groups in qRT-PCR was complied with one-way ANOVA by Prism 7 software.

RESULTS

A Compendium of the Eligible Microarrays in This Study

Under our inclusion criteria, a gene expression compendium was finally built using six independent studies (GSE18341, GSE104214, GSE17355, GSE102016, GSE16409, and GSE2411) from the GEO database. In GSE18341, juvenile (21 day) and adult (16 week) C57BL/6 mice were treated with inhaled LPS (20 mL of 0.1 mg/mL) for 30 min in a sealed aerosol chamber and then spontaneous breathing. Untreated age-matched controls acted as comparison groups. In GSE104214, the experimental group was treated with 5 μg of intranasal LPS at 8 and 24 h. Controls received intranasal PBS. As for GSE17355, injury peaked at day 4 and is almost completely resolved by day 10 in wild type C57BL/6 mice in LPS model of ALI. Total RNA was isolated from mouse lung at times 0, 1, 4, and 10 days following LPS treatment of wild type. In GSE102016, mice were intranasally treated with LPS (20 μg/mouse) to induce pulmonary inflammation and all lung samples were collected after 72 h. In GSE16409, samples were from control or LPS exposed (1.5, 6, and 12 h post exposure). Twenty-four C57/B6 male mice were randomized to control and LPS groups in GSE2411. A total of 47 LPS-challenged and 29 normal mice were included in this study (Table 1). Then, we extracted and annotated these six GEO data, yielding a compilation of 2030 genes from 76 samples totally. The entire work flow in this study is depicted in Figure S2.

Overlap of Differentially Expressed Genes Among Datasets

A comprehensive analysis across studies was performed by computing the DEGs per dataset and assessing the overlap of the significant results. As shown in Fig. 1A, 3400, 1392, 653, 194, 94, and 13 DEGs were generated in GSE18341, GSE104214, GSE17355, GSE102016, GSE16409, and GSE2411, respectively. Besides, seven overlapped genes (Ifi44, Tnip1, Oasl1, Casp4, Ccl12, Zbp1, and Cxcl13) were tightly associated with LPS among at least five datasets (Fig. 1A). Among them, caspase-4 (Casp4) acts as an innate immune receptor for cytoplasmic LPS [21]. And the inflammatory cytokines and chemokines Ccl12 and Cxcl13 have been reported to participate in the pathogenic mechanism of acute lung injury (ALI) [22]. These findings implied that these three genes in the five datasets were closely related with the LPS host response in mice.

Fig. 1
figure 1

Comprehensive analysis of six GEO compendium. (A) Significant ALI-associated genes in six independent studies of the ALI compendium. The upset plot showed the DEG count in each group. The black dots and lines between the dots represented DEGs in which GEO dataset. The bars represented number of DEGs in the group. (B) Venn diagram of DEGs identified from the meta-analysis using Fisher’s method, the vote counting method, and a fixed-effect model. (C) Heat map representation of the top 30 upregulated DEGs across different microarrays identified from the meta-analysis (row-wise comparison). The color of boxes varying from blue to red represented the combined effect size value of DEG in each sample.

Full size image

Comprehensive Analysis of GEO Microarrays

Next, we performed meta-analyses through NetworkAnalyst 3.0 platform. The merged data of this analysis is tabulated in Table S2. After employing three meta-analysis methods, we validated 1703, 1316, and 996 DEGs by the Fisher method, FEM, and vote counting, respectively. Of those, 958 DEGs were identified via all three methods (Fig. 1B; Table S3). Thereinto, mRNA levels of 470 (49.1%) DEGs elevated and 488 (50.9%) DEGs were lower in LPS group in comparison with the control. Two heat map visualization of the top 30 up- and downregulated DEGs across the different studies is displayed in Fig. 1C and Figure S3, respectively. In the overlapping DEGS across the six microarray cohorts, Ebi3 was the top upregulated gene followed by F10, and Fmo3 was the most prominently downregulated gene.

Functional Enrichment and GSEA Analyses of DEGs

To exploit potential features of DEGs, we analyzed GO function and KEGG pathway enrichment. The top 5 most enriched KEGG pathways of DEGs were involved in osteoclast differentiation, fc gamma R-mediated phagocytosis, MAPK signaling pathway, fluid shear stress and atherosclerosis, and leishmaniasis (Fig. 2). Moreover, we were able to obtain a global perspective of the changes in gene expression patterns In GO enrichment analysis (Table S4). In MF term, DEGs were centered on “ubiquitin-like protein ligase binding,” “coenzyme binding,” and “GTPase activity.” As for the BP category, the core DEGs were significantly enriched related to “leukocyte proliferation,” “cytokine-mediated signaling pathway,” and “response to oxidative stress.” In addition, DEGs were enriched in the CC category focused on “membrane region,” “actin cytoskeleton,” and “NADPH oxidase complex” (Fig. 3).

Fig. 2
figure 2

Top 5 KEGG pathway analysis results of DEGs. The color of links represented different KEGG categories. Red and green dots showed up- and downregulated DEGs.

Full size image
Fig. 3
figure 3

GO enrichment analysis of DEGs. Significantly GO enriched terms in biological process (A), molecular function (B), and cellular component (C). Each dot represented one DEG. Blue and red bars displayed decreasing and increasing z-score. The description of GO terms was listed on the bottom right.

Full size image

In GSEA analysis, the gseGO outcomes showed that DEG expressions were mainly involved in response to a chemical, the immune system process, and response to an external stimulus (Fig. 4A). As for gseKEGG results, DEGs were associated with the TNF, IL-17, and C-type lectin receptor signaling pathway (Fig. 4B). Of these, the signaling pathways of cytokine-mediated, TNF, and IL − 17 were reported to accumulate in the lung epithelial cells within acute inflammation and serious viral infection [23]. Since ALI is characteristic of severe inflammation of lungs resulting from uncontrolled inflammatory immune response, the outcomes implied that these DEGs might play essential roles in the pathogenesis and progression of ALI/ARDS. Given the functional enrichment of these DEGs, we delved further into the results in an integrative meta-analysis.

Fig. 4
figure 4

GSEA analysis of DEGs in gseGO (A) and gseKEGG (B). The size of dot showed count number of enriched genes in each term. The color of dot displayed adjusted p-value.

Full size image

Analysis of PPI Networks of DEGs

We searched for hub genes via NetworkAnalyst which might be useful as biomarkers and therapeutic targets for ALI. The “first-order” PPI network for ALI was difficult to interpret and navigate due to it yielded 3055 nodes and 7107 connection edges. To better present the PPI network, “zero order” network analysis was conducted possessing four subnetworks (at least 5 nodes), including one big subnetwork (subnetwork1 contained 223 nodes and 423 connection edges) and three smaller ones (subnetwork2–4; Table S5). We used Cytoscape software to better visualize the subnetwork1 (Fig. 5). DEGs within the subnetwork1 were ranked by the betweenness, and the top 7 DEGs (Stat1, Syk, Jak3, Rac2, Ripk1, Traf6, and Mapk3) were regarded as hub genes whose node degrees ≥ 10. Among them, Mapk3 was the most highly ranked key gene (degree = 12; betweenness = 3910.56) among the downregulated DEGs. Meanwhile, Stat1 was the centermost upregulated gene (degree = 23; betweenness = 10,125.24) followed by Syk (degree = 21; betweenness = 2883.29).

Fig. 5
figure 5

Module visualization and identification of hub genes. Visualization of the gene co-expression network of the DEGs was generated using Cytoscape software. The color of dot represented combined ES of DEGs, which varied from blue to red. Seven hub DEGs were highlighted by diamond dots while others were represented by circles.

Full size image

Verification of Functional Roles Using an Established Mouse Model

To next confirm the candidate 7 core DEGs (Stat1, Syk, Jak3, Rac2, Ripk1, Traf6, and Mapk3) screened with this network-based approach, we employed an established ALI mice model through intratracheal instillation of LPS and utilized qRT-PCR validation. The results revealed that the relative mRNA levels for Stat1, Syk, Jak3, Rac2, Ripk1, and Traf6 were increased, while Mapk3 was lower expressed in the LPS group than the control (Fig. 6A–G; Table S6). Except those in Traf6 (p-value = 0.1544) and Mapk3 (p-value = 0.2156), LPS vs control comparisons in other genes showed the significant differences (p-value < 0.05). Overall, the qRT-qPCR results were in accord with our integrative meta-analysis, suggesting the critical role of the 7 key DGEs that might play in the pathogenic mechanism of ALI/ARDS. We subsequently tested whether these DEGs had an intrinsic relationship in ALI via analyzing the STRING database. We noticed that Traf6 had the most significant connections with other genes followed by Stat1 (Fig. 6H).

Fig. 6
figure 6

(A–G) Expression levels of seven core genes (Stat1, Syk, Jak3, Rac2, Ripk1, Traf6, and Mapk3) in the lung of established mouse model. GAPDH was used as an internal control gene. Experiments were performed at least in triplicates. ****, ***, **, *, and ns represent p-value < 0.0001, < 0.001, < 0.01, < 0.05, and no difference, respectively. (H) Functional links among core genes determined by STRING database.

Full size image

DISCUSSION

ALI and ARDS are correlated with short-term and long-term syndromes, including physical and cognitive impairments [24, 25], as well as poor clinical conditions outcomes with mortality rates as high as 35–55% [26], whereas the underlying molecular mechanism of ALI/ARDS is far from being understood. In the present study, to reconnoiter the promising gene changes and key mediators that occur in ALI, we built a transcriptional compendium of the genes concerning LPS-induced ALI from six microarray studies that covered 76 tissues in total. Ultimately, 958 overlapping DEGs were confirmed and found to be enriched in several known terms of GO and KEGG pathways related to inflammatory responses and necroptosis. A PPI network-based analysis was used to evaluate the relationship between these genes. This analysis suggested these genes interrelated with four significant subnetworks. We detected 7 hub genes (Stat1, Syk, Jak3, Rac2, Ripk1, Mapk3, and Traf6) with at least 10 node degrees. Among them, Stat1 was the centermost over-expressed gene, followed by Syk. Mapk3 was the most clearly downregulated gene. We verified mRNA levels of these seven genes between LPS-challenged mice and normal mice. These outcomes demonstrated that Stat1, Syk, Jak3, Rac2, Ripk1, and Traf6 were markedly increased in an established aging model, as well as Mapk3 was decreased.

The innate immune response exerts a profound impact on the pathogenesis of ALI [27]. Miscellaneous signal transduction pathways partake in mediating lung inflammatory responses, such as the JAK/STAT, NF-κB, and MAPK signaling pathways [28]. Among them, the JAK/STAT pathway are critical determinants of initiating adaptive immune responses and constraining inflammatory responses eventually [29]. Janus tyrosine kinase 3 (JAK3), one of the JAK family members, is an attractive selective target for the treatment of immune-mediated disorders [30]. Inhibiting JAK3 reduces the hyperproduction of cytokine/chemokine [31]. STAT1 is an important component of inflammation and tumorigenesis. Recent studies have demonstrated that the downregulation of STAT1 might inhibit lung cell apoptosis in ALI [32]. As for spleen tyrosine kinase (SYK), it modulates inflammatory signaling in various immune cells via classical immunoreceptors [33]. LPS results in increasing Syk expression in neutrophils and gamma delta (γδ) T cells directed to inflammatory process of ALI [34].

Neutrophils, dendritic cell (DC)s, and macrophages are associated with multiple immunological processes and tissue injury by initiating inflammation during ALI [35]. As a member of the Rho GTPase subfamily, RAC2 acts as a necessary regulator of neutrophil degranulation. In alveolar macrophages and neutrophils, lung injury in response to immune complex deposition is reliant on Rac2 [36]. MAPK3 is a negative moderator of DC and is necessary for the induction of T cell anergy towards an inflammatory phenotype [37]. Tumor necrosis factor receptor–associated factor 6 (TRAF6) is a decisive adaptor regulating Toll-like receptors (TLRs), whose polyubiquitination leads to mediate pro-inflammatory cytokines production in ALI [38]. Moreover, necroptosis, a cell death form modulated by the RIPK1-RIPK3-MLKL signaling, has been implicated in the close relationship between necroptosis and ALI/ARDS [39, 40]. As cellular signaling molecules, receptor-interacting protein kinases (RIPKs) are critical for homeostatic signaling in many disease processes [41]. The activation of RIPK1/3-dependent necroptosis would result in ALI/ARDS and even death [42].

In summary, as mentioned above, Stat1, Syk, Jak3, Rac2, Ripk1, and Traf6 were the central upregulated genes in our network analysis. The expression of those DEGs was higher in LPS-induced pathological injury in the lung than control counterparts. The enhanced DEGs may be critical factors in attenuating LPS-stimulated ALI/ARDS. These observations coincided with the earlier reports proving our analyses are reliable and practical. On the other hand, some DEGs like Ebi3 and F10 have been described previously in many other illnesses [43,44,45], but their regulatory functions in ALI/ARDS have not been fully known. Additionally, this configuration of network analysis may serve as a novel tool to screen for corresponding potential signatures in an attempt to fill the gaps of knowledge about properties and pathogenesis of ALI. To address the limitations in this study, further research using knockout gene mice for each DEGs is indispensable and in urgent need. It will facilitate us to get better acquainted with its role of LPS in aggravating ALI/ARDS.