Construction of shared gene signature between rheumatoid arthritis and lung adenocarcinoma helps to predict the prognosis and tumor microenvironment of the LUAD patients

Introduction: Rheumatoid arthritis (RA) is a common chronic autoimmune disease with high incidence rate and high disability rate. One of the top complications is cancer, especially lung adenocarcinoma (LUAD). However, the molecular mechanisms linking RA and LUAD are still not clear. Therefore, in this study, we tried to identify the shared genetic signatures and local immune microenvironment between RA and LUAD and construct a clinical model for survival prediction. Methods: We obtained gene expression profiles and clinical information of patients with RA and LUAD from GEO and TCGA datasets. We performed differential analysis and Weighted Gene Co-expression Network Analysis (WGCNA) to discover the shared genes between RA and LUAD. Then, COX regression and LASSO analysis were employed to figure out genes significantly associated with survival. qRT-PCR and Western blot were utilized to validate the expression level of candidate genes. For clinical application, we constructed a nomogram, and also explored the value of RALUADS in characterizing immune infiltration features by CIBERSORT and xCell. Finally, responses to different drug therapy were predicted according to different RALUADS. Results: Our analysis identified two gene sets from differentially expressed genes and WGCNA gene modules of RA and LUAD. Filtered by survival analysis, three most significant shared genes were selected, CCN6, CDCA4 and ERLIN1, which were all upregulated in tumors and associated with poor prognosis. The three genes constituted RA and LUAD score (RALUADS). Our results demonstrated that RALUADS was higher in tumor patients and predicted poor prognosis in LUAD patients. Clinical nomogram combining RALUADS and other clinicopathological parameters had superior performance in survival prediction (AUC = 0.722). We further explored tumor immune microenvironment (TME) affected by RALUADS and observed RALUADS was closely related to the sensitivity of multiple immune blockades, chemotherapy and targeted drugs. Conclusion: Our findings suggest that there are shared physiopathologic processes and molecular profiles between RA and LUAD. RALUADS represents an excellent prognosis predictor and immune-related biomarker, which can be applied to select potential effective drugs and for LUAD patients with RA.


Introduction
Rheumatoid arthritis (RA) is one type of chronic autoimmune diseases (AIDs) characterized by progressive joint degeneration and extra-auricular inflammation, leading to systemic destruction, decreased quality of life and even permanent disability (Smolen et al., 2016).According to epidemiological data, the global age-standardized prevalence rate of RA is 224.25, ranking the first in the major AIDs and demonstrated a significantly increasing trend from 1990 to 2019 (Cao et al., 2023).Although the etiology of this common AID is still uncertain, increased antibodies and immune-related substances such as rheumatoid factor and TNFα, is considered one of the most important features of RA.Therefore, drugs targeting the immune system, mainly disease-modifying antirheumatic drugs (DMARDs) is the standard treatment for RA with nanocarriers containing targeted drugs being investigated vigorously (Lin et al., 2020;Zhang et al., 2021;Smolen et al., 2023).However, in spite of the significantly improved response rate with the progression of RA treatment, increasing incidence of comorbidity and related mortality have been reported, including cardiovascular diseases (CVDs), infections, cancers etc. (England et al., 2018;Lortholary et al., 2020;Ytterberg et al., 2022;Chen et al., 2023).The underlying reasons may be systemic and chronic inflammation, genetic factors, immunosuppressant drugs and conventional risk factors such as high blood lipids level and other imbalance caused by disturbed metabolism (Szekanecz et al., 2011;Ytterberg et al., 2022).
Recent epidemiological surveys showed that the standardized incidence ratios of RA to develop cancer is 3.99 (95% CI = 3.40-4.65)(Zhou et al., 2022), with cancer types ranging from different hematological malignancies which rank the first, to other various solid tumors, such as lung cancer, melanoma, cervical cancer etc. (Zhang et al., 2022).Of note, lung adenocarcinoma (LUAD) was shown to be the most common lung cancers in AIDs (Jacob et al., 2020) and pulmonary involvement is a common extraarticular manifestation of RA (Sparks, 2019).Immune dysregulation is considered as a crucial contributor to the development of LUAD of RA patients.A bioinformatic study suggested that CD8A, GZMA, and PRF1 were related to CD8 + T cell in RA and positively associated with 33 tumors (Zhao et al., 2022).What's more, MAP4K3 which interacts with and activates PKC, resulted in the activation of IKK/ NF-kB in human T cells (Chuang et al., 2016).Except for T cells, Zhang etc. showed that extracellular ADP disrupted tissue homeostasis and function to trigger RA by recruiting neutrophils (Zhang et al., 2022).Nonetheless, the mechanisms under the association between RA and LUAD are still unclear yet.
Hence, in this study, we employed bioinformatic tools to demonstrate the shared mechanisms of RA and LUAD from the perspective of gene expression and clinical profiles, and immune infiltration characteristics.In order to apply to the clinical settings, we constructed a nomogram model and predict potential therapeutic drugs for LUAD, which we hope may help improve the early and effective management and treatment of LUAD.

Data download and processing
In order to obtain RA expression data, we applied the keyword "rheumatoid arthritis" in the Gene Expression Omnibus (GEO) dataset (GEO, https://www.ncbi.nlm.nih.gov/geo/).We search according to the following criteria: (Smolen et al., 2016): the dataset must be complete and correct, (Cao et al., 2023) samples contained in the dataset must be enough; (Lin et al. , 2020) samples were collected from joint lesions rather than blood etc.Finally, we chose GSE236924, an expression profiling by array, which is normalized using RMA and contains 36 RA samples and 7 normal samples in total.The mRNA expression data of LUAD patients with corresponding clinical information were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/)database with 541 cancer patients and 59 normal patients contained.During analysis, patients with missing information were excluded.In addition to the above two discovery cohorts, we selected another LUAD dataset, GSE229705, from GEO with complete gene expression profile data and clinical information as validation cohort, which contained paired tumor tissues and adjacent normal tissues from 123 patients.Detailed clinical information of two LUAD cohorts are shown in Supplementary Tables S1, S2.

Bulk gene expression data analysis
After importing the datasets and group information, we first conducted principal component analysis (PCA) of all samples in each dataset to ascertain that samples with different phenotypes demonstrate different genotypes.Then, after quality control, we used R package "Limma" and "DESeq2" to perform differential analyses for RA array dataset and LUAD RAN-seq dataset respectively.Differentially expressed genes (DEGs) were computed by setting the parameters of (|log2FC| > 1, adjusted P-value <0.05).DEG results were demonstrated as volcano plots.Further, enrichment analyses, including Kyoto Encyclopedia of Genes and Genomes (KEGG), gene ontology (GO) and Gene Set Enrichment Analysis (GSEA) were performed based on "clusterProfiler", "org.Hs.eg.db", "enrichplot" packages" to examine the significant pathways involved in DEGs.Finally, we identified the shared differentially expressed gene signatures via Venn gram.

WGCNA analysis
We used the Weighted Gene Co-Expression Network Analysis (WGCNA) tool to evaluate gene expression profiles in RA and LUAD with "WGCNA" package (Langfelder and Horvath, 2008).After removing the missing values and outliers, we constructed the co-expression networks related with clinical traits according to the best soft threshold, and chose the core modules with the highest Pearson coefficients.Hub genes in both RA and LUAD core modules were then identified and the shared gene signatures were discovered via Venn gram.

Construction and validation of RALUADS
For the 22 shared DEGs, we performed univariate Cox regression analysis and identified one gene with P-value <0.05.For the 62 shared WGCNA genes, we first applied the least absolute shrinkage and selection operator (LASSO) analysis for dimension reduction.After model fitting and cross validation, 29 candidate genes were finally selected according to lamda.1se.Then, we successively performed univariate and multivariate Cox regression analysis to screen out two potential genes most related with overall survival (OS).Based on the above results, we established a three-gene prognostic signature termed RA and LUAD signature (RALUADS).RALUADS for each LUAD patient is calculated by the following formula: Expression refers to the expression level of the selected gene, and Coefficient is the coefficient of the selected gene in the univariate Cox regression model.LUAD patients were divided into two groups according to the median RALUADS.Kaplan-Meier curves were drawn to assess the OS differences between different RALUADS groups.COX regression was analyzed with "survival" package and visualized with "forestplot" package.LASSO regression was performed via "glmnet" package.

Construction and evaluation of clinical model
We established a practical clinical model using "nomogram" package after screening out the meaningful clinical features via multivariate COX regression analysis.Then, we compared the prognostic performances of the model with other variables and in different survival time using receiver operating characteristic curve (ROC) with "timeROC" package.Additionally, calibration curve was drawn with "rms" package.

Pan-cancer analysis of three RALUADS genes
We used the online analysis tool TIMER 2.0 (http://timer.cistrome.org/)to systemically evaluate the differential expression level of a given gene in normal and tumor tissues across diverse cancer types in the "Diff Exp" module.The results were automatically created by the website.

Analysis of immune infiltration features
We used CIBERSORT (Chen et al., 2018) and xCell (Aran et al., 2017) to estimate the immune infiltration properties of TME associated with RA and LUAD.CIBERSORT is a professional analytical tool for evaluation of the abundance of different immune cells according to the reference data and gene expression profile.xCell is a novel gene signature-based in silico method for cell immunophenotyping with a total of sixty-four immune cells and three kinds of immune scoring systems.

Analysis of common drug sensitivity
The R package "oncoPredict" is specialized to predict drug sensitivity value according to the gene expression data.By calculating IC50 values, we make linear regression model to indicate the drug response efficacy to chemotherapy corresponding to different RALUADS.In order to examine the relationship between RALUADS and immunotherapy response, immunophenoscore (IPS) based on the expression of MHC molecules, markers of immunomodulators, effector cells and suppressor cells was obtained from The Cancer Immunome Atlas (https://tcia.at/)(Charoentong et al., 2017).Four types of IPS, including IPS, IPS-CTLA4 blocker, IPS-PD-1/PD-L1/PD-L2 blocker, IPS-CTLA4 & PD-1/PD-L1/PD-L2 blocker, were calculated from the TCGA-LUAD cohort.

Cell culture and clinical samples
Human bronchial epithelium cell line (BEAS-2B) and human lung adenocarcinoma cell lines (A549, H1975 and PC9) were purchased from Procell (Procell Life Science&Technology Co., Ltd.).All cells were cultured under standard conditions (37 °C, 5% CO 2 ).All cell lines were authenticated by the short tandem repeat DNA profiling test and checked for absence of mycoplasma contamination.

RNA isolation and quantitative real-time PCR (qRT-PCR)
Total RNA was extracted using TRIzol Reagent (Invitrogen, United States) according to the manufacturer's instructions.NanoDrop was used to detect RNA concentration by A260/ A280 ratio.We performed cDNA synthesis and conducted qRT-PCR with PrimeScript RT reagent kit (EZBioscience, China), and SYBR Green PCR reagent (EZBioscience, China).The reaction was incubated at 95 °C for 10 min followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min.ACTB was used as an internal control.The primer sequences were exhibited in Supplementary Table S4.Data were analyzed using the 2 −ΔΔCT relative quantification method.

Statistical analysis
All bioinformatics analyses were conducted on R software (version 4.3.1)and GraphPad Prism (Version 7.0.0).Independent sample t-tests were employed for normally distributed continuous variables to compare variables between two groups.Continuous variables that were not normally distributed were tested with Wilcoxon test.Correlation analysis between BMRGI and drug IC50 was performed with Spearman method.Statistical significance was defined with p < 0.05.

Unveiling unique features of RA
First, in order to identify genes involved in the onset and development of RA, we examined microarray profiles from 36 RA joint tissues and 7 non-RA joint tissues as control.Principle component analysis (PCA) showed distinct gene expression patterns of RA joint tissues compared with the normal ones (Figure 1A).For differential analysis (|log2FC| > 1, FDR <0.05), we identified 610 differentially expressed genes (DEGs) between RA and normal tissues (Figure 1B).Of note, CYP4B1, DEPP1, HILPDA, GADD45A and other ten genes were significantly upregulated in RA, while SEC61A1 were significantly downregulated.To reveal the functions of DEGs, KEGG analysis showed that DEGs were mainly enriched in "cytokine-cytokine receptor interaction", "chemokine signaling pathway", "NF-kappa B signaling pathway" and etc. Enrichment of "Rheumatoid arthritis pathway" supported that DEGs truly represented the genetic characteristics of RA (Figure 1C).The result of GO analysis showed that, in terms of Biological Process, DEGs were enriched in "muscle system process", "muscle tissue development", "muscle cell differentiation", "muscle organ development", "muscle contraction" and etc. in terms of cellular component, DEGs were enriched in "contractile fiber", "myofibril", "sarcomere", "actin cytoskeleton", "collagen-containing extracellular matrix" and etc. and in terms of molecular function, DEGs were enriched in "actin binding", "receptor ligand activity", "actin filament binding", "glycosaminoglycan binding", "sulfur compound binding" and etc., the above of which were all related to muscle and skeletal function (Figure 1D).We further conducted GSEA analysis and selected some of the most significant changed pathways to make plots.The upregulated pathway "MTORC1 SIGNALING" was reported to be associated with cell growth and further promote lung tumor progression.The downregulated pathways, "apoptosis" and "Glycolysis" might be related with cell proliferation and altered metabolism (Figures 1E-G).

Analysis of shared expression networks of RA and LUAD
In order to explore shared molecular mechanisms of RA and LUAD more deeply, we conducted Weighted Gene Co-Expression Network Analysis (WGCNA) for the two datasets respectively according to the Section 2. For RA, by setting the soft threshold as 8 and dividing the patients into disease and normal groups (Supplementary Figures S1A), we obtained 27 modules of genes significantly related to the characteristics of RA (Figures 3A, B), of which the green module showed the highest correlation with RA (R = 0.78, P-value <0.001).Genes in the green module were highly correlated with gene significance of RA status (R = 0.75, P-value <0.001) (Figure 3C).For LUAD, the soft threshold was set as 7 and patients were divided into tumor and control groups (Supplementary Figures S1B).WGCNA results revealed that there were 22 modules significantly associated with LUAD (Figures 3C,  D), and the brown module was the most significant one (R = 0.54, P-value <0.001).Also, genes in the brown module were highly correlated with gene significance of LUAD (R = 0.71, P-value <0.001) (Figure 3F).

Construction of the common gene signature of RA and LUAD
Based on the findings of differential analysis and WGCNA, we tried to build a signature shared by RA and LUAD.Using Venn tool, we first showed 22 shared DEGs between RA and LUAD (Figure 4A), including ACTC1, ANKRD1, C10orf71, CCN6, CSF3, and FABP4 etc.Then, we performed univariate Cox analysis (Figure 4B) and found that only one gene, CCN6, was the candidate gene to affect overall survival (OS).For genes in two WGCNA modules, Venn plot showed 62 common genes shared by both the characteristics of RA and LUAD (Figure 4C).To select the most significant genes, we first carried out LASSO analysis for dimension reduction (Figure 4D; Supplementary Figures S1C) and found 29 genes.Then, univariate Cox analysis was performed on these genes (Supplementary Figures S1D) and 12 of them were showed to be significantly associated with survival (P-value <0.05).Further, we built multivariate Cox regression model (Figure 4E; Supplementary Table S4) and finally identified two genes, CDCA4 and ERLIN1, to be risk factors affecting OS (CDCA4: HR (95%CI) 1.0104 (1.0004-1.0205),P-value 0.042; ERLIN1: HR (95%CI) 1.0069 (1.0013-1.0125),P-value 0.015).Therefore, we concluded that CCN6, CDCA4, ERLIN1 were three candidate genes with prognostic value in LUAD patients affected by RA at the same time.The three genes all showed significantly higher expression level in tumor tissues compared with the paired normal ones in TCGA dataset (Supplementary Figures S2A-C) and demonstrated strong prognostic value that increased expression is associated with poor overall survival (Supplementary Figures S2D-F).We validated the expression level of the three genes in one normal bronchial epithelium cell line (BEAS-2B) and three human LUAD cell lines (A549, H1975 and PC9).qRT-PCR and Western blot showed that CCN6, CDCA4 and ERLIN1 were significantly upregulated in LUAD cells in terms of mRNA (Figures 4F-H) and protein (Figure 4I) levels.Pathological slides also revealed higher expression of ERLIN1 in LUAD tissues compared with normal lung tissues (Supplementary Figures S2G).Additionally, we conducted pan-cancer analysis and found that these genes had differential expression levels in various tumors compared with the normal ones (Supplementary Figures S2H-J), implying that they might play different roles in different cancers.By combining the expression value and coefficients of the selected genes according to Methods Part, we constructed RA and LUAD prognostic score, which we termed as RALUADS.RALUADS was significantly higher in patients with tumor than in normal people (Figure 4J).According to the median of RALUADS, we divided tumor patients into high RALUADS group and low RALUADS group.The kaplan-meier analysis showed that OS of high RALUADS group was significantly lower than that of low RALUADS group (Figure 4K).Another LUAD cohort from GSE229705 validated the prognostic value of RALUADS.Apart from the higher level of RALUADS in tumor tissues (Figure 4L), we conducted logistic regression and found that higher RALUADS predicted LUAD progression (OR (95%CI) 1.232 (1.0115-1.4528),P-value 0.06) and recurrence (OR (95%CI) 1.3500 (1.0882-1.6115),P-value 0.02).We further explored characteristics of patients in high RALUADS group compared with those in low RALUADS group.After differential expression analysis (|log2FC| > 1), we obtained 5389 upregulated genes and 1658 downregulated genes.KEGG analysis revealed that high RALUADS was positively related to pathways including "Systemic lupus erythematous", "Neutrophil extracellular trap formation", "Spliceosome", "cAMP signaling pathway" and "MicroRNAs in Cancer" (Supplementary Figures S3).Neutrophil extracellular traps have been reported to be a bidirectional mediator of tumor progression (Adrover et al., 2023) and microRNA is a group of small non-coding RNAs which function as sponges and might exert pro-tumor effects (Alahdal and Elkord, 2023).On the other hand, high RALUADS was negatively associated with "Cellular senescence", "p53 signaling pathway" and "Transcriptional misregulation in cancer" (Supplementary Figures S3B).Besides, GSEA analysis demonstrated multiple pathways enriched in high RALUADS group, one of which is KRAS signaling pathway, a crucial "switch" regulating cell proliferation (Supplementary Figures S3C).

Evaluation of the clinical value of RALUADS
In order to apply RALUADS into the clinical settings, we incorporated clinicopathological parameters, including gender, tumor stage, age and smoking amount and performed  5A).Although gender and age were not significant in the COX model, they were commonly considered as crucial in the development and progression in RA and LUADS.Therefore, we combined RALUADS, stage, gender and age to establish a new nomogram model to visualize the contribution of these clinicopathological parameters to OS in different years (Figure 5B).In order to evaluate the accuracy and robustness of the model, we performed ROC curve analysis and the results showed that area under the curve (AUC) of 1-year, 2-year and 3-year were respectively 0.722, 0.710 and 0.711 (Figure 5C).The prognostic performance of various clinical indices in predicting OS of LUAD was also assessed.The ROC curve revealed that compared with other traditional indicator, our nomogram model had the best performance in survival prediction (Figure 5D).In addition, we quantified the performance in terms of calibration, which showed the agreement between predictions and observations (Figure 5E).

Integrated analysis of correlation between RALUADS and tumor immune microenvironment in LUAD
Tumor immune microenvironment (TME) plays a crucial role in shaping the behavior of tumor cells and interaction between tumor and stroma.Hence, we employed two analytical tools, CIBERSORT and xCell to illustrate cellular components of TME respectively.On the one hand, in the analysis of CIBERSORT, compared with low RALUADS group, the levels of M0, M1, activated memory CD4 T cells were significantly higher in high RALUADS group, while naïve B cells, plasma cells, Tregs, resting mast cells, resting dendritic cells (DCs) and activated NK cells were significantly decreased, indicating an activated TME in LUAD (Figure 6A).On the other hand, in the analysis of xCell, the levels of Th1, Th2 and M1 were significantly higher in high RALUADS group compared with the low one, while NKT cells, DCs, CD8 + Tcm cells, M2, plasma cells, CD8 + cells, CD4 + memory T cells, CD4 + Tcm cells and CD4 + naïve T cells were significantly decreased (Figure 6B).The two methods showed similar results in the changes of various immune cell proportions between high and low RALUADS groups.For example, the elevated level of M1 in the high RALUADS group revealed its pro-inflammatory environment and the descended level of plasma cells indicated that decreased antibody production might impair the anti-tumor ability of NK cells via antibody dependent cellular cytotoxicity.Additionally, results of xCell showed decease of a group of memory cells, implying that rapid anti-tumor response might be dampened in the lung tissue (Gebhardt et al., 2023).Then, we analyzed the levels of immune activation related molecules and immune checkpoint molecules in the TME.The radar plot showed that multiple TNF family molecules, IL family molecules and chemokines were significantly positively correlated with RALUADS (Figure 6C).What's more, the expression level of several immune checkpoints including LAG3, HAVCR2 and PDL1 (CD274) were upregulated in high RALUADS group while BTLA, CD40LG etc. were downregulated (Figure 6D).

Discussion
AIDs are a group of diseases that usually involve the whole body with insidious clinical manifestations.RA is one of the most common chronic AIDs that presents as joint damage and extra-auricular lesions, causing heavy global burden, such as movement disability and multi-organ dysfunction.Moreover, there are studies showing the linkage between RA and different types of cancers (Szekanecz et al., 2006;Mao et al., 2020), especially lung adenocarcinoma.In recent years, the prevalence of LUAD has grown rapidly, becoming the main pathological entity in lung cancer and still presenting the low survival rate (Barta et al., 2019;Van Hal and Diab Garcia, 2021).Although there are several investigations revealing the possible mechanisms between RA and LUAD, the most distinct one is immune dysregulation.Given the fact that RA is a systemic inflammatory disease and LUAD has also its unique and heterogenous TME, we sought to explore more deeply into the shared immune mechanism between RA and LUAD.To address this issue, we first conducted differential analysis and WGCNA in RA and LUAD datasets respectively.Enrichment analysis showed that DEGs in RA participated in several tumor related pathways, such as metabolism and proliferation, while DEGs in LUAD participated in various immune processes, indicating the interaction between the microenvironment of RA and LUAD.WGCNA identified two modules significantly related to RA and LUAD.Furthermore, we subjected the overlapping genes of differential analysis and WGCNA to LASSO analysis and COX regression analysis to screen out genes that most represent prognostic value.Finally, three genes were identified as shared genes of RA and LUAD, CCN6, CDCA4 and ERLIN1, whose higher expression level in LUAD cells compared with normal epithelial cells were confirmed by qRT-PCR, WB and IHC.CCN6, CDCA4 and ERLIN1 were used to establish RALUADS.The score was significantly higher in tumor patients and correlated with poor prognosis in both the discovery cohort and validation cohort.
CCN6, or WISP3, or Wnt1-inducible signaling pathway protein 3, belongs to CCN family and is structurally characterized by a glycosylated part and four conserved cysteine-rich domains (Holbourn et al., 2008).Members of CCN family, including CCN1-6, are evolutionarily conserved and participate in regulating different pathophysiological processes, including cell proliferation, adhesion, angiogenesis, ECM modeling, migration, tumor growth (Jun and Lau, 2011).Studies have highlighted that expression of CCN6 is higher in RA synovium and fibroblast-like synoviocytes compared with osteoarthritis and normal synovial tissue (Cheon et al., 2004).What's more, a recent study found that gender related lncRNA XIST could bind to GATA1, leading to CCN6 upregulation and driving RA pathogenesis by promoting SF proliferation and angiogenic activity (Yu et al., 2023).Although there is still no evidence that relates CCN6 to LUAD, CCN6 have been shown to play a crucial role in gastrointestinal cancers (Thorstensen et al., 2003), breast cancers (Djomehri et al., 2020) and chondrosarcoma metastasis to lung (Tzeng et al., 2018).CDCA4, Cell Division Cycle-Associated Protein 4, is a member of E2F family of transcription factors having a role cell cycle regulation.Pan cancer analysis revealed that a great number of tumors highly expressed CDCA4, which was associated with poor survival and different immune infiltration characteristics (Fang et al., 2022).In particular, CDCA4 was found to be a marker for poor OS in patients with LUAD (Tan et al., 2022).Although there are not evidences showing direct association between CDCA4 and RA, previous studies have shown that transcription factors in E2F family were enriched in RA (Takeshita et al., 2019).ERLIN1, also named for ER Lipid Raft Associated 1, is a part of a protein complex which mediates degradation of inositol 1,4,5trisphosphate receptors in the endoplasmic reticulum, important for cellular cholesterol homeostasis.However, ERLIN1 functions beyond lipids metabolism.In breast cancers, ERLIN1 was discovered to be targeted by estrogen/MYC/miR-26 axis and promote cell growth (Tan et al., 2014).Moreover, in pancreatic adenocarcinoma, higher expression of ERLIN1 was correlated with poor survival and lower CD8+T cell infiltration (Chen et al., 2022).Taken together, the three genes in RALUADS were closely related to poor cancer survival.As some of them have not been reported to participate in the pathogenesis of RA or LUAD, we are the first to reveal the possible novel roles of these genes in LUAD with RA.
Based on RALUADS, we constructed a nomogram model with three other clinicopathological parameters and demonstrated that it had accurate and robust performance in predicting OS of LUAD.Furthermore, we investigated TME affected by RALUADS and found that immune infiltration was characterized by increased level of M1 and decreased levels of plasma cells in high RALUADS group.Also, immune chemotaxis and stimulation were significant in high RALUADS group while various immune checkpoints were positively or negatively correlated with RALUADS.Given the abundant and complex immune components in LUAD TME, we found that RALUADS could be a potential biomarker to predict immune infiltration and immunotherapy response.IPS tended to be higher in low RALUADS group and consistently, patients with lower RALUADS would be expected to benefit more from anti-CTLA4 and anti-PD1/PDL1/PDL2 immunotherapy.In contrast, by examining the differences in the therapeutic effects of chemotherapy and targeted medications in LUAD patients, IC50 of anticancer drugs was lower in the higher RALUADS group, including Cisplatin, Docetaxel and Gefitinib etc., indicating that patients with higher RALUADS might be more sensitive to chemotherapy and targeted drugs.Previous studies have put forward various models regarding the features and prognosis of LUAD.As extracellular matrix plays a crucial role in tumor progression and invasion, Wang et al. (2021) constructed a score with three integrin genes and metastasis-related microenvironmental pathways were enriched in high score group.What's more, Zhang et al. (2023) established a signature with seventeen basement membrane related genes and patients with high risk had higher tumor mutation burden, lower immune score and poorer prognosis, which was consistent with our results that high RALUADS was associated with low IPS.Another study revealed the relationship of immune cell death and TME of LUAD and defined nine damage-associated molecular pattern related genes to be related to active immune activity and response to immunotherapy (Wu et al., 2023).However, a majority of these researches, focusing on LUAD, explored only one type of molecules, which might be inadequate to depict the panorama of LUAD.In comparison, we investigated the complex environment of LUAD affected by RA and screened out three genes from different gene family with diverse functions via differential gene analysis, WGCNA and survival analysis.This contributed to demonstrate the comprehensive profile of LUAD, including immune characteristics, cell death and cell metabolism etc.
In our research, we employed comprehensive analysis tools to identify the genes with the highest relationship in RA and LUAD.Three genes, CCN6, CDCA4 and ERLIN1, exhibited significant correlation and constituted RALUADS.Notably, patients with higher RALUADS showed poor prognosis, suggesting that the signature could be a promising biomarker in mortality prediction of NSCLC in MS patients.Our investigation further demonstrated the regulatory effects of RALUADS on immune infiltration and expression of immune activation molecules and immune checkpoints.In additional, RALUADS was significantly related to the effects and sensitivity of multiple therapeutic drugs.Our findings suggested that RALUADS could be applied in the clinical settings to help stratify LUAD patients and further guide the treatment strategies (Figure 8).
Our study has some limitations.First, the expression of CCN6, CDCA4 and ERLIN1 have not been validated by in vivo experiments.Additionally, multi-omics data may be needed to   Atlas (TCGA, https://portal.gdc.cancer.gov/)database.Dataset of validated cohort can be accessed through GSE229705.The detailed statistical data in the study are included in the article/ Supplementary Material.

FIGURE 1
FIGURE 1 Unveiling unique features of RA. (A) PCA of 36 RA and 7 normal samples in GSE236924 datasets (B) Volcano plot of DEGs between RA and NC.(C) KEGG enrichment analysis of DEGs in RA compared with NC. (D) GO enrichment analysis of DEGs in RA compared with NC. (E-G) GSEA analysis of RA compared with NC.RA, rheumatoid arthritis; PCA, principal component analysis; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology.Note: NC, normal tissues.

FIGURE 2
FIGURE 2 Identification of genetic patterns of LUAD.(A) PCA of 549 LUAD and 41 normal samples in TCGA-LUAD datasets (B) Volcano plot of DEGs between LUAD and NC.(C) KEGG enrichment analysis of DEGs in LUAD compared with NC. (D) GO enrichment analysis of DEGs in LUAD compared with NC. (E-G) GSEA analysis of LUAD compared with NC.LUAD, lung adenocarcinoma; PCA, principal component analysis; DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology.Note: NC, normal tissues.

FIGURE 3
FIGURE 3 Shared expression networks of RA and LUAD.(A,D) Gene clustering dendrogram.(B,E) Correlation heatmap between modules and clinical traits.(C,F) Linear correlation between module membership and gene significance.

FIGURE 4
FIGURE 4 Construction of the common gene signature of RA and LUAD.(A) Venn plot of shared DEGs between RA and LUAD.(B) Forest plot showing the univariate Cox regression analysis of shared DEGs.(C) Venn plot of shared significant module genes between RA and LUAD.(D) LASSO coefficient profiles.(E) Forest plot showing the multivariate Cox regression analysis of selected shared significant module genes.(F-H) mRNA expression of CCN6, CDCA4, ERLIN1 in different cell lines by qRT-PCR.(I) Protein expression of CCN6, CDCA4, ERLIN1 in different cell lines by Western blot.(J) Expression level of RALUADS in LUAD and normal tissues in TCGA-LUAD cohort.(K) Kaplan-Meier survival curves stratified by RALUADS.(L) Expression level of RALUADS in LUAD and normal tissues in validation cohort.

FIGURE 5
FIGURE 5 Evaluation of the clinical value of RALUADS.(A) Forest plot showing multivariate Cox regression analysis of common clinicopathological parameters and RALUADS.(B) Nomogram integrating significant clinicopathological parameters and RALUADS.(C) ROC curves of nomograms predicting 1-year, 2year and 3-year survival rates.(D) ROC curves for clinicopathological parameters, RALUADS and nomogram.(E) Calibration curves for nomograms predicting 1-year, 2-year and 3-year survival rates.

FIGURE 6
FIGURE 6 Integrated analysis of correlation between RALUADS and tumor immune microenvironment in LUAD.(A,B) Differences in infiltration level of various types of immune cells between high and low RALUADS groups by analysis of CIBERSORT and xCell.(C) Correlation between RALUADS and immune activation related molecules.(D) Correlation between RALUADS and immune checkpoint molecules.Note: Treg, regulatory T cell; Th1, type 1 helper T cell; DC, dendritic cell; Tcm, central memory T cell; Tem, effector memory T cell.

FIGURE 7
FIGURE 7 Identification of sensitive drugs in LUAD patients with RA. (A-D) Correlation analysis between RALUADS and response to immunotherapy.(E-L) Correlation analysis between IC50 value of chemotherapy or targeted drugs and RALUADS.

FIGURE 8
FIGURE 8Schematic diagram of the shared signature of RA and LUAD.