COPD gene expression data collected from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE38974, GSE76925, and GSE106986 datasets were downloaded. The GSE106986 dataset (14 COPD samples and 5 control samples) was derived from GPL13497 platform of the Agilent-026652 Whole Human Genome Microarray 4x44K v2. The GSE76925 dataset (111 COPD samples and 40 control samples) was derived from the GPL10558 platform of Illumina HumanHT-12 V4.0 expression beadchip. GSE38974 dataset (23 COPD samples and 9 control samples) was derived from the GPL4133 platform of the Agilent-014850 Whole Human Genome Microarray 4x44K G4112F. All samples from the 3 different databases were obtained from human lung tissue.

Each dataset was background corrected and normalized using the “limma” package and converted to gene symbols referencing the probe annotation files of the probe names in each dataset. A metadata cohort comprising of merged GSE38974 and GSE76925 cohorts was created for further analysis. Batch variability between platforms was eliminated using the combat function of the “SVA” package (18). The GSE106986 dataset was used as the validation cohort. Similarly, 134 COPD and 49 normal samples were analyzed for DEGs using the R package “limma” with DEGs selected based on thresholds of P < 0.05 and |log2 FC| > 1.

The Metascape (http://metascape.org) database allowed for further understanding of the DEG biological significance using Gene Ontology enrichment analysis. The Molecular Signatures Database (MSigDB) Hallmark Gene Sets and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway were used for pathway enrichment analysis. Meanwhile, the most significant terms which arose in statistical analysis were selected for visualization.

Two machine-learning algorithms, the support vector machine recursive feature elimination (SVM-RFE) and least absolute shrinkage and selection operator (LASSO) were used in this study to screen for significant prognostic variables. SVM-RFE represents a widely used supervised machine-learning protocol for classification and regression that is applied using the “e1071” package. The SVM-RFE algorithm was used to identify genes with higher discriminative power (19). LASSO was performed using the “glmnet” package in R and represents a regression analysis algorithm that applies regularization for variable selection. Using LASSO, we were able to identify genes significantly associated with COPD and normal samples. We applied two algorithms in the metadata cohort and utilized the GSE106986 dataset to analyze overlapping genes from the two algorithms to further validate the expression levels of candidate diagnostic biomarkers.

CIBERSORT is a deconvolution algorithm that quantifies immune cell infiltration (22 various cell types) in COPD gene expression profiles (20). The “corrplot” R package was used to carry out visualization and correlation analysis of 22 types of infiltrating immune cells. Differences in infiltrating immune cells between COPD and healthy samples were visualized using boxplots drawn with the “ggplot2” package in R.

Pearson correlation analysis allowed for in-depth scrutiny of relationships between diagnostic biological markers and infiltrating immune cells. “ggplot2” of the R package was used to visualize the results of the analysis.

Six-week-old C57BL/6 mice was used to establish the COPD model. The details procedure was as described by He et al. (21). Cigarette smoke extract (CSE) was prepared by dissolving the smoke of non-filter Furong cigarette in to PBS (21). The histological changes and immune infiltration of lung tissues were then detected (21, 22). Bronchoalveolar lavage fluid (BALF) was collected. Lung tissue was cut into small pieces and digested into single cell suspension by Collagenase A (4 mg/ml in RPMI1640) for further detection.

After treatment, lung tissues were collected. Paraffin embeded lung tissues were cut into slices for histological staining. HE staining was performed according to the manufacturer’s instructions (Solarbio). The mean alveolar septal thickness (MAST), mean linear intercept (MLI), and destructive index (DI) were calculated. Immunohistochemistry staining was conducted using primary antibodies for SLC27A3 (Proteintech), STAU1 (Affinity), CD57 (Affinity), and CD8a (Affinity). Immunofluorescence staining was performed using primary antibodies for CD19 (Affinity) and CD27 (Santa), and secondary antibodies including Cy3-labeled goat anti-rabbit IgG (Invitrogen) and FITC-labeled goat anti-mouse IgG (Abcam).

BEAS-2B cells were cultured in DMEM cultural medium (Servivebio) with 10% fetal bovine serum (EVERY GREEN) in an incubator (37°C and 5% CO₂). COPD cell model was induced by 5% CSE for 24 h (23). Small interfering RNAs targeting SLC27A3 (si-SLS27A3) and STAU1 (si-STAU1), as well was negative control RNA (si-NC) were transfected into BEAS-2B cells by Lipo 3000 (Invitrogen), respectively. CCK-8 assay kit (Wanleibio) was used to detect viability of BEAS-2B cells according to the manufacturer’s instructions. The optical density was measured by a microplate reader (BIOTEK).

Total protein was extracted from lung tissues and BEAS-2B cells by protein extraction kit (Wanleibio). The concentration of protein was detected by BCA kit (Wanleibio). The protein was separated by SDS-PAGE and transferred to PVDF membrane (Millipore). The membrane was blocked by non-fat milk and then incubated with primary antibodies for SLC27A3 (Proteintech), STAU1 (Affinity), ACPL2 (Bioss), RABL4 (Biorbyt), and β-actin (Wanleibio), respectively. After 12-h of incubation at 4°C, the membrane was washed, incubated by the secondary antibody, and visualized by CEL (Wanleibio). β-actin serves as an internal control.

Total RNA was extracted from lung tissues and BEAS-2B cells by TRIpure (BioTeke). The RNA was reverse transcribed into cDNA using BeyoRT™ II M-MLV RNase H- (Beyotime). Real-time PCR assay was conducted on Exicycler 96 system (BIONEER) using SYBR Green (Solarbio). The results was calculated using 2^-△△CT method. The primer sequences are as follows. Mus musculus ACPL2: F, AATCGCTTCTTGGTGCTG; R, CTACGCTTGGAATGTTGC. Mus musculus RABL4: F, GAAATGCTGGATAAGTTGTG; R, GAGGGAGATGCCTGAAGT. Mus musculus SLC27A3: F, GCATTGTGGGCTGCTTGG; R, GGGCTGGTTGACGAGGTAT. Mus musculus STAU1 F, GTAAAGAAACCAGGAGACG; R, CTGCTGATGGCTAAGATAA.

Mus musculus β-actin: F, CTGTGCCCATCTACGAGGGCTAT; R, TTTGATGTCACGCACGATTTCC. Homo sapiens ACPL2: F, ATGGAGCACTTCAAGGTAA; R, AGCAGAGTAGAGGGCAAA. Homo sapiens RABL4: F, GGAATGGATTTGGTGGTGAA; R, GAGATGCCTGGAGCCTGTGA. Homo sapiens SLC27A3: F, GTTCGGATGGCAAATGAGG; R, TGTACCGGGCAGTTGTGAG. Homo sapiens STAU1 F, ATCCGATTAGCCGACTGG; R, ACTTGAGTGCGGGTTTGG. Homo sapiens β-actin: F, GGCACCCAGCACAATGAA; R, TAGAAGCATTTGCGGTGG.

The concentration of IL-6, IL-1β, and TNFα was measured by ELISA. ELISA kits (Wanleibio) for IL-6, IL-1β, and TNFα were used in this study and conducted according to the manufacturer’s instructions. The optical density was measured by a microplate reader (BIOTEK).

All statistical analyses were performed in R (version 3.6.1). The nonparametric Kruskal-Wallis test was used to perform intergroup comparisons of continuous variables. The degree of efficacy of each diagnostic biomarker was assessed using receiver operating characteristic (ROC) curves. Pearson correlation was used to analyze the relationship between the expression of diagnostic biomarkers and infiltrating immune cells. Comparison between two groups or among multiple groups were analyzed by student t-test or one-way ANOVA, respectively. Statistical significance was identified based on P < 0.05. Furthermore, bioinformatics analysis almost run for two months.

Differential expression analysis between 134 COPD samples and 49 normal samples in the metadata (GSE38974 and GSE76925) cohort was carried out utilizing the “limma” R package. Of the 80 identified DEGs, 4 genes were significantly upregulated and 76 genes were significantly downregulated (Figure 1A). The heatmap of DEGs is depicted in Figure 1B.

GO and pathway analyses were performed to identify the biological functions of DEG using the Metascape database. DEGs were primarily involved with telomerase regulation, cellular amino acid metabolic processes, multicellular organic homeostasis, MHC class I-mediated peptide antigen presentation and antigen processing in addition to various other biological processes (Figure 2A). Moreover, these DEGs were also mainly enriched in the S phase, TriC/CCT association with target proteins during biosynthesis and nucleotide metabolism pathways. Figure 2B shows the relationships between the enriched terms.

Candidate diagnostic biomarkers were screened by two different algorithms. We utilized the LASSO logistic regression algorithm to identify and COPD related 5 feature variables from DEGs (Figure 3A). The SVM-RFE algorithm was used to identify a subset of 42 features in the determined DEGs (Figure 3B). 4 overlapping diagnostic related genes (ACPL2, RABL4, SLC27A3, and STAU1) of these two algorithms were finally selected (Figure 3C). The GSE106986 dataset was used to validate the accuracy of this method as well as the expression levels of the four candidate diagnostic biomarkers. ACPL2 and RABL4 expressions levels were not significantly different between COPD and normal samples (Figures 4A, B). However, the expression levels of SLC27A3 and STAU1 in COPD samples were notably raised in COPD samples in contrast to controls (Figures 4C, D; all P < 0.05). To further test the diagnostic efficacy of SLC27A3 and STAU1, we validated them using metadata and the GSE106986 dataset, respectively. The AUCs of SLC27A3 and STAU1 in the metadata cohort were 0.734 and 0.745, respectively (Figures 5A, B). Furthermore, the AUCs of SLC27A3 and STAU1 in the GSE106986 dataset were 0.900 and 0.971, respectively (Figures 5C, D), indicating that both SLC27A3 and STAU1 have high diagnostic values.

With the CIBERSORT algorithm, we firstly calculated the proportion of immune cell infiltration in COPD and normal tissues. The results showed that the degree of infiltration of resting mast cells (P = 0.001), M1 macrophages (P = 0.049), and memory B cells (P = 0.020) in COPD tissues were notably raised in COPD samples in contrast to controls. Conversely, control samples demonstrated a higher proportion of infiltration of eosinophils (P = 0.002), activated mast cells (P = 0.006), resting NK cells (P = 0.004) and plasma cells (P = 0.003) compared to those found in COPD tissues (Figure 6A). Furthermore, we calculated the correlation between the 22 types of infiltrating immune cells (Figure 6B). Plasma cells and resting NK cells both individually correlated negatively with resting mast cells but positively with activated mast cells. Moreover, activated mast cells correlated positively with memory B cells.

SLC27A3 was significantly positively correlated with memory B cells (r=0.218, P=0.003), CD4 memory resting T-cells (r=0.245, P<0.001), CD8 T-cells (r=0.276, P<0.001), follicular T-helper cells (r=0.353, P<0.001), resting NK cells (r=0.406, P<0.001), eosinophils (r=0.411, P<0.001), activated mast cells (r=0.598, P<0.001) and plasma cells (r=0.619, P<0.001), and significantly negatively correlated with naive B cells (r=-0.657, P<0.001), CD4 memory activated T-cells (r=-0.578, P<0.001), resting mast cells (r=-0.388, P<0.001), resting dendritic cells (r=-0.289, P<0.001), activated NK cells (r=-0.274, P<0.001) and gamma delta T-cells (r=-0.169, P=0.022; Figure 7A). STAU1 was significantly positively correlated with regulatory T-cells (r=0.157, P=0.033), CD4 memory resting T-cells (r=0.174, P=0.019), CD8 T-cells (r=0.237, P=0.001), memory B cells (r=0.244, P<0.001), follicular T-helper cells (r=0.327, P<0.001), eosinophils (r=0.369, P<0.001), resting NK cells (r=0.468, P<0.001), activated mast cells (r=0.548, P<0.001) and plasma cells (r=0.575, P<0.001), and significantly negatively correlated with naive B cells (r=-0.590, P<0.001), CD4 memory activated T-cells (r=-0.513, P<0.001), resting mast cells (r=-0.422, P<0.001), resting dendritic cells (r=-0.230, P=0.002), activated NK cells (r=-0.299, P<0.001), gamma delta T-cells (r=-0.193, P=0.009) and M2 macrophages (r=-0.156, P=0.035; Figure 7B).

Subsequently, we established a COPD mouse model using GSE. The body weight of mouse in COPD group was significantly lower than that in control group (Figure 8A). HE staining showed that the MAST was shorter, and the MLI was larger in COPD group compared with that in control group (Figures 8 B–E). We then observed the immune cell infiltrating in COPD group. The proportion of leukocyte, and the percentage of neutrophil and macrophage were larger in COPD group (Figure 8F). Simultaneously, the expression levels of IL-6, IL-1β, and TNF-α were all increased in bronchoalveolar lavage fluid (BRLF) and lung tissues from COPD group (Figures 8 G–I). In addition, the expression of immune cell markers including CD8a, CD57, CD19, and CD20 were upregulated in lung tissues from mouse in COPD group (Figures 8J, K).

GSE-induced COPD models were established both in vivo and in vitro. Immunohistochemistry staining, real-time PCR and western blot assay showed that the mRNA and protein levels of SLC27A3 and STAU1 were upregulated in lung tissues of mice from COPD group (Figures 9 A–F). In accordance with these results, SLC27A3 and STAU1 were upregulated in GSE-treated BEAS-2B cells (Figures 9 G–J). Besides, ALCPL2 was also upregulated in COPD models, while RABL4 was not significantly changed (Figure S1).

Finally, we observed the effect of SLC27A3 and STAU1 knockdown on CSE-treated BEAS-2B cells. There siRNAs sequences were designed to downregulate SLC27A3 and STAU1, respectively. The expression of SLC27A3 and STAU1 were significantly downregulated by transfection of corresponding siRNAs (Figures 10 A–D). CSE treatment decreased the viability and increased the expression inflammatory cytokines (IL-6, IL-1β, and TNF-α) in BEAS-2B cells, which were attenuated by knockdown of SLC27A3 or STAU1 (Figures 10 E–H).