Abstract

Background

   The prognosis of lung adenocarcinoma (LUAD) is poor. Infection with
   coronavirus disease 2019 (COVID-19) may further worsen the outcome of
   LUAD. This study utilized the immune model and the COVID-19 receptor
   signal to identify the potential immune structure affecting the
   prognosis of COVID-19 and LUAD.

Methods

   A prognostic model was established and verified. The correlation
   between immune cells and risk score was examined through a variety of
   immune calculation methods. Gene set variation analysis (GSVA) was used
   to explore the correlation between the immune signaling pathway, risk
   model, COVID-19 binding receptor (CO19ORS) signal, and different
   clinicopathological factors.

Results

   The analysis showed that the prognosis of patients was better in the
   low-risk group versus the high-risk group. The tertiary lymphoid
   structure dominated by T and B cells (TLS1) can improve the prognosis
   of patients in the low-risk group. Interestingly, the CO19ORS was
   enriched only in females and aged >65 years. The age group >65 years is
   closely related to the tertiary lymphatic structure of the newborn
   (TLS2), while the female sex is closely related to the TLS2 and TLS1
   signature. The two groups exhibited a high level of
   inflammation-related signal distribution. In the near future, I will
   collect LUAD and COVID-19 related organizations to verify the changes
   of 8 risk protein.

Conclusion

   TLS1 structure may improve the prognosis of patients with LUAD and
   SARS-COV-2 (Severe acute respiratory syndrome coronavirus 2). This
   unexpected discovery provides new insight into the comprehensive
   treatment of patients with LUAD and SARS-COV-2.

   Keywords: LUAD, SARS-COV-2, TLS1, Immune prognostic model, COVID-19
   binding receptor, GSVA

1. Introduction

   According to the latest data published by the International Cancer
   Research Center in 2020, lung cancer is one of the most common
   malignancies worldwide. Moreover, it is associated with the highest
   rates of morbidity and mortality among all cancers [[33]1]. Non-small
   cell lung cancer (NSCLC) is the primary type of lung cancer. LUAD is a
   predominant histological type of NSCLC. Owing to the early development
   of metastatic disease, and despite the use of molecular pathology
   testing and targeted therapies, the prognosis of patients with LUAD is
   poor (average 5-year survival rate: <20%) [[34][2], [35][3], [36][4]].
   Since 2019, the coronavirus disease 2019 (COVID-19) pandemic poses a
   serious threat to the health and life of individuals worldwide. The
   infection is caused by a severe acute respiratory syndrome coronavirus
   2 (SARS-COV-2) and is mainly characterized by fever and inflammation of
   the lungs [[37]5,[38]6]. Therefore, it is important to consider the
   impact of infection with SARS-COV-2 on the prognosis of patients with
   LUAD.

   Furthermore, accumulating evidence suggests that the immune
   microenvironment significantly influences the progression of tumors and
   infection with SARS-COV-2, as well as the effectiveness of
   immunotherapy [[39][7], [40][8], [41][9], [42][10]]. Currently, the
   therapeutic targeting of programmed cell death 1 (PD1) and programmed
   cell death 1 ligand 1 (PDL1) reduces tumor volume and improves the
   prognosis in multiple types of cancer [[43][11], [44][12], [45][13]].
   However, these treatments do not play a determinative role in
   immunotherapy for NSCLC, with effective treatment observed in only a
   limited subset of patients [[46]14]. Recent studies found that tertiary
   lymphoid structure (TLS) is an aggregate of immune T and B cells
   responding to immune stimulation. Mature TLS are often linked to B cell
   development, and the presence of TLS in tumors is also associated with
   prolonged survival in patients with various types of cancer [[47]15].
   Moreover, the potential impact of TLS on SARS-COV-2 has not been
   investigated thus far. Therefore, the evaluation of multiple immune
   mechanisms and development of new immune-related therapeutic targets is
   necessary to improve the prognosis of patients with LUAD and COVID-19.

   In this study, we used a variety of databases to detect significant
   alterations in the expression of immune function genes, and established
   a risk model closely related to the prognosis of patients with LUAD.

2. Materials and methods

2.1. Data source and preprocessing

   The RNA-sequencing data of 535 tumors and 59 normal samples, and the
   corresponding data from a clinical environment, were downloaded from
   the Cancer Genome Atlas (TCGA) database. The clinical features of
   interest were age, sex, tumor stage, TNM classification, survival time,
   and vital status. Transcriptome profiling of 364 patients with LUAD was
   performed using the [48]GSE42127, [49]GSE50081, and [50]GSE115002 data
   sets of the Gene Expression Omnibus (GEO) database. The [51]GSE42127
   and [52]GSE50081 data sets were used for validation, while
   [53]GSE115002 was used to identify significant differences in gene
   expression using the GEO2R online tools. In addition, the proteomic
   data of 102 normal and 109 tumor samples were obtained from the
   Clinical Proteomic Tumor Analysis Consortium (CPTAC) database.

2.2. Screening of immune-related differentially expressed genes (IMDEGS)

   According to the limma package and Wilcoxon test of R (R Vision 4.0.3),
   the differentially expressed genes (DEGs) were analyzed. The screening
   thresholds for TCGA and GEO databases were |log fold change [FC]| >1
   and false discovery rate <0.05 or P-values <0.05. Additionally,
   P-values <0.05 and |log FC| >0.5 were used as the identification
   criteria for significant proteins in the CPTAC database.

   Firstly, 1811 immune-related genes were obtained from the tumor immune
   evaluation resource (TIMER) website. Next, the intersection of
   significant genes and proteins in multiple databases and among
   immune-related genes was determined. Finally, 81 IMDEGS were obtained.

2.3. Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG)
analyses

   We used the Database for Annotation, Visualization, and Integrated
   Discovery (DAVID) to perform GO and KEGG pathway analyses online. These
   analyses were conducted to further investigate the molecular function,
   cellular composition, biological process, and signaling pathway of
   these IMDEGS.

2.4. Construction of the immune-related prognostic model

   Initially, we used the “survival” package in R to identify IMDEGS
   related to prognosis. Subsequently, we implemented least absolute
   shrinkage and selection operator (LASSO) regression analysis based on
   the “glmnet” and “survival” packages for dimensionality reduction. We
   also calculated the analysis risk coefficient through multivariate Cox
   regression analysis.

   The risk score was computed using the following formula: risk score =
   (expr IMDEG1 × Coef IMDEG1) + (expr IMDEG2 × Coef IMDEG2) + … + (expr
   IMDEGN × Coef IMDEGN), where Coef is the coefficient obtained through
   the LASSO Cox regression, and expr is the expression of the signature
   genes. The receiver operating characteristic (ROC) curve in the
   “survival Roc” package was used to analyze the clinical optimal risk
   threshold, which was used as the basis for the classification of
   patients into high- and low-risk groups.

2.5. Correlations between risk models and clinical features

   Based on the “limma” package, differences in the risk scores of groups
   with different clinicopathological factors were analyzed. The
   “survival” package was used to investigate differences in survival
   between the high- and low-risk groups according to different
   clinicopathological factors. During the same period, univariate and
   multivariate Cox regression and ROC curves were used to analyze the
   impact and value of the risk score and multiple clinical indicators on
   the prognosis of patients. The calibration chart and consistency index
   (C-index) were determined using the “RMS” package in R to calibrate and
   evaluate the nomogram of various combined models.

2.6. Immune landscape exploration

   The immune cell infiltration for each sample in the TCGA-LUAD data set
   was assessed through various calculation methods, including TIMER
   [[54]16], CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL, and
   EPIC. In addition, Spearman correlation (R) and “limma” packages were
   used to examine the association between risk score and immune
   infiltration.

2.7. HALLMARK enrichment and GO analysis in different risk groups

   Cancer HALLMARK enrichment and GO analysis was performed in different
   risk groups by gene set enrichment analysis (GSEA;
   [55]http://www.gsea-msigdb.org/gsea/index.jsp) and using DAVID tools
   ([56]https://david.ncifcrf.GOv/home.jsp). A P-value <0.05 denoted
   statistically significant differences.

   Signal gene set related to thymocyte selection associated high mobility
   group box protein (TOX), TLS1-2, and CD8 + T cell depletion.

   Newborn TLS related 12 chemokines (TLS2_Signature), T, and B cells
   clustered in TLS-related nine genes (TLS1_Signal) [[57][17], [58][18],
   [59][19]].

   The “ACUTE_VS_CHRONIC_LCMV_PRIMARY_IARY_INF_CD8_TCELL_UP” gene set can
   be considered as CD8 + T cell depletion related signal gene set, one
   hundred previously reported TOX-related signal genes were used and well
   verified in lung squamous cell carcinoma [[60][20], [61][21],
   [62][22]]. Seven genes were associated with COVID-19 binding receptor
   (CO19ORS) signaling [[63][23], [64][24], [65][25], [66][26]].

2.8. Gene set variation analysis (GSVA)

   Significant differences in the enrichment of signaling pathways in the
   high- and low-risk groups were analyzed using the GSVA package and
   multiple immune-related gene sets. The screening criteria were |log FC|
   >0.3 and P-values <0.05.

3. Results

3.1. Identification of immune-related core genes

   RNA sequencing data of 535 LUAD tissues (tumor group) from TCGA
   database ([67]https://portal.gdc.cancer.gov/) and 59 non-cancerous
   tissues (healthy group) were used to perform differential expression
   analysis. Through this analysis, we obtained 6776 DEGs ([68]Fig. 1a).
   Thereafter, data from the [69]GSE115002 data set and CPTAC database
   were analyzed for DEGs and proteins; we obtained 5479 genes and 2234
   proteins ([70]Fig. 1b and c). Approximately 1811 immune-related genes
   were downloaded from the TIMER database; the DEGs and proteins detected
   in TCGA, CPTAC, and GEO databases overlapped. Finally, the analysis
   yielded 81 immune-related core genes ([71]Fig. 1c).

Fig. 1.

   [72]Fig. 1
   [73]Open in a new tab

   Identification of the differentially expressed immune genes (IMDEGS)
   between the three databases. Volcanic diagram of DEGs based on the
   comparison of Tumour/Normal group from TCGA database (a–b).
   Differential protein expression levels in different groups of the CPTAC
   database were determined by volcano map (b). Differential genes of
   multiple databases (TCGA, [74]GSE115002 and CPTAC) and the
   immune-related genes were taken intersection (c). d-e, IMDEGs were
   submitted to GO (d) and KEGG (e) enrichment.

3.2. GO enrichment and KEGG pathway analyses

   GO and KEGG pathway analyses were performed using the DAVID website to
   further recognize the feature of IMDEGS; the results were visualized
   using R. In the GO analysis of IMDEGS, the results revealed enrichments
   in biological process (immunoreactivity, positive regulation of
   angiogenesis, hypoxia, inflammatory, immune response, positive
   regulation of endothelial cell proliferation, negative chemotaxis,
   chemokine mediated signaling pathway, and chemotaxis), cellular
   composition (extracellular space, extracellular domain, extracellular
   body, protein, and extracellular matrix), and molecular function
   (retinoic acid-binding, chemokine activity, retinal binding,
   retinol-binding, transporter activity, heparin-binding, growth factor
   activity, signal receptor binding, and cytokine activity) ([75]Fig.
   1d).

   The KEGG analysis demonstrated that IMDEGS were mainly enriched in axon
   guidance, rheumatoid arthritis, cytokine receptor interaction,
   chemokine signaling pathway, phagosome, and the tumor necrosis factor
   (TNF) signaling pathway. The results showed that IMDEGS may be closely
   associated with immune reactions ([76]Fig. 1e).

3.3. Construction and analysis of the immune-related prognostic model

   Univariate Cox regression analysis was performed to investigate 26
   prognostic IMDEGS ([77]Fig. 2a). Risk genes and coefficients were
   determined using LASSO regression analysis and multivariate Cox
   regression; finally, eight risk genes were obtained ([78]Fig. 2b and
   c). The expression levels of these eight genes were used to construct
   the immune-related prognostic model. The prognostic risk score was
   calculated using the following formula: risk score = (−0.185 × expr LIF
   receptor subunit alpha [LIFR]) + (−0.197 × expr Fetal growth
   restriction [FGR]) + (−0.127 × expr WAP four-disulfide core domain 2
   [WFDC2]) + (0.164 × expr semaphorin 4B [SEMA4B]) + (0.728 × expr
   Fibroblast growth factor 2 [FGF2]) + (0.151 × expr C-C motif chemokine
   ligand 20 [CCL20]) + (−0.211 × expr growth differentiation factor 10
   [GDF10]).

Fig. 2.

   [79]Fig. 2

   [80]Fig. 2
   [81]Open in a new tab

   The technological process of constructing the immune-related risk
   model. The prognosis-related IMDEGS were performed by univariate
   analysis (a). The figures of the LASSO method (most minor absolute
   shrinkage and selection operator) and coefficient profiles (b).
   Multivariate Cox regression analyses were performed on the forest plots
   (c). Optimal thresholds points were obtained by the operating
   characteristic (ROC) curve (d). the distribution of the risk score (e).
   The survival status of patients was plotted in different risk groups
   and scores (f–h). The analysis of the time-dependent ROC curve was
   performed for patients in 1,2,3 years (i). Kaplan-Meier curves of the
   correlative models of immune in TCGA (training set) (j), [82]GSE42127
   and [83]GSE50081 (testing set 1, 2) (k–l), The distribution of
   different genes in high/low-risk groups was displayed as
   box-and-whisker plot (m). Association between eight prognosis-related
   genes with risk score (n).

   According to the analysis performed using the “survival rock” package
   in R ([84]Fig. 2d) the best-predicted risk coefficient of the risk
   model was 1.025; this coefficient was used as the basis for the
   classification of patients into high- and low-risk groups. Also, our
   data showed that the mortality rate was higher in the high-risk group
   than the low-risk group ([85]Fig. 2e–h). Kaplan–Meier survival curves
   for TCGA training data set indicated that the survival rate was higher
   in the low-risk group versus the high-risk group ([86]Fig. 2i). In the
   training set, the areas under the curve (AUCs) of the time-dependent
   ROC curve were 0.697, 0.712, and 0.716 for 1, 2, and 3-year OS
   predictions, respectively ([87]Fig. 2j). These findings confirmed that
   the immune-related risk model was sufficiently accurate to predict the
   prognosis of patients with LUAD. After performing a GEO data set
   analysis, we also conducted a validation study using the validation
   groups established by the same risk model and scoring formula to
   validate the stability of the model. The predictive efficacy of the
   model was validated using validation set1 and set2 by Kaplan–Meier
   analysis ([88]Fig. 2k and l). The expression levels of LIFR, FGR,
   WFDC2, CD79A, and GDF10 was higher in the low-risk groups than the
   high-risk group. There was a negative association with the risk score,
   contrary to the expression of SEMA4B and CCL20 in different risk groups
   and risk scores ([89]Fig. 2m, n).

3.4. Associations between risk score and clinicopathological factors

   The interrelationships between the risk score and clinical features of
   patients with LUAD were evaluated. The risk score was associated with
   sex (P = 0.003), clinical stage (P = 7.7e-07), T classification
   (P = 0.0016), and N classification (P = 0.00012) ([90]Fig. 3a–f).
   Meanwhile, according to the survival analysis of different risk groups
   based on the clinical characteristics, we found that the high-risk
   group (age: >65, ≤65 years; sex: male and female; stage: I–II, M0,
   T1–2, T3–4, and N1-2) was characterized by a low survival rate
   ([91]Fig. 3g).

Fig. 3.

   [92]Fig. 3

   [93]Fig. 3
   [94]Open in a new tab

   The correlation of the various clinical features with the Risk Score.
   Spearman was performed the association between risk score and different
   clinical parts (a–f). The difference of survival between high and
   low-risk groups in other clinical characteristics (g).

   Univariate and multivariate Cox regression models identified that TMN
   stage and risk scores were significantly associated with prognosis
   ([95]Fig. 4a and b). The ROC curve results showed that the predictive
   performance of the model for the rate of overall survival was
   excellent. [96]Fig. 4c shows the AUC values of the various clinical
   features and risk scores: age (AUC = 0.471), sex (AUC = 0.558), TMN
   stage (AUC = 0.669), T classification (AUC = 0.576), M classification
   (AUC = 0.496), N classification (AUC = 0.632), and risk score
   (AUC = 0.697). Subsequently, we constructed a nomogram using the risk
   scores and TMN stages for the prediction of the patient ([97]Fig. 4d,
   e). The C-index of the nomogram was 0.69 (95% confidence interval:
   0.6468–0.73312), and a calibration plot showed optimal agreement
   compared with the actual observation.

Fig. 4.

   [98]Fig. 4
   [99]Open in a new tab

   Effect of risk model and various clinicopathological features on the
   prognosis of patients. The immune-related prognostic model and Clinical
   features were analyzed via the Univariate, multivariate Cox, and
   receiver operating curve (ROC) (a–c). Based on the immune-related model
   and more prominent clinical factors, a nomogram was built for
   predicting 1-, 2- and 3-year survival rates of LUAD, and the prognostic
   accuracy of the nomogram was shown by calibration plots (d–e).

3.5. Composition of immune cells in the risk model

   We estimated the immune cell composition of TCGA-LUAD through multiple
   calculation methods from the TIMER website ([100]Fig. 5a) to analyze
   the association between risk score and immune cells. The risk score
   exhibited a generally negative correlation with B cells and CD8 T
   cells. There were significant differences in B, CD8, and myeloid
   dendritic cells between the low- and high-risk groups; B cell
   infiltration was the most prominent among them ([101]Fig. 5b–e). The
   distribution differences in the immune score, matrix score, and risk
   score, were calculated by x cell under different risk groups. Based on
   the data, we concluded that risk score had a significant negative
   association with the immune score and stroma score, which were higher
   in the low-risk group versus the high-risk group ([102]Fig. 5f).

Fig. 5.

   [103]Fig. 5

   [104]Fig. 5
   [105]Open in a new tab

   Research of different risk groups of immune infiltration. Spearman
   correlation analysis was computed between risk scores and immune
   infiltration, using multiple calculation methods from TIMER Database
   (a–e). Distribution of the microenvironment score and composition of
   immune cells in different risk groups (f). (_A: TIMER,_B: CIBERSORT,
   _C: CIBERSORT-ABS, _D: QUANTISEQ, _E: MCPCOUNTER, _G: XCELL, _H: EPIC).

3.6. Distribution of immune-related genes in different risk groups

   The genes related to B cells (CD20, CD23, marginal zone B and B1 cell
   specific protein [MZB1], joining chain of multimeric IgA and IgM
   [JCHAIN], immunoglobulin lambda-like [IGLL0]), CD8 cells (interleukin 2
   [IL2], interleukin 4 [IL4], CD8, C-X-C motif chemokine receptor 4
   [CXCR4], C-X3-C motif chemokine receptor 1 [CX3CR1]), CD4 T cells
   (CD4), follicular dendritic cells (CD21), and primary T cells (TCF7)
   were highly expressed in the low-risk group ([106]Fig. 6a). We also
   found that hypoxia inducible factor 1 subunit alpha (HIF1A) was highly
   expressed in the high-risk group and positively associated with the
   risk score ([107]Fig. 6 b). CD8 cells, B cells, TOX, and TLS1
   signatures were higher in the low-risk group versus the high-risk
   group. These signals were negatively correlated with the risk score
   ([108]Fig. 6c, d).

Fig. 6.

   [109]Fig. 6
   [110]Open in a new tab

   Distribution level of tertiary lymphatic structure (TLS2, TLS1) in
   different risk groups. (a–b) The difference of the T, B, and dendritic
   cell-related genes in the high/low-risk group. (c) Based on GSVA
   analysis and GO data set in GSEA, (d) the distribution levels of CD8
   cells, B cells, TLS1, TLS2, CD8+T cell Exhaustion, and Tox-related
   signal pathways in high and low-risk groups and their correlation with
   the
   riskscore.(GO_1:B_CELL_HOMEOSTASIS,GO_2:B_CELL_PROLIFERATION,GO_3:B_CEL
   L_MEDIATED_IMMUNITY,GO_4:B_CELL_DIFFERENTIATION,GO_5:B_CELL_ACTIVATION,
   GO_6:B_CELL_RECEPTOR_SIGNALING_PATHWAY,GO_7:B_CELL_ACTIVATION_INVOLVED_
   IN_IMMUNE_RESPONSE,GO_8:CD8_POSITIVE_ALPHA_BETA_T_CELL_ACTIVATION,GO_9:
   CD8_POSITIVE_ALPHA_BETA_T_CELL_DIFFERENTIATION).

3.7. GO analysis of different risk groups

   The high-risk group exhibited enrichments in the immune-related
   pathways of cellular signal transduction, the defense response of
   Gram-negative bacteria, chemotaxis, bacterial defense response,
   neutrophil chemotaxis, and hypoxia response ([111]Fig. 7a). The
   low-risk group exhibited enrichments in the immune-related pathways of
   the antibacterial immune response, humoral reaction, complement
   activation, receptor-mediated endocytosis, classical pathway,
   leukotriene biosynthesis process, B cell receptor signaling pathway,
   and positive regulation of B cell activation ([112]Fig. 7b).

Fig. 7.

   [113]Fig. 7
   [114]Open in a new tab

   Enrichment analysis of highly expressed genes in different risk groups.
   GO analysis of the highly expressed genes of high/low-risk groups,
   respectively (a–b). KEGG enrichment analysis of upregulated genes of
   the high and low-risk groups, respectively (c–d).

3.8. Analysis of the KEGG pathways by GSEA

   Thereafter, we conducted a specific pathway enrichment analysis using
   GSEA for different risk groups according to the KEGG data set. The
   low-risk group was enriched in arachidonic acid metabolism, linolenic
   acid metabolism, and production of intestinal immune network ([115]Fig.
   7d). The high-risk group was enriched in DNA replication, p53 signaling
   pathway, cell cycle, glycolysis, galactose metabolism, gluconeogenesis,
   oocyte meiosis, and renal cell carcinoma ([116]Fig. 7c).

3.9. Association of CO19RS signaling with immunity

   Immune and stromal scores, TLS2, TLS1, and TOX signals were highly
   distributed in the high-CO19ORS expression group, and were directly
   proportional to the expression of CO19ORS ([117]Fig. 8a–d). For the
   clinicopathological factors, the distribution level of the CO19ORS
   signal in the age and sex groups differed; the highest concentration
   was observed in females and aged >65 years ([118]Fig. 8e and f). The
   expression of CO19ORS-related genes (AXL, epidermal growth factor
   receptor [EGFR], and neuropilin 1 [NRP1]) was higher in females versus
   males. The levels of other related genes (AXL receptor tyrosine kinase
   [AXL], angiotensin converting enzyme 2 [ACE2], and CD209) were higher
   in patients aged >65 years versus those aged ≤65 years ([119]Fig. 8k,
   l).

Fig. 8.

   [120]Fig. 8

   [121]Fig. 8
   [122]Open in a new tab

   Correlation between COVID-19 binding receptor, immune response and
   different pathological parameters. (a) Distribution of immune and
   matrix scores in high and low-COVID-19 receptor group. (b) Correlation
   between COVID-19 binding receptor and immunity and matrix score. (c)
   Distribution of immune-related signals and risk scores in high and low
   COVID-19 binding receptor groups. (d) Correlation between
   immune-related signals, risk score, and high and low COVID-19 combined
   receptor score. e-j Distribution level of COVID-19 binding receptors
   under different clinical and pathological factors. (e) age. (f) gender.
   (g) M. (h) N. (I) stage. (j) T. k-l, Expression levels of COVID-19
   binding receptors in different age and gender groups, (k) gender. (l)
   age.

3.10. Differences in the enrichment of TOX, TLS2, and TLS1 signatures between
groups

   The highest TLS2, TLS1, and TOX signals were observed in those aged >65
   years with high expression of CO19ORS (Group 1) and females with high
   expression of CO19ORS (Group 3) ([123]Fig. 9a–d). In Group 1, the TLS2
   and TOX signals were proportional to the CO19ORS signal ([124]Fig. 9e).
   In those aged ≤65 years with low expression of CO19ORS (Group 2), the
   CO19ORS signal was directly proportional to the TOX signal ([125]Fig.
   9f). In Group 3 and Group 4 (males with low expression of CO19ORS),
   there was a positive correlation between the signals of the CO19ORS,
   TLS2, TLS1, and TOX ([126]Fig. 9g and h).

Fig. 9.

   [127]Fig. 9
   [128]Open in a new tab

   Using GSVA analysis, the enrichment levels of TOX, TLS1, and TLS2
   signals in different groups Age>65-high-COVID-19 binding receptor
   (Group1) V Age ≤ 65-Low-COVID-19 binding receptor(Group2), and
   FEMALE-high-COVID-19 binding receptor(Group3) V MALE-Low-COVID-19
   binding receptor (Group4), a-b Group1v2. (a) heatmap. (b) Deviation
   diagram. c-d Group3v4. (c) heatmap. (d) Deviation diagram. e-h, The
   correlation degree between TLS2, TOX, TLS1 and COVID-19 binding
   receiver signals under different groups. (e) Group1. (f) Group2. (g)
   Group3. (h) Group4.

3.11. Inflammatory signaling pathways

   Based on the GSVA analysis and HALLMARK data set, we found that
   inflammation-related pathways, including the IL6/JAK/STAT3 signaling
   and inflammatory response, were highly expressed in Group 1 ([129]Fig.
   10a–c) and Group 3 ([130]Fig. 10d–f), compared with Groups 2 and 4.

Fig. 10.

   [131]Fig. 10
   [132]Open in a new tab

   Based on GSVA analysis, HALLMARK-related signal pathways are enriched
   in different groups Group1) V Group2, Group3 V Group4. a-c,
   Inflammation related signaling pathways were distributed at Group1. (a)
   Volcan, (b) Heatmap, (c) Deviation diagram. d-f, Inflammation related
   signaling pathways were distributed at Group3. (d) Volcan, (e) Heatmap,
   (f) Deviation diagram.

4. Discussion

   Lung adenocarcinoma is the most common pathological subtype of NSCLC
   (30–40%) [[133]27]. NSCLC poses a significant public health challenge,
   and is the leading cause of cancer-related death worldwide. The study
   and development of various therapeutic have suggested that the immune
   damage is decisive factors in the development of cancer. With the
   application of cancer immunotherapy, it has been shown that some
   immunomodulators also play an increasingly important role in the
   treatment of LUAD [[134][28], [135][29], [136][30], [137][31],
   [138][32]]. SARS-COV-2 current causes severe lung disease globally
   [[139]7]. This ongoing phenomenon undoubtedly complicates the treatment
   of LUAD. Therefore, we strive to improve the prognosis of patients with
   COVID-19 and LUAD through the development of new immunotherapies.

   In this study, we obtained 81 core IMDEGS based on the intersection of
   immune-related genes and DEGs reported in TCGA, GEO, and CPTAC
   databases. The IMDEGS exhibited enrichments in immune response,
   positive regulation of angiogenesis, hypoxia response, inflammatory
   response, negative chemotaxis, positive regulation of endothelial cell
   proliferation, chemokine mediated signaling pathway, chemotaxis, etc.
   According to the KEGG pathway analysis, the IMDEGS were predominantly
   enriched in rheumatoid arthritis, cytokine receptor interaction,
   phagosome, TNF signaling pathway, and the chemokine signaling pathway.
   The results showed that IMDEGS may be closely related to immune
   response. Eight genes significantly associated with prognosis were
   selected to create an immune-related prognostic model. As demonstrated
   by univariate and multivariate Cox regression analyses, as well as
   LASSO regression analysis, the model showed good prognostic accuracy.
   Moreover, we found that the high-risk group had a lower overall
   survival rate versus the low-risk group; of note, the classification
   into risk groups was based on the optimal threshold of the training
   (TCGA) and the validation ([140]GSE42127, [141]GSE50081) data sets. In
   addition, the relationship and prognostic value of the risk score with
   different clinical indicators was well verified. The nomogram – a
   statistical model based on the risk model and stage (the most prominent
   prognostic factors) – provided a more individualized prediction. In
   summary, the results demonstrated that the established immune-related
   model was stable and accurate. Meanwhile, we analyzed differences in
   the biological functions between the different risk groups. In the
   high-risk group, enrichment was noted in C-C motif chemokine receptor 6
   (CCR6) chemokine receptor binding, Gram-negative bacteria defense
   response, bacterial defense response, neutrophil chemotaxis, and
   hypoxia response. In the low-risk group, enrichment was observed in the
   signal of the receptor-mediated endocytosis, positive regulation and
   phagocytosis of B cell activation recognized leukotriene biosynthesis,
   immune response, complement activation, classical pathway, B cell
   receptor signaling pathway, antibacterial humoral response, and
   immunoglobulin receptor binding. Through the GO analysis of different
   risk groups, we established that the related paths of B cell and
   humoral immunity were ubiquitous in the low-risk group. We also found
   that the hypoxic environment may inhibit B cell and humoral response in
   the low-risk group. The high expression group significantly enriched
   hypoxia related pathways and genes. Recent studies have also suggested
   that hypoxia may affect the development of B cells and humoral response
   [[142]33].

   Using the TIMER database, we examined the correlation between immune
   cell infiltration and risk score in LUAD. We found that B cells, CD8
   cells, and TLS1-related signals were highly enriched in the low-risk
   group. Moreover, immune genes related to B cells, CD4 cells, CD8 cells,
   and dendrites were also highly expressed in the low-risk group,
   indicating that the low-risk group also had TLSs [[143][34], [144][35],
   [145][36], [146][37]]. Recent studies have found that extensive B cell
   infiltration, or accompanied by the TLS system of the CD20 + B cells in
   sarcoma, melanoma, renal cell carcinoma, and lung cancer contribute to
   the prognosis of patients [[147]19,[148]35,[149]38,[150]39]. Previous
   studies have shown that CD8 T cells play an essential role in the
   anti-tumor immune response and are positively correlated with the
   overall survival rate of patients with various malignant tumors
   (including lung cancer) [[151]40,[152]41]. Interestingly, the TLS1, TOX
   and CD8 T cell depletion signaling pathways were highly enriched in the
   low-risk group. T cell depletion leads to T cell dysfunction, which
   occurs in the tumor microenvironment. This process is often responsible
   for the inability of the immune system to eliminate the tumor
   [[153]42,[154]43]. Several investigations have reported that TOX
   regulates the critical factor of CD8 cell depletion [[155][44],
   [156][45], [157][46]]. This finding suggests that TLS1 counteracts the
   impact of CD8 cell depletion on the prognosis of patients.

   Among the eight genes included in the risk model, LIFR, FGR, WFDC2,
   GDF10, and CD79A were upregulated. In the low-risk group, SEMA4B and
   CCL20 were downregulated. FGF2 did not show statistically significant
   differences in any group. As a serum tumor marker, the upregulation of
   LIFR can inhibit the proliferation and invasion of lung cancer cells,
   and promote their apoptosis [[158]47]. WFDC2 plays essential roles in
   the early diagnosis of lung cancer and monitoring of treatment
   efficacy, and is significantly correlated with the TNM stage and
   prognosis of LUAD. Overexpression of WFDC2 can repress the metastatic
   potential of prostate cancer by inhibiting the EGFR signal
   [[159]48,[160]49]. SEMA4B restrains the growth of NSCLC in vitro, ex
   vivo, inhibited by hypoxia and HIF1A, promoting the invasion of NSCLC
   [[161][50], [162][51], [163][52]]. The downregulation of FGF2 inhibits
   the invasion of breast cancer, osteosarcoma, and LUAD [[164][53],
   [165][54], [166][55], [167][56]]. The overexpression of CCL20 promotes
   the progression of lung cancer and breast epithelial cells.
   Furthermore, it recruits regulatory T cells through forkhead box
   O1/CCAAT enhancer-binding protein/nuclear factor-kB (FOXO1/CEPB/NF-kB)
   signal transduction to promote drug resistance in colorectal cancer
   cells [[168][57], [169][58], [170][59]]. The downregulation of GDF10
   promotes the proliferation of prostate cancer and the invasion of
   hepatocellular carcinoma [[171]60,[172]61]. It has been shown that high
   CD79A expression was beneficial to the prognosis of cervical cancer
   [[173]62]. In addition, FGR expression is associated with decreased
   cytotoxic T-lymphocyte infiltration and poor prognosis in colorectal
   cancer. Other studies have demonstrated that FGR plays a role in the
   myeloid cell-mediated inflammatory response in vivo [[174]63,[175]64].

   Ordinary model construction is based on gene sequencing; immune-related
   genes are evaluated based on the intersection of differential genes and
   proteins. The survival results obtained from multiple data sets were
   consistent. Additionally, the model was associated with a good ROC
   curve and C-index, which confirmed its stability and accuracy. We
   downloaded the immune invasion data through various immune calculation
   methods to ensure the accuracy of immune cell scores for different risk
   groups. Moreover, the enrichment analysis of a variety of
   immune-related gene sets further verified the results of immune cell
   infiltration.

   The gene set of CO19ORS signal is composed of several potential and
   cellular receptors that overlap with those of COVID-19, including ACE2,
   AXL, EGFR, LDLR, NRP1, CD209, and CLEC4M [[176][23], [177][24],
   [178][25], [179][26]]. The signal was mainly noted in females and
   patients aged >65 years.

   There was a marked difference in the CO19ORS between the age and sex
   groups. Females and patients aged >65 years were linked to high
   enrichment in COVID-19 receptor signals. It is likely that such
   patients are susceptible to COVID-19. Also, there was a significant
   correlation between the COVID-19 signaling pathway and immune score,
   TLS2, TLS1, and TOX signals. These results suggested that females and
   patients aged >65 years may exhibit higher immunoreactivity. In Group
   1, the receptor signal was positively correlated with the TOX and TLS2
   signatures. In Group 2, the receptor was only related to the TOX
   signal. Inflammatory signals were mainly concentrated in Group 1
   compared with Group 2. In recent reports, the mortality and reinfection
   rates were higher in patients with COVID-19 aged >65 years versus those
   aged ≤65 years [[180]65,[181]66]. In Groups 3 and 4, the TOX, TLS2, and
   TLS1 signatures were positively correlated with the COVID-19 receptor
   signaling pathway. In Groups 3 and 4, inflammatory-related signals were
   mainly concentrated in Group 3. In numerous reports, SARS-COV-2 was
   associated with mild symptoms and a low mortality rate in females
   versus males [[182][67], [183][68], [184][69], [185][70], [186][71],
   [187][72]]. Hence, TLS1 may play an important role in the prognosis of
   female patients. The adaptive immune response dominated by T and B
   cells plays a positive role in resisting SARS-COV-2 infection
   [[188][73], [189][74], [190][75]]. The inflammatory reaction often
   increases the damage caused by SARS-COV-2 to patients
   [[191]76,[192]77]. Meanwhile, TLS1 is closely related to B cells and
   CD8 cells in TLS [[193]18,[194]19]. TLS2 is closely related to newborn
   TLS, and is accompanied by inflammatory reaction
   [[195]17,[196]78,[197]79]. Therefore, in Group 1, the close
   relationship of TLS2 with the CO19RS may increase the damage induced by
   SARS-COV-2 in patients aged >65 years. In Group 2, the close
   relationship of TLS1 with the CO19RS may reduce the
   inflammation-related damage and prevent the damage caused by SARS-COV-2
   infection.

5. Conclusion

   According to the present immune-related model, TLS1 may improve the
   prognosis of patients with LUAD. The differences observed in the
   correlation and enrichment level of CO19ORS, TLS2, and TLS1 in
   different age and sex groups further imply that TLS1 is associated with
   resistance to SARS-COV-2 infection. Therefore, sufficient construction
   of T cell- and B cell-based TLSs may improve the prognosis of patients
   with LUAD and COVID-19.

Author contribution statement

   Kang Sun: Performed the experiments; Wrote the paper.

   Zhiqiang Zhang; Dongqin Wang; Hongyu Ma: Analyzed and interpreted the
   data.

   Jing Zhang: Conceived and designed the experiments.

   Chaoqun Lian: Conceived and designed the experiments; Contributed
   reagents, materials, analysis tools or data.

Funding statement

   Chaoqun Lian was supported by the Key Natural Science Project of Anhui
   Provincial Education Department (KJ2021A0773), National Innovation
   Program for College Students (202210367076), Anhui Provincial
   Undergraduate Innovative Training Program (S202110367030). Jing Zhang
   was supported by the Key Natural Science Project of Anhui Provincial
   Education Department (KJ2020A0578). Dongqin Wang was supported by the
   Program for Graduate Research of Bengbu Medical College (Byycxz22016).

Data availability statement

   The authors do not have permission to share data.

Declaration of interest's statement

   The authors declare no conflict of interest

Contributor Information

   Jing Zhang, Email: jade.zhangjing@bbmc.edu.cn.

   Chaoqun Lian, Email: lianchaoqun@bbmc.edu.cn.

Abbreviation

   COVID-19
          Coronavirus disease 2019

   CO19ORS
          COVID-19 binding receptor

   CPTAC
          Clinical Proteomic Tumor Analysis Consortium

   DEGs
          Differentially expressed genes

   GEO
          Gene Expression Omnibus

   GO
          Gene Ontology

   Group 1
          Aged >65 years with high expression of CO19ORS

   Group 2
          Aged ≤65 years with high expression of CO19ORS

   Group 3
          Female with high expression of CO19ORS

   Group 4
          Male with low expression of CO19ORS

   IMDEGS
          Immune-related differentially expressed genes

   KEGG
          Kyoto Encyclopedia of Genes and Genomes

   LUAD
          Lung adenocarcinoma

   NSCLC
          Non-small cell lung cancer

   SARS-COV-2
          Severe acute respiratory syndrome-coronavirus-2

   TCGA
          The Cancer Genome Atlas

   TIMER
          Tumor Immune Estimation Resource

   TLS
          Tertiary lymphoid structures

   TLS1
          Tertiary lymphoid structure signal dominated by T and B cells

   TLS2
          Tertiary lymphatic structure of the newborn

References