Exploring the relationship between the interleukin family and lung adenocarcinoma through Mendelian randomization and RNA sequencing analysis

3.1 Six interleukin factors in plasma are associated with the occurrence of LUAD

A total of 189 datasets related to interleukin factors were collected from the IEU database, comprising 79 interleukin factors. (Table S1). The heterogeneity analyses (Table S3), and horizontal pleiotropy analyses (Table S4) showed that the heterogeneity and horizontal pleiotropy were not significant. The causal effects of each protein on LUAD were calculated. Among these 79 proteins, elevated IL6RB (OR: 1.35, 95% CI 1.03 − 1.77, p-value = 0.03), IL27RA (OR: 1.17, 95% CI 1.03 − 1.34, p-value = 0.02), IL22RA1 (OR: 1.21, 95% CI 1.02 − 1.43, p-value = 0.03), and IL16 (OR: 1.14, 95% CI 1.02 − 1.28, p-value = 0.02) in plasma were risk factors for LUAD. Conversely, the plasma concentrations of IL18R1 (OR = 0.92, 95% CI 0.86–0.99, p-value = 0.04) and IL11RA (OR = 0.73, 95% CI 0.57–0.92, p-value = 0.007) were protective factors against LUAD. (Fig. 1). The effects of other interleukin factors on LUAD did not reach statistical significance. (Table S5). Leave one out test showed that the MR results obtained by us are reliable, but rs35026308 is strongly affected when calculating the effect of IL27RA on LUAD, which may cause unrobust results. (Figure S1).

Fig. 1figure 1

The effects of interleukin factors on lung adenocarcinoma

The figure showed six interleukin factors with significant causal effects on LUAD onset, the association effects between each interleukin factor and LUAD were calculated by three methods. Inverse variance weighted effect is considered as the main estimate effect. IL11RA and IL18R1 are considered protective factors for lung adenocarcinoma (shown in blue), whereas IL16, IL6RB, IL22RA1, and IL27RA are considered risk factors (shown in red).

3.2 IL18R1, IL11RA and IL16 were downregulated in LUAD tissue

Differential expression analyses were also conducted on the aforementioned 6 mRNAs across the 4 different datasets. The results revealed that in all four datasets, IL18R1, IL11RA, and IL16 were significantly downregulated in tumor tissues. For IL6RB, IL22RA1, and IL27RA, the differential expression results varied across different cohorts. (Fig. 2A–D).

Fig. 2figure 2

Differential analyses and survival validation of six interleukins/interleukin receptors via transcriptomics. AD Differential expression of the six interleukins/interleukin receptors in the (A) TCGA-LUAD cohort. B GSE31210 cohort. C GSE116959 cohort. D GSE102511 cohort. The red column indicates the expression level of tumor tissue, and the blue column indicates the expression level of normal tissue. The p value of the difference between tumor and normal tissue expression was above the column. E Survival differences in the TCGA-LUAD cohort based on IL11RA expression. F Differences in survival in the GSE31210 cohort stratified according to IL11RA expression. G Survival differences in the GSE72094 cohort based on IL11RA expression. H Survival differences in the TCGA-LUAD cohort based on IL18R1 expression. I Differences in survival in the GSE31210 cohort stratified according to IL18R1 expression. J Differences in survival in the GSE72094 cohort stratified according to IL18R1 expression. K Survival differences in the TCGA-LUAD cohort stratified according to IL16 expression. L Differences in survival in the GSE31210 cohort stratified according to IL16 expression. M Differences in survival in the GSE72094 cohort stratified according to IL16 expression

3.3 High expression of IL11RA indicates a better prognosis

Survival analysis revealed that, in LUAD patients across three different cohorts, patients with higher expression of IL11RA had better overall survival than did those with lower expression of IL11RA. (Fig. 2E–G) In the GSE31210 and GSE72094 cohorts, patients with higher IL18R1 expression and higher IL16 expression exhibited improved prognosis, while there was no difference in survival among patients in the TCGA-LUAD cohort. (Fig. 2H–M). Then stratified analyses were employed to take into account the effects of other factors on overall survival of patients with lung adenocarcinoma, the results showed a robust association between IL11RA and overall survival, whereas the association between IL16/ IL18R1 and overall survival in patients with lung adenocarcinoma was disrupted by multiple factors. (Figures S2, S3, and S4).

3.4 Western blotting

According to the above results, only the expression pattern of IL11RA and its relationship with prognosis were validated in multiple cohorts, so we further explored the role of IL11RA in lung adenocarcinoma.

To validate the downregulation of IL11RA not only at the RNA level but also at the protein level in tumor tissues, Western blotting (WB) was employed to assess the protein expression of IL11RA in LUAD cell lines and HBEC cell lines. The results revealed that, in three lung adenocarcinoma cell lines (A549, H1299, and H441), the protein levels of IL11RA were lower than those in HBE cells. However, compared with HBE cells, the PC-9 cell line showed no significant difference in IL11RA expression. (Fig. 3A). The full blot images are shown in Figure S5.

Fig. 3figure 3

Differential analyses and survival analyses of driver gene negative (DG-) patients. A Strip diagram and bar diagram of semi-quantitative results of the western blot illustrating the differences in the protein expression of IL11RA in HBEC cells and four different lung adenocarcinoma cell lines (H1299, H441, A549, PC9). B Differential expression of IL11RA in the GSE31210 DG- cohort. C Differential expression of IL11RA in the FAH cohort. D Survival differences in the GSE31210 DG- cohort based on IL11RA expression. E Differences in survival in the GSE72094 DG- cohort stratified according to IL11RA expression. F Differences in survival in the FAH cohort stratified by IL11RA expression

3.5 The relationship between IL11RA and prognosis was conserved in driver gene-negative patients

We compared the expression of negative driver genes between tumor tissue and normal tissue, and the results showed that IL11RA was still significantly downregulated in driver gene-negative patients. (Fig. 3B and C). The difference in survival between the high-low IL11RA expression group was compared in driver gene-negative patients: EGFR/KRAS-negative patients from the GSE72094 cohort, EGFR/KRAS/ALK-negative patients from the GSE31210 cohort, and pan-driver gene-negative patients from the FAH cohort. Remarkably, in all three cohorts, patients with high IL11RA expression continued to exhibit significantly better overall survival. (Fig. 3D–F).

3.6 Patients with more malignant LUAD tended to have lower expression of IL11RA

Differential analyses were also conducted between patients with different disease stages from both the TCGA-LUAD and FAH cohorts. In the TCGA-LUAD cohort, IL11RA expression was significantly downregulated in stage III patients compared to stage I and II patients. (Fig. 4A). Moreover, T2-stage patients exhibited lower IL11RA expression than T1-stage patients, and T3-stage patients had even lower IL11RA expression than T2-stage patients. (Fig. 4B). N2-stage patients had lower IL11RA expression compared to N1- and N0-stage patients. (Fig. 4C). However, no significant differences in IL11RA expression were observed between M0-stage patients and M1-stage patients. (Fig. 4D).

Fig. 4figure 4

Differential expression of IL11RA among patients at different stages. AD Differential expression of IL11RA among patients with different clinical stages/T stages/N stages/M stages in the TCGA-LUAD cohort. EH Differential expression of IL11RA among patients with different clinical stages/T stages/N stages/M stages in the FAH cohort

According to the FAH cohort from our center, stage IV patients exhibited significantly lower IL11RA expression than stage I patients. (Fig. 4E). Additionally, T4-stage patients exhibited significantly lower IL11RA expression compared to T3-, T2-, and T1-stage patients. (Fig. 4F). Similarly, N3-stage patients had significantly downregulated IL11RA expression compared to N2- and N0-stage patients. (Fig. 4G). Although there were no significant differences in IL11RA expression between M0-stage patients and M1-stage patients, a significant negative correlation was observed between IL11RA expression and the number of circulating tumor cells in patients. Therefore, further assessment of its correlation with metastatic risk is warranted. (Figs. 4H and 5A).

Fig. 5figure 5

Cell infiltration analysis results. A The correlation between IL11RA expression and the number of circulating tumor cells in the FAH cohort. B Correlations between IL11RA expression and stromal score, immune score, and ESTIMATE score in the TCGA-LUAD/GSE31210/GSE72094/FAH cohorts. C The correlation between IL11RA expression and 22 immune cell infiltration scores in the four cohorts. D The correlation between IL11RA expression and 64 immune cell and stromal cell infiltration scores in the four cohorts

We then performed Chi-square tests/Fisher’s exact test to compare the staging of patients with high or low IL11RA expression. There was no significant difference in the distribution of patients with high or low expression of IL11RA at the M0 and M1 stages, (p-value = 0.17). The same results were obtained in the FAH cohort (p-value = 0.29). (Table S6). However, it was observed that the T-stage distribution of patients with high or low expression of IL11RA was significantly different, whether in TCGA-LUAD (p = 4.3E-03) or FAH cohort (5.9E-03) (Table S6), which implied that higher IL11RA expression was associated with lower T staging.

3.7 Cell infiltration analyses

To explore how IL11RA impacts LUAD development, we computed immune scores and various cell infiltration scores for tumor patients from the TCGA-LUAD, GSE31210, GSE72094, and FAH cohorts. The estimated scores revealed a positive correlation between IL11RA expression and both stromal and immune scores across all four cohorts, resulting in decreased tumor purity (Fig. 5B).

The CIBERSORT results revealed significant heterogeneity in immune cell infiltration across the different cohorts. Notably, CD4 +memory-activated cells, M0 macrophages, and neutrophils were consistently negatively correlated with IL11RA expression (Fig. 5C).

The Xcell results indicated a consistent positive correlation between IL11RA expression and CD4 +T cell, fibroblast, mast cell, and iDC infiltration levels but a consistent negative correlation with epithelial cell and keratinocyte infiltration levels. These findings shed light on how IL11RA may influence the immune microenvironment in LUAD (Fig. 5D).

3.8 Co-expression and enrichment analyses

With respect to the four LUAD cohorts, we computed Pearson correlation coefficients between the expression levels of each gene and the expression level of IL11RA. We identified 538 genes that showed a positive correlation with IL11RA expression across all four cohorts. Further functional analyses through GO and KEGG pathway enrichment revealed that genes positively correlated with IL11RA expression were enriched in pathways related to GTPase-mediated signal transduction, such as the regulation of the RAS signaling pathway and cell–matrix adhesion. (Fig. 6A–C).

Fig. 6figure 6

IL11RA co-expression gene and enrichment analyses. A The intersection of genes positively correlated with IL11RA expression in the four cohorts. B GO enrichment results for genes positively correlated with IL11RA expression. C KEGG enrichment results for genes positively correlated with IL11RA expression. D The intersection of genes negatively correlated with IL11RA expression in all four cohorts. E GO enrichment results for genes negatively correlated with IL11RA expression. F KEGG enrichment results for genes negatively correlated with IL11RA expression

Conversely, 263 genes exhibited a negative correlation with IL11RA expression. These genes were enriched primarily in processes related to glycosylation and the HIF-1 signaling pathway. These findings provide insights into the potential biological mechanisms and pathways associated with IL11RA in LUAD. (Fig. 6D–F).

3.9 Single-cell sequencing analyses

To further explore the expression profile and function of IL11RA, we merged single-cell data from the GSE131907, GSE171145, and GSE189357 cohorts. After performing PCA dimension reduction and clustering, these cells were categorized into 29 clusters. (Fig. 7A). Clustering analyses revealed significantly greater IL11RA expression in normal tissues than in tumor tissues. (Fig. 7B and C).

Fig. 7figure 7

Results of single-cell RNA sequencing data clustering. A Uniform manifold approximation and projection (UMAP) plot of a total of 154,476 cells colored according to different clusters. B UMAP plot colored according to tissue type. C The expression of IL11RA in different tissues. D UMAP plot colored according to different cell types. E The expression of IL11RA in different cell types. F UMAP plot of 5348 fibroblasts classified by tissue type. G IL11RA + fibroblasts in tumor and normal tissues

After annotation, these cells were classified into 10 major cell categories. (Fig. 7D). IL11RA was expressed at higher levels in lung alveolar cells than in cancer cells, which is consistent with the findings obtained through Western blotting analyses. Furthermore, fibroblasts were the primary cell population expressing IL11RA (Fig. 7E). Given the known association between IL11RA and fibrosis, we performed subgroup analyses on 5348 fibroblasts extracted from the dataset. The results showed that there was a greater proportion of IL11RA+fibroblasts in normal tissues (30%) than in tumor tissues (13.5%) (Fig. 7F and G).

Trajectory analyses divided the fibroblasts into 3 states. (Fig. 8A). The differentiation proceeded from state 1 to state 2 and state 3. (Fig. 8B). IL11RA+fibroblasts were primarily located in state 1, representing an initial fibroblast state. (Fig. 8C and D). The proportion of IL11RA-expressing fibroblasts significantly decreased upon differentiation toward both ends. Thus, IL11RA+fibroblasts may represent a unique subgroup of fibroblasts and potentially serve as a protective factor for patients.

Fig. 8figure 8

Trajectory analyses of fibroblasts. A The trajectory analyses categorize fibroblasts into three states. B The evolutionary process of fibroblasts in different states. C IL11RA + and IL11RA- fibroblasts in each state. D The differential expression of IL11RA in fibroblasts across different states. E Results of GO enrichment analyses of genes differentially expressed between IL11RA + and IL11RA- fibroblasts. F KEGG enrichment analysis results for genes differentially expressed between IL11RA + and IL11RA- fibroblasts

Subsequently, we conducted differential analyses between IL11RA+and IL11RA- fibroblasts, which revealed 227 differentially expressed genes. Functional enrichment analyses revealed that these DEGs were enriched primarily in processes related to epithelial cell differentiation and extracellular matrix construction. (Fig. 8E). KEGG analyses revealed enrichment in several classical pathways. These findings provide insights into the potential role of IL11RA+fibroblasts in lung adenocarcinoma and their association with protective mechanisms and relevant biological processes (Fig. 8F).

Comments (0)

No login
gif