Identification and validation of candidate genes dysregulated in alveolar macrophages of acute respiratory distress syndrome

Ethics statement

This study was approved by the ethics committee of Shanghai Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine. Written informed consent was obtained from all patients and volunteers.

Data collection and study design

Two mRNA expression datasets of ARDS (GSE89953 and GSE116560) (Morrell et al., 2019; Morrell et al., 2018) were downloaded from the gene expression omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) database. GSE89953 included 26 alveolar macrophage samples isolated from bronchoalveolar lavage fluid of ARDS patients and paired 26 peripheral blood monocyte samples. As for GSE116560, 68 alveolar macrophage samples were collected from ARDS patients and 30 of 68 were obtained on ARDS onset day 1. The above datasets were detected using GPL6883 platforms and natural scale quantile normalization was performed.

As the initial response to injury, the exudative phase of ARDS is characterized by alveolar macrophage activation inducing the release of pro-inflammatory mediators and chemokines that promote the accumulation of monocytes and neutrophils. Therefore, alveolar macrophages which were derived and differentiated from peripheral blood monocytes were considered as the main contributor to alveolar epithelial and endothelial damage. To clarify the molecular mechanism underlying alveolar macrophage activation and accumulation in the early phase of ARDS, mRNA expression data of GSE89953 was used to screen and identify the DEGs in alveolar macrophages of ARDS. Moreover, we utilized GSE116560 to conduct prognostic analysis of the hub genes obtained from the GSE89953 dataset.

DEGs screening and identification

To investigate the underlying molecular mechanism of macrophage activation and infiltration into pulmonary alveoli in the early stage of ARDS, mRNA expression dataset (GES89953) was downloaded and differential mRNA expression in alveolar macrophage relative to peripheral blood monocyte was acquired using the limma package in R language. P < 0.01 and —log₂fold change (FC)—≥1 was set as threshold values to define DEGs.

Principal component analysis

Principal component analysis (PCA) is a multivariate statistical method that identifies patterns and classifies the factors which influence a certain phenomenon. As a technique widely used in medical field, PCA decomposes a set of correlated variables into a set of uncorrelated variables named as principal components (PCs) and reduces dimensionality by discarding the less important principal components. In this study, we used the Principal Component Analysis tool of Origin software to perform PCA basing on the mRNA expression matrix of the top 300 DEGs across the 26 alveolar macrophage and paired peripheral blood monocyte samples of ARDS patients. Since 52 samples were included, 52 PCs were generated.

Pathway and process enrichment analysis

Metascape (http://metascape.org) is an efficient, free, web-based and user-friendly analysis tool for experimental biologists to comprehensively analyze and interpret OMICs-based studies (Zhou et al., 2019). In this study, Metascape was chosen to conduct process and pathway enrichment analysis for the top 300 DEGs. For this purpose, the Gene Ontology (GO) terms for biological process and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were enriched using the Metascape analysis tool. P-values <0.01, a minimum count of 3 and enrichment factor >1.5 were chosen as cutoff values and remaining significant terms were grouped into clusters based on Kappa-statistical similarities. Specifically, Kappa score of 0.3 was applied as the threshold and terms with a similarity >0.3 were considered as a cluster. The term which was most statistically significant within a cluster was selected to represent this cluster.

Next, Gene Set Enrichment Analysis (GSEA) (Subramanian et al., 2005) was performed basing on expression data of GSE89953 to explore the potential biological processes and pathways involved in the activation and accumulation of macrophages in ARDS. GO and KEGG gene sets were chosen as priori gene sets for functional enrichment analysis.

Weighted gene co-expression network analysis and module identification

Weighted gene co-expression network analysis (WGCNA) was performed using iDEP.91 web application (https://bioinformatics.sdstate.edu/idep/) (Ge, Son & Yao, 2018) and applied to construct gene co-expression networks based on mRNA expression matrix (GSE89953). The processes of WGCNA mainly included several steps. Firstly, a gene-gene similarity network was established and, afterwards, divided into network modules basing on expression similarity of group genes. Thirdly, the correlation between modules and phenotypic traits was established and, finally, the “driver” genes in modules were identified. β value of 14 was chosen as soft threshold on the basis of scale independence and mean connectivity plots.

Protein–protein interaction analysis

A protein–protein interaction network of the top 300 DEGs was constructed through the STRING database (http://string-db.org/) (Szklarczyk et al., 2019) with a confidence score >0.7 and the PPI network was visualized by Cytoscape (Version 3.7.2) software. The module analysis of the PPI network was performed using the Molecular Complex Detection (MCODE) of Cytoscape and functional enrichment analysis of the obtained modules were carried out using Metascape tool. Moreover, hub genes were screened out using cytoHubba tool of Cytoscape calculating by MCC method and the top 10 hub genes were chosen for further analysis.

Prognostic analysis of the hub genes

mRNA expression matrix of GSE116560 was downloaded from GEO database and the 30 alveolar macrophage samples gathered on ARDS onset day 1 were selected for prognostic analysis. Ventilation-free days (VFDs) was a common description for clinical outcomes of ARDS patients. In this study, VFDs > 0 was defined as subjects who were alive and liberated from mechanical ventilation within 28 days and VFDs = 0 was subjects who died or were persistently supported on mechanical ventilation. For these 30 ARDS patients whose alveolar macrophage samples were available on the first day, the VFDs ranged from 0 to 25 and the average value of 12 was set as the cutoff value. The patients whose VFDs exceeded 12 were defined as individuals with better clinical outcomes and those with VFDs < 12 were accepted as poor prognosis individuals. T test was used to investigate the association between transcript levels of the hub genes and the clinical outcomes and P < 0.05 was defined as statistically significant.

Cell culture and treatment

The human monocytic cell line THP-1 was purchased from the American Type Culture Collection (Manassas, VA, USA) and cultured in RPMI-1640 medium with 10% fetal bovine serum, 100 U/ml penicillin and 100 µg/ml streptomycin at 37 °C in an incubator with a humidified atmosphere containing 5% CO₂.

1 × 10⁶ THP-1 cells were seeded into 6-well plate and differentiated with PMA (Sigma-Aldrich, St Louis, MO, USA) for 48 h in a concentration of 50ng/ml. Cell adhesion, spreading and decreased nucleocytoplasmic ratio, as the main features for macrophage differentiation, were visualized by phase-contrast microscopy. Surface markers of CD14 (monocyte marker) and CD68 (macrophage marker) during differentiation were detected by flow cytometry. To measure the response of PMA-differentiated macrophages to LPS (Escherichia coli 055:B5; Sigma-Aldrich, St Louis, MO, USA), cells were exposed to LPS (1 µg /ml) for 12 h after a recovery period (an incubation of 12 h in fresh RPM-1640 removing PMA and 10% FBS) and levels of TNFα in cell supernatants were quantified using Human TNFα ELISA Kit (Sangon Biotech, Shanghai, China) according to the manufacturer’s instructions. After the acquisition of a macrophage-like phenotype, cells were treated with LPS (1 µg /ml) for 0, 1, 3, 6, 12 and 18 h, respectively, and macrophage model of LPS-induced ARDS was successfully established.

Primary alveolar macrophages isolation

BALF samples were collected from six mechanically ventilated patients diagnosed as ARDS according to the Berlin Definition (ARDS Definition Task Force et al., 2012). Samples were collected within 3 days after-diagnosis and, meanwhile, four BALF samples were obtained from healthy volunteers as negative controls. Magnetic-activated cell sorting was performed and negative selection for alveolar macrophages was achieved by incubating cells collected from the BALF samples with antibody-conjugated microbeads specific for the following markers: CD3 (T cells), CD15 (neutrophils), CD19 (B cells), CD235a (red blood cells), CD229f (eosinophils, basophils) and CD326 (epithelial cells).

RNA isolation and qRT-PCR

Total RNA extraction from cell samples was performed using TRIzol reagent (Thermofisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions. Primers were synthesized by Sangon Biotech and the sequences were shown in Table S1. The reverse-transcription was performed using PrimeScript™ RT Master Mix (Takara, Tokyo, Japan). The mRNA expression levels of 6 hub genes (CCL13, FPR3, PLAU, CX3CR1, CXCL16 and CCR2) in cells were detected by qPCR using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) and calculated as described (Yu et al., 2016). GAPDH mRNA was used as the internal reference and the relative mRNA expressions were calculated using the 2^−ΔΔCT method.

Western blot

Cells were lysed using RIPA lysis buffer (Beyotime, Shanghai, China) in the presence of Cocktail protease inhibitor (Abmole, USA) and Cocktail phophatase inhibitor (Abmole, USA). Lysate protein concentrations were determined by BCA Protein Assay Kit (Beyotime, Shanghai, China). Equivalent protein from different samples was separated by 10% SDS-PAGE protein electrophoresis and the separated protein samples were transferred onto the Immobilon PVDF membranes. Then, the membranes were incubated with 5% skim milk to block the non-specific sites. Rabbit polyclonal anti-human FPR3 antibody (orb608040, Biorbyt, UK), rabbit polyclonal anti-human CCR2 antibody (NBP1-48337, Novus Biologicals, USA) and HRP-conjugated secondary antibody (1:2000, Proteintech, IL, USA) were used. GAPDH (1:5000, Proteintech, IL, USA) was simultaneously used as an internal control. Signals were detected and visualized with Immobilon Western Chemiluminescent HRP Substrate (Merck Millipore, Darmstadt, Germany).

Flow cytometry analysis

THP-1-derived macrophages or primary alveolar macrophages collected from ARDS patients were resuspended in 50 µl of PBS. Cell suspension was incubated with 2 µl of APC anti-human CD14 antibody (Biolegend, USA) , FITC anti-human CD68 antibody (Biolegend, USA) and PE anti-human CCR2 antibody (Biolegend, USA) in dark at 4 °C for 30 mins. Flow cytometry analysis was carried out with BD LSRFortessa™ X-20 cell analyzer (BD Biosciences, USA).

Chemotaxis assay

THP-1 cells were differentiated into macrophages using PMA as described above. siRNA-FPR3, siRNA-CCR2 and siRNA-control were purchased from Fubio Biological Technology (Shanghai, China) and transiently transfected into PMA-differentiated macrophages. At 18 h post-transfection with siRNA-FPR3, siRNA-CCR2 or siRNA-control, macrophages were added into the upper chamber of Transwell plate (24-well plates, 8.0 µm pores, Corning, USA) in 100 µl serum-free medium containing LPS (1 µg /ml) and allowed to migrate toward CCL2 (100 ng/ml) (300-39, Peprotech, USA), fMLP (10⁻⁸M) (F3506, Sigma-Aldrich, St Louis, MO, USA) or culture medium in the lower chamber. After a 3 h incubation at 5%CO₂ and 37° C, non-migrated macrophages (upper chamber) and migrated macrophages (lower chamber) were collected and counted.

Statistical analysis

All data were presented as mean ± SD. Statistical analysis was performed with GraphPad software version 6.01. Comparison between 2 groups was carried out using unpaired, 2-tailed, Student’s t test. When there were more than 2 groups, one-way AVOVA was performed, followed by Bonferroni correction when necessary. Differences were considered to be statistically significant when the p-value was less than 0.05.

Screening and identification of differentially expressed genes in alveolar macrophages of ARDS patients

To identify the DEGs in alveolar macrophages of ARDS patients, mRNA expression data (GSE89953) was downloaded from GEO database. Variance stabilization and quantile normalization was performed for the raw microarray data and screening and identification of the DEGs in ARDS alveolar macrophages compared to peripheral blood monocytes was accomplished using the limma package of R language. According to the thresholds set (P < 0.01 and —log2FC— ≥ 1), 385 up-regulated and 325 down-regulated genes were found respectively (Fig. 1A). Subsequently, we selected the top 300 DEGs ranked in ascending order of adjusted p-value for further analysis. Basing on these 300 DEGs, heatmap was generated using the online TBtools software (Chen et al., 2020) and revealed a remarkable disparity of transcriptional pattern between alveolar macrophages and peripheral blood monocytes in ARDS patients (Fig. 1B).

Principle component analysis (PCA)

The PCA results indicated that 3 of the 52 PCs correspond to 81.2% of the overall variance. More specifically, principal component 1 and principal component 2 accounted for 74.6% and 4.1% of overall variance, respectively. As shown in Fig. 1C, a clear segregation of gene expression was revealed between alveolar macrophages and peripheral blood monocytes in ARDS patients.

Pathway and process enrichment analysis

Process and pathway enrichment analysis were performed using Metascape analysis tool to predict the potential biological processes and pathways participating in ARDS onset. As shown in Fig. 2A, the most enriched GO terms for biological process were ‘myeloid leukocyte activation’, ‘cell chemotaxis’, ‘regulation of inflammatory response’, ‘positive regulation of hydrolase activity’ ,‘lymphocyte activation’ , ‘positive regulation of defense response’, ‘leukocyte differentiation’, ‘response to wounding’, ‘organic hydroxy compound metabolic process’, ‘phagocytosis’ and ‘regulation of cytokine production’. Meanwhile, KEGG analysis indicated that the top 300 DEGs significantly enriched in pathways of ‘Phagosome’, ‘Apoptosis’, ‘B cell receptor signaling pathway’, ‘Chemokine signaling pathway’ and ‘PPAR signaling pathway’.

GSEA indicated that up-regulated genes in alveolar macrophages statistically enriched in GO and KEGG gene sets of ‘Regulation of Non Canonical Wnt Signaling Pathway’, ‘Monocyte Chemotaxis’, ‘Adenylate Cyclase inhibiting G protein Coupled Receptor Signaling Pathway’, ‘Macrophage Migration’, ‘ERK1 and ERK2 Cascade’, ‘Negative Regulation of Apoptotic Signaling Pathway’, ‘Cytokine-Cytokine Receptor Interaction’, ‘Regulation of Autophagy’, ‘PPAR Signaling Pathway’, ‘Complement and Coagulation Cascades’ and ‘JAK-STAT Signaling Pathway’, which was consistent with the results obtained from Metascape tool.

Gene co-expression networks identification and enrichment analysis of modules

Lowly expressed genes were filtered out and 819 remaining genes were used in WGCNA analysis. According to the criterion of scale free topology fix index reaching 0.9 and an appropriate mean connectivity value, the soft threshold value was set as 14. Finally, eight co-expression gene modules were identified by the dynamic tree cut method and each module was marked by a color. The turquoise was the largest module with 422 genes and the pink was the smallest with 23 genes. The clustering dendrogram of genes was displayed in Fig. 2B.

Biological process and pathway enrichment analysis of the eight modules were performed using Metascape tool as previously described and a substantial proportion of the resulting terms shared in common with those obtained from previous analysis. As shown in Fig. 2C, genes of module 1, 2, 6 and 8 mainly involved in immune response, leukocyte activation and response to external stimulus. Besides, module 3 probably participated in cholesterol mechanism. As for module 4, genes significantly enriched in neutrophil chemotaxis and migration. Meanwhile, KEGG analysis indicated a significant enrichment in pathways of ‘complement and coagulation cascades’, ‘cytokine-cytokine receptor interaction’, ‘Cell adhesion molecules’, ‘MAPK signaling pathway’ and ‘apoptosis’ (Fig. 2D).

Construction of PPI network and identification of hub genes

To further explore the promising molecular complexes responsible for the aberrant activation of alveolar macrophages in ARDS, the PPI network was constructed basing on the top 300 DEGs using the STRING database and visualized by Cytoscape software (Fig. 3A). Moreover, MCODE tool was applied to identify functional modules and the top 5 ranked clusters were showed in Fig. 3B.

Afterwards, we performed functional enrichment analysis for these promising functional clusters using Metascape. For cluster 1, the most enriched terms were ‘adenylate cyclase-modulating G protein-coupled receptor signaling pathway’ , ‘G alpha (i) signaling events’, ‘cellular calcium ion homeostasis’, ‘response to tumor necrosis factor’ and ‘Neuroactive ligand–receptor interaction’ (Fig. 3C). A potential involvement of cluster 2 was revealed in oncogenic MAPK signaling, blood coagulation, apoptosis and cellular response to organonitrogen compound. Meanwhile, it was predicted that genes within cluster 3 implicated in ‘protein processing’, ‘cellular component disassembly’ and ‘myeloid leukocyte activation’. Besides, enrichment analysis indicated a possible involvement in endocytosis process for cluster 4 and in neutrophil activation and cell adhesion process for cluster 5, respectively.

Basing on the constructed PPI network, the hub genes which may play a key role in ARDS pathogenesis were identified by MCC method and the top 10 hub genes were chosen for subsequent research and showed in Fig. 3D. Noticeably, the 10 hub genes were coincident with the genes within cluster 1 to a large extent indicating a crucial role of cluster 1 in ARDS onset and progress through Gi-protein-coupled receptor signaling and phosphatidylinositol-calcium second messenger activation.

Prognostic value of the hub genes

To explore the predictive capability of the 10 hub genes, we used t test to analyze the differences of transcript level between better (VFDs > 12) and poor prognosis group (VFDs < 12). Analysis of mRNA expression profiles from GSE116560 indicated that none of the 10 hub genes showed statistically significant difference of mRNA expression between better and poor prognosis group. However, CCL13 (p = 0.0558), CCR2 (p = 0.1578) and FPR3 (p = 0.1611) showed a mild trend of higher transcript level in poor prognosis subjects (Fig.3E). Meanwhile, a slight but not statistically significant over-expression of PLAU was observed in better outcome individuals (p = 0.2631) (Fig.3E). These findings implied that CCL13, CCR2, FPR3 and PLAU may play a crucial role in the onset and progression of ARDS, though no statistically significant result was concluded.

Establishment of cell model of LPS-induced ARDS

For the purpose of differentiation into a macrophage phenotype, THP-1 cells were incubated with PMA (50 ng/ml) for 2 days and observed by phase-contrast microscopy. Compared to untreated cells, PMA treatment remarkably enhanced the adherence of THP-1 cells and reduced the nucleocytoplasmic ratio and these changes were viewed as main features for macrophage differentiation (Fig. S1A). Flow cytometry analysis revealed that, compared to THP-1 cells with higher expression of surface marker CD14, THP-1-derived macrophages differentiated by PMA expressed higher levels of both CD14 and CD68 on surface (Fig. S1B). Mature THP-1-derived macrophages secrete TNFα in response to pro-inflammatory stimuli such as LPS (Park et al., 2007). To assess the pro-inflammatory response of THP-1-derived macrophages to LPS, cells were treated with LPS (1 µg/ml) for 12 h and the secretion of TNFα was evaluated by ELISA. Data indicated that fully-differentiated macrophages pre-treated with PMA secreted higher levels of TNFα in response to LPS (P < 0.05 versus LPS-stimulated THP-1 cells), whereas non-differentiated THP-1 cells showed negligible TNFα secretion in the presence of LPS stimulation (Fig. S1C). These data verified that THP-1 cells differentiated with PMA (50 ng/ml) can be used as a model of macrophages. After optimized THP-1 differentiation into macrophages, cells were challenged with LPS (1 µg/ml) for a varying time of 0, 1, 3, 6, 12 and 18 h and macrophage model of LPS-induced ARDS was successfully established.

Up-regulated expression of FPR3 and CCR2 in macrophage model of LPS-induced ARDS

Next, we sought to validate the expression trend of the hub genes obtained from mRNA microarray assay. Using GeneCards database (http://www.genecards.org/), we selected 6 hub genes (CCL13, FPR3, PLAU, CX3CR1, CXCL16 and CCR2) for further study mainly due to their potential involvement in chemotaxis-related signaling pathways. Afterwards, qRT-PCR was performed and the expression levels of CCL13, FPR3, PLAU, CX3CR1, CXCL16 and CCR2 in macrophage model of LPS-induced ARDS were detected. As shown in Fig. 4A, LPS stimulation induced significantly enhanced FPR3 expression in a time-dependent manner, which was coincident with the previous microarray data. Meanwhile, we also observed that CCR2 mRNA expression was dramatically elevated after LPS stimulation, though this trend was contradictory with that obtained from microarray assay. Afterwards, Western blot analysis was performed to examine the protein levels of FPR3 and CCR2. A substantial up-regulation of CCR2 protein on macrophages after LPS addition was detected. Likewise, LPS stimuli dramatically enhanced FPR3 protein level in a time-dependent manner (Fig. 4B). Moreover, flow cytometry indicated a significantly elevated expression of CCR2 on surface of PMA-differentiated macrophages after LPS stimulation (Fig. 4C). However, FPR3 was hardly detected on surface of THP-1-derived macrophages no matter LPS was administered or not (Fig. S2).

Elevated expression of FPR3 and CCR2 in primary alveolar macrophages of ARDS patients

To further confirm the expression levels of FPR3 and CCR2 in primary alveolar macrophages of ARDS patients, BALF samples were collected from 6 ARDS patients and 4 healthy volunteers and primary alveolar macrophages were isolated using magnetic-activated cell sorting. qRT-PCR indicated that the expression of CCR2 was higher in primary AMs of ARDS patients than in those of healthy volunteers (P < 0.05 versus healthy volunteers). Likewise, an elevated FPR3 expression was detected in primary AMs from ARDS patients compared with healthy individuals (P < 0.05 versus healthy volunteers) (Fig. 4D). Moreover, we also observed a significant increase of FPR3 and CCR2 protein levels in primary AMs of ARDS patients (Fig. 4E). Besides, flow cytometry analysis revealed that expression of surface marker CCR2 was significantly elevated in primary AMs of ARDS patients (Fig. 4F) and FPR3 was barely detected (Fig. S3), which was in line with our previous findings in macrophage model of LPS-induced ARDS.

Collectively, these data confirmed the ectopic-expression of FPR3 and CCR2 in both macrophage model of LPS-induced ARDS and primary alveolar macrophages collected from ARDS patients. These observations suggested a possible participation of FPR3 and CCR2 in the initial inflammation response of alveolar macrophages in ARDS.

FPR3 and CCR2 down-regulation attenuated LPS-induced macrophage chemotaxis toward CCL2

Using in vitro transwell assay described above (Fig. 5A), we examined the ability of macrophage chemotaxis toward CCL2. LPS stimulation induced significantly enhanced macrophage migration toward CCL2 than toward culture medium (P < 0.05 versus PBS-treated macrophages) (Fig. 5B), which may correlate with LPS-induced over-expression of CCR2 on macrophages (Figs. 4A–4C). This was confirmed by an inhibitory effect of siRNA-induced CCR2 silence on LPS-triggered macrophage chemotaxis toward CCL2 (P <  0.05 versus LPS-stimulated macrophages) (Fig. 5B). Therefore, our findings suggested that CCL2/CCR2 axis may take part in the monocyte/macrophage activation and migration into alveolar spaces in early ARDS. Further, we assessed the biological effect of FPR3 on LPS-induced macrophage migration toward chemokine CCL2. As shown in Fig. 5B, LPS-induced increase in macrophage chemotaxis was strongly attenuated by siRNA-FPR3 transfection (P <  0.05 versus LPS-stimulated macrophages). However, neither FPR3 nor CCR2 silence affected the LPS-induced macrophage migration toward chemoattractant fMLP which is also known as a common ligand for FPR2 and FPR1 inducing chemotaxis (Fig. 5C).

Considering the possible link between FPR3 and CCR2 shown in Fig. 3B, pearson correlation analysis was carried out with mRNA expression profiles of GSE89953 and GSE116560. Interestingly, we noticed a strong association of CCR2 with FPR3 mRNA expression in alveolar macrophages of ARDS patients (Fig. 5D). Moreover, we conducted FPR3 silence by transient siRNA transfection and observed that LPS-induced over-expression of CCR2 on macrophages was significantly inhibited (Fig. 5E), while FPR3 protein level was not obviously affected by siRNA-CCR2 transfection (Fig. S4).