Identification of an Alveolar Macrophage-Related Core Gene Set in Acute Respiratory Distress Syndrome

Introduction

Sepsis is a syndrome in which the body’s immune system overreacts to infection. There may be several life-threatening complications. One of them is acute respiratory distress syndrome (ARDS).^1,2 ARDS is an acute diffuse lung injury with high incidence. It has been shown that the pathophysiological basis of ARDS includes an excessive and protracted systemic inflammation,³ but the exact mechanisms remain unknown. The main risk factors that may cause ARDS include drinking, smoking, air pollution, hypoproteinemia, and diabetes mellitus.^4–6 Currently, urgent treatment strategies mainly include low tidal mechanical ventilation, but the outcome is very poor.³ Although some progress has been made in treatment options in recent years, the mortality rate remains high.⁷

Previous studies have identified some genes used as biomarkers for ARDS. The high expression of interleukin (IL)-33 can elevate the expression of matrix metallopeptidases (MMP) 2 and 9 in patients with acute lung injury.⁸ SLC2A6 may be a critical biomarker for predicting survival of sepsis patients.⁹ V-set immunoglobulin-domain-containing 4 (VSIG4) may be involved in lung injury by inducing phagocytosis.¹⁰ Activation of ERK1/2 by NGF1B depends on MAPKKK c-Raf, eventually inducing inflammation. The NGF1B in the cytoplasm may be involved in the occurrence and development of ARDS.¹¹ Down-regulation of the long non-coding RNA GAS5 decreases the expression of angiotensin-converting enzyme 2 (ACE 2), leading to increased levels of microRNA (miR) miR-200c-3p, which promotes the progression of ARDS.¹² MYC and STAT3 may be the key regulatory genes in the underlying dysfunction of sepsis-induced ARDS.¹³ In addition, a previous study showed that gene set variation index can be used the potential diagnostic tool for sepsis.¹⁴ And gene set variation index may be considered as the biomarker of bacterial and fungal sepsis.¹⁵ However, we found that a majority of previous studies mostly focused on single genes or molecules for ARDS. There are only a few studies involving in a set of genes related to ARDS.

Several studies confirmed that alveolar macrophage transcriptional programs are associated with ARDS.^16–18 We hypothesized that these alveolar macrophage-related ARDS-specific transcriptional programs may be also reflected in blood cells. To investigate our hypothesis in the present study, we identified genes upregulated both in alveolar macrophages and blood cells and constructed a corresponding protein–protein interaction network associated with host response to ADRS. We found that eight genes formed a core gene set, and that the core gene set variation score may be used as a circulating biomarker for ARDS.

Materials and Methods

Data Collection and Processing

The data was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/): GSE116560,¹⁶ GSE89953,¹⁹ GSE76293²⁰ and GSE32707.²¹ GSE116560 and GSE89953 were based on GPL6883. GSE116560 was obtained from alveolar macrophages of 68 ARDS samples, while GSE89953 was obtained from peripheral blood monocytes of 26 ARDS samples. The two datasets were combined and then batch effects were removed using the sva routine in R.²² From GSE76293, based on GPL570, the gene profiles of neutrophils from 12 ARDS samples and 12 healthy controls were extracted for subsequent analysis. From GSE32707, based on GPL10558, the gene profiles of whole blood cells were used to perform the subsequent analysis, from 18 ARDS samples and 34 healthy controls. The data were normalized using the “normalizeBetweenArrays” function of the limma routine²³ in R. The probe set was converted into gene symbols, according to the annotation information of the three platforms. The workflow of the study is shown in Figure 1.

Figure 1 Workflow chart of the study.

Principal Component Analysis (PCA) and Screening Differentially Expressed Genes (DEGs)

The ggbiplot routine²⁴ in R was applied to perform the PCA.²⁵ In GPL6883 datasets (GSE116560 and GSE89953), the DEGs in alveolar macrophages were screened in comparison with peripheral blood monocytes. These DEGs may reflect the response of host monocyte-macrophages to ARDS.¹⁸ In GSE76293, DEGs in ARDS neutrophils were screened compared to healthy controls, since neutrophils are the most frequent nucleated cell in whole blood. In GSE32707, we screened the DEGs in whole blood between ARDS patients and healthy controls. A total of the DEGs were screened using the limma routine in R. The cut-off criteria are a false discovery rate (FDR)-adjusted P < 0.05 and |log2 (fold-change)| > 1. In the study, the common upregulated genes were defined as circulating ARDS alveolar macrophage-related genes (cARDSAMGs).

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway Enrichment Analysis

To explore the potential biological functions of the cARDSAMGs, GO and KEGG pathway enrichment analysis were performed using the clusterProfiler routine in R.²⁶ GO and KEGG networks were clustered using the plug-in ClueGO²⁶ in Cytoscape software.²⁷ Differences associated with an FDR-adjusted P < 0.05 was considered statistically significant.

Analysis of Core Genes

Using STRING database²⁸ (https://string-db.org/) and KEGG pathway enrichment analysis, a protein–protein interactions (PPIs) network was constructed for the cARDSAMGs. If genes were involved in ARDS-related pathways of interest, we defined them as core genes in the network. The GOsemsim routine inR²⁹ was applied to analyze the semantic similarities among GO terms for the cARDSAMGs, applying the following formula: Pi = -lg (p) * |log2 (fold-change)|, weight (W) value = semantic similarities * mean value of Pi.

Gene Set Variation Analysis (GSVA) and Receiver Operating Characteristic (ROC) Curve Analysis

The GSVA routine inR³⁰ was applied to calculate the core gene set variation analysis (CGSVA) score for individual samples to assess the concerted functional behavior of the core gene set. To further evaluate the potential diagnostic value of the CGSVA score for ARDS, an ROC curve analysis was performed using the pROC routine in R.³¹

Results

Multiple Genes Identified as cARDSAMGs

PCA showed that global gene expression patterns could distinguish ARDS from controls in GSE116560, GSE89953 and GSE76293, but not in GSE32707 (Figure S1). There were 8522 DEGs in GSE116560 and GSE89953, 6064 in GSE76293, and 7713 in GSE32707 (Figure 2A). A set of 60 DEGs was overexpressed in all the datasets, defined as cARDSAMGs (Figure 2B). Moreover, we obtained the p.value and rank of 60 cARDSAMGs in above three datasets and the ranking based on the p.value (Table S1). PCA showed that the expression patterns of the cARDSAMGs could not distinguish ARDS from healthy controls (Figure S2). However, we found that cARDSAMGs can better distinguish ARDS from control compared with the global gene expression patterns (Figures S1 and S2).

Figure 2 Multiple genes identified as cARDSAMGs. (A) Manhattan plot for DEGs. The three most significantly overexpressed genes are marked in yellow and labeled with their names. (B) Venn diagram to identify upregulated DEGs common to datasets from the three platforms. Abbreviation: DEGs, differentially expressed genes.

Biological Significance of cARDSAMGs

GO enrichment analysis showed that changes in expression involved mainly genes related to metabolism- and inflammation-related biological processes (Figure 3A). KEGG result showed that DEGs were mainly involved in pathways related to ARDS, such as Notch the tumor necrosis factor (TNF) and the tumor necrosis factor (TNF) signaling pathways, as well as complement and coagulation cascades (Figure 3B). In addition, we performed the enrichment analysis using the 60 top ranked up-regulated genes in each of the datasets. As shown in Figure S3, we found that the biological processes (BPs) and KEGG pathways were different in each of datasets. But most of them were related to immunity, inflammation and phagocytosis. Furthermore, we randomly selected 60 genes in the three datasets and thereby calculated the enriched GO term. We found that the enriched GO terms mostly were not significantly related with inflammation, immunity and phagocytosis, which further evidenced that the gene set-cARDAMGs was not simply a randomly overlapping set (Figure S4). The ClueGo analysis indicated that DEGs were clustered mainly in biological processes, such as bone resorption, collagen catabolic process, regulation of cartilage development, positive regulation of ATP biosynthetic process, negative regulation of IL-2 production and positive regulation of animal organ morphogenesis (Figure 3C). No significant clustering of KEGG pathways was observed, based on a threshold of FDR-adjusted P < 0.05.

Figure 3 Biological significance of cARDSAMGs. (A) The significantly enriched biological processes, (B) KEGG pathways and (C) GO network in circulating ARDS involving alveolar macrophage-related genes in the three datasets. Different colors of nodes represent different functional groups. Abbreviations: ARDS, acute respiratory distress syndrome; KEGG, Enrichment in Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Core Genes in the PPIs Network of ARDS

We constructed the PPIs network of cARDSAMGs based on the STRING database. There were 25 nodes and 189 interaction pairs (Figure 4A). Previous studies showed that some of the biological processes and pathways were related to ARDS, including immunity, inflammation and phagocytosis.^32,33 Therefore, we selected the following ARDS-related pathways as pathways of interest: Notch, TNF and IL-17 signaling pathways, as well as complement and coagulation cascades. Subsequently, we identified eight cARDSAMGs as core genes involved in these pathways. The core genes included PTCRA, JAG1, C1QB, ADAM17, C1QA, MMP9, VSIG4 and TNFAIP3. In the complement and coagulation cascades, antigen-antibody complexes act on C1QA and C1QB. They further act on VSIG4, thereby inducing phagocytosis.^34,35 In the TNF signaling pathway, TNF acts on TNF receptor 1 (TNFR1), ultimately inducing tissue remodeling, autoimmune pathology, neutrophil recruitment and immunity to extracellular pathogens.^36,37 In the Notch signaling pathway, ADAM17 and JAG1 act on the cell surface receptor Notch, which in turn affects PTCRA. In the TNF signaling pathway, TNF binds to TNFR1, which signals to surface receptors, intracellular signaling and remodeling of extracellular matrix.^38,39 In the IL-17 signaling pathway, the IL-17 family signals activate downstream pathways that include NF-kappaB, MAPKs and C/EBPs, which further induces the expression of chemokines, cytokines and antimicrobial peptides^40,41 (Figure 4B). The semantic similarities among the GO terms of the eight core genes are ranked in Figure 4C, and the W value of core genes, which may suggest the importance of molecules in this dysfunction, are displayed in Figure 4D.

Figure 4 Core genes in the PPIs network of ARDS. (A) PPIs network for alveolar macrophage-related genes in circulating acute respiratory distress syndrome. Ellipses represent genes; diamonds, pathways; and font size, relative degree. (B) Core genes and pathways of interest. (C) Semantic similarities among GO terms for the eight core genes. (D) Weight value of the eight core genes in the three datasets. Abbreviations: PPIs, protein-protein interactions; GO, Gene Ontology.

CGSVA Score May Be Used as a Circulating Biomarker for ARDS

As shown in Figure S5, the C1QA, C1QB, MMP9, PTCRA and VSIG4 may be serve as circulating diagnostic markers for ARDS. Their areas under the ROC curves (AUCs) were higher than 0.7 in the three datasets. In addition, the CGSVA score was higher in the ARDS group than in the control group in all datasets (Figure 5A). The expression of the core genes is shown in Figure 5A. AUC for CGSVA score was higher than 0.7 in all three datasets, which suggests that the score may be useful as a diagnostic marker for ARDS (Figure 5B). The score in macrophages from the alveolar lavage fluid (GSE116560 and GSE89953) and whole blood (GSE76293 and GSE32707) was higher in ARDS individuals than in healthy controls (Figure 5C).

Figure 5 CGSVA score may serve as a circulating biomarker for ARDS. (A) a, GSVA scores of the core gene set in the three datasets. Expression of core genes in (b) GSE76293, (c) GSE116560 & GSE89953 and (d) GSE32707. (B) ROC curves of the gene set GSVA score in the three datasets. (C) The GSVA scores in alveolar whole blood (GSE76293 and GSE32707) and lavage fluid (GSE116560 and GSE89953) were higher for patients with ARDS than for controls. Abbreviations: GSVA, gene set variation analysis; ROC, receiver operating characteristic; ARDS, acute respiratory distress syndrome.

Discussion

In this study, the differences in gene expression were explored between whole blood and alveolar lavage in ARDS samples and healthy controls. We selected four gene expression datasets from three platforms, where we found a common set of 60 DEGs upregulated in ARDS patients compared to healthy controls. Some among them have been previously informed to be related to ARDS. For example, immunohistochemical results showed that CD163 is overexpressed in ARDS patients.⁴² The highly expressed IL-33 can drive expression of MMP9 and MMP2 by activating STAT3. The process mainly happens in alveolar macrophages for acute lung injury patients induced by lipopolysaccharide.⁸ Mutations in the ALOX5 gene appear to be most likely cause of inter-individual differences in ARDS progression.⁴³

Subsequently, we performed a function enrichment analysis for cARDSAMGs. The results showed that these cARDSAMGs were involved mainly in biological processes and pathways related to immunity, inflammation, phagocytosis and nucleic acid turnover. Some of these pathways have been related with ARDS in previous studies. Production of IL-1 and low concentrations of anti-inflammatory cytokines, such as IL-10 and IL-1 receptor antagonist, in bronchoalveolar lavage fluid of early ARDS patients have been closely related to poor prognosis.⁴⁴ Numerous studies have linked inflammation and phagocytosis to ARDS.^32,33,45,46

Moreover, we constructed a protein–protein interaction network to analyze the relationship among genes and pathways. We focused on TNF, IL-17 and Notch signaling pathways, as well as complement and coagulation cascades. In complement and coagulation cascades, complement activation may induce the occurrence of phagocytosis,⁴⁷ which plays an essential role in ARDS.⁴⁸ TNF acts as an inflammatory factor, up-regulating the expression of surface receptors, promoting intracellular signaling and remodeling of the extracellular matrix.⁴⁹ Some biomarkers are closely related to ARDS, including surfactant-related proteins and cytokines (IL-1, −2, −6, −8, −10 and −15, and TNF-α), as well as neutrophil activation markers (MMP9, leukotriene B4 and ferritin).³² In the IL-17 signaling pathway, IL-17A and IL-17B act through surface receptors IL-17-17RA and IL-17RC to induce immunity and recruit neutrophils.^49,50 Neutrophil-derived mediators induce epithelial cell death by oxidation of soluble Fas ligand.^33,45

Eight genes corresponding to the selected pathways were identified as core genes: PTCRA, JAG1, C1QB, ADAM17, C1QA, MMP9, VSIG4 and TNFAIP3. VSIG4 is the main part of phagocytic system: it rapidly removes C3-modified particles.¹⁰ MMP9 is a secreted zinc metallopeptidase.⁵¹ MMP9 mediates protein degradation in the extracellular matrix of alveolar epithelial cells, mainly in intercellular junction proteins and proteins that anchor cells to the basement membrane.⁵² In addition, MMP9 is enriched mainly in macrophages and neutrophils,^53,54 and involved in the occurrence and development of emphysema and asthma.^55,56 TNFAIP3 was identified as a gene whose expression is rapidly induced by TNF. The protein encoded by this gene has been shown to inhibit nuclear factor (NF)-kappa B activation as well as TNF-mediated apoptosis.⁵⁷ PTCRA encodes a single channel type 1 membrane protein in mature T cells. This protein forms a T-cell pre-receptor complex together with CD3 and T cell receptor beta chain (TCRB) complexes, thereby regulating the early development of T cells.^58–60 C1QB may be associated with the immune response after injury.⁶¹

In addition, CGSVA score may be used as a biomarker for ARDS patients. Furthermore, the CGSVA scores associated with ARDS both in whole blood and alveolar macrophages were significantly higher in ARDS patients than in healthy controls. This may suggest that these alveolar macrophage-related ARDS-specific transcriptional programs are also reflected in blood cells.

Our study presents several limitations. Firstly, our predictions are based on bioinformatic analyses, and therefore further experimental verification is needed. Secondly, most of alveolar macrophage genes have a dynamic pattern, but we only select genes that are significantly up-regulated in all time point without the consider of time changes. Thirdly, the whole transcriptome sequencing data of macrophages in whole blood and alveolar lavage fluid are difficult to find, currently. Whether these 8 cARDSAMGs may be used as diagnostic markers needs to verify in a larger independent clinical sample. Fourthly, it is shown that the gene set has a diagnostic performance in each of datasets in this study. Therefore, we believe that this gene set is specific for ARDS to some extent. But whether it is specific enough to ARDS, any inflammation or infection, still needs to be further clarified. Finally, we need to further explore the upstream regulators of core genes to reveal in more detail the imbalance network in ARDS.

Conclusions

We identified a core gene set (PTCRA, JAG1, C1QB, ADAM17, C1QA, MMP9, VSIG4 and TNFAIP3) that may help predict ARDS. The ARDS alveolar macrophage-related CGSVA score may be helpful for the diagnosis of ARDS.