Gender specific eRNA TBX5-AS1 as the immunological biomarker for male patients with lung squamous cell carcinoma in pan-cancer screening

Data preparation and patient samples

The RNA expression file and related clinical information on carcinomas were downloaded from the TCGA database via University of California, Santa Cruz (UCSC) Xena (https://xenabrowser.net/). In total, 34 types of solid cancers and the Genotype-Tissue Expression (GTEx) project were included for supplementing normal samples of human organs. The list of cancer types is provided in the Supplemental Files (Table S1). Besides, 364 patients with NSCLC from GSE37745 and GSE50881 were also involved in the study. GSE37745 and GSE50881 were combined via batch analysis. Predicting Specific Tissue Interactions of Genes and Enhancers (PreSTIGE) predicted the information on eRNAs and their potential targets (Corradin et al., 2014; Vučićević et al., 2015).

The validation cohort consisted of 30 male patients with squamous cell carcinoma (SCC) and 64 patients with adenocarcinoma (ADC) from the Department of Thoracic Surgery, Shandong Provincial Hospital. The ethics review committee of Shandong Provincial Hospital approved this study (SWYX: No.2021-260). It was performed based on the Declaration of Helsinki and Good Clinical Practice guidelines, as defined by the International Conference on Harmonization.

Screening of eRNA with prognostic significance and functional enrichment analysis

The Kaplan–Meier (KM) survival curves were used to pick valuable eRNAs in specific cancer types. Pearson correlation coefficients were used to evaluate the correlation between eRNAs and protein-coding genes (PCGs). PCGs that met the standards were considered lncRNA-related PCGs (p < 0.001; |Pearson correlation coefficient| > 0.4). The screened PCGs were used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on R packages colorspace, stringi, ggplot2, enrichplot, clusterProfiler, and DOSE. The proteomics data were obtained from The Cancer Proteome Atlas (TCPA) database to explore the potential downstream targets.

Assessment of infiltrating immune cells

A computational framework CIBERSORT was used to assess the infiltration degree of immune cells. A total of 22 types of immune cells were evaluated: naïve/memory B cells, plasma cells, CD8+ T cells, naïve CD4+ T cells, resting/activated memory CD4+ T cells, follicular helper T cells, regulatory T cells (Treg cells), gamma delta T cells, resting/activated NK cells, monocytes, M0/M1/M2 macrophages, resting/activated dendritic cells, resting/activated mast cells, eosinophils, and neutrophils.

Quantitative real-time polymerase chain reaction assay

Total RNA was extracted using TRIzol (Lot A2A0209, Accurate Biotechnology, Hunan Province, China) reagent and assessed using a Nanodrop 2000 system (Thermo Fisher Scientific, Waltham, MA, USA). A reverse transcription kit (A2A1386) was obtained from Accurate Biotechnology (Human) Co., China. The sequence of 18S, TBX5, and TBX5-AS1 primers is listed in Table S2. The real-time polymerase chain reaction (PCR) was performed on a LightCycler 480 II (Roche, Basel, Swiss), using the SYBR Green system (Lot A2A1436, Accurate Biotechnology, Hunan Province, China).

Prediction of response to ICI treatment

Immunophenoscore (IPS) was applied to The Cancer Immunome Atlas (TCIA, https://tcia.at/) to quantify tumor immunogenicity and predict patients’ response to ICI treatment (Charoentong et al., 2017). The T-cell immune dysfunction and exclusion (TIDE) was used to predict whether patients could benefit from ICIs based on the cytotoxic T lymphocyte (CTL) status (Fu et al., 2020b). The preparation files were obtained from UCSC and normalized.

The construction of competing endogenous RNA (ceRNA) network and prediction of combination region of TBX5-AS1 to AR

We used Starbase (https://starbase.sysu.edu.cn/) to predict potential target miRNA of TBX5-AS1 and miRDB (http://mirdb.org/) for AR. Besides, the catRAPID (http://service.tartaglialab.com/page/catrapid_group) was performed to explore the possibility of TBX5-AS1 combining AR.

Statistical analysis

All R packages were run under the environment of R 3.6.1. GraphPad Prism (version 8.0) and SPSS 24.0 (IL, USA) were used to perform statistical analyses. A p value < 0.05 indicated a statistically significant difference.

Contradictory effects of TBX5-AS1 in pan-cancer survival analyses

Screening of eRNAs with prognostic significance was initially performed in an NSCLC cohort (Tables 1 and 2). The Venn diagram showed that TBX5-AS1 might play a key role in the prognosis of patients with ADC and SCC (Fig. 1A), despite its equally strong correlation to target TBX5 in both NSCLC subtypes (R > 0.9 and p < 0.0001). Interestingly, TBX5-AS1 had opposite effects on the prognosis (Figs. 1B and 1C). TBX5-AS1 seemed to act as a risk factor for patients with SCC and became a protective factor in patients with ADC.

For further pan-cancer analysis, the prognostic values of TBX5 and TBX5-AS1 were analyzed in 31 other solid cancers from the TCGA database via hazard ratio (HR) calculation and KM p value. The data were scattered in four quadrants, split by p = 0.05 and HR = 1. In ACC and uterine corpus endometrial carcinoma (UCEC) cohorts, the expression level of TBX5-AS1 was associated with worse outcomes. The same was true for TBX5 in kidney renal clear cell carcinoma, UCEC, and bladder urothelial carcinoma (Fig. 1D). On the contrary, patients with skin cutaneous melanoma having high expression of TBX5 might possess a better prognosis.

Specific connection between TBX5/TBX5-AS1 and gender

In the LUAD cohort, the expression levels of TBX5 and TBX5-AS1 were influenced by gender (Fig. 2A). Similarly, the expression levels of TBX5 and TBX5-AS1 in male patients with SCC were significantly lower than those in female patients (Fig. 2B). However, there was no difference of TBX5-AS1 and TBX5 expression based on age grouping (Fig. S1). The analysis of organ gene expression based on the GTEx database showed that TBX5 and its eRNA were mainly enriched in the heart and lungs of the normal population (Fig. S2). After distinguishing gender as a variable, TBX5-AS1 and TBX5 were expressed differently in the lung of male and female patients (Figs. 2C & 2D). Different distributions of TBX5-AS1 between male and female patients occurred in the skin, brain, and adipose tissue. A similar distribution of TBX5 also existed in the skin, adipose tissue, breast, and adrenal gland. That indicated eRNAs with targets might have tissue and gender specificity (Fig. S2). eRNAs might played different roles in male and female. Besides, the pan-cancer screen predicted that this expression difference between female and male did not appear in other cancer cohorts of TCGA database (Figs. 2E & 2F).

Survival analysis of male and female patients in NSCLC and cohort validation

Given that the gender could cause TBX5 and TBX5-AS1 lung research bias, the KM curves were plotted in male and female patients separately. Interestingly, the significant results of TBX5 and TBX5-AS1 converged at the male LUSC cohort (p = 0.006 for TBX5-AS1 and 0.040 for TBX5; Fig. 3A). The prognosis of the rest of the groups, female LUSC, male LUAD, and female LUAD, were not significantly associated with TBX5 or TBX5-AS1 (Figs. 3B–3D).

Moreover, 32 female patients with ADC, 32 male patients with ADC, and 30 male patients with SCC were included to form a validation cohort. The survival analysis proved no significant difference between low- and high-expression groups of TBX5-AS1 in the female and male ADC cohorts (Fig. 3E). The prognosis was better in the low-expression group of TBX5-AS1 compared with the high-expression group in male patients with SCC, and the quantitative reverse transcriptase (qRT)-PCR verified that a strong correlation existed between eRNA and its homologous target once again (r = 0.9576 and p < 0.0001; Fig. 3F).

Furthermore, the gender specificity of TBX5 and TBX5-AS1 in patients with NSCLC was validated by the combined cohorts based on GSE37745 and GSE50881. The TBX5-AS1 (p = 0.024) and TBX5 (p = 0.014) exerted a prognostic biomarker effect only in male patients with SCC (Fig. 4A). For female patients with NSCLC and those with SCC, the result was not statistically significant (Figs. 4B–4D).

Correlation between TBX5-AS1 and AR signaling pathway

Considering the vast impact of gender differences on the role of TBX5-AS1, a correlation analysis between the expression levels of TBX5-AS1 eRNA and AR mRNA was performed in LUSC and LUAD samples. For male patients, the correlation was evident: R = 0.55 in LUSC and 0.52 in LUAD (both p < 0.0001; Fig. 4E). Furthermore, the GSEA results indicated that the androgen-mediated pathway (GO_ANDROGEN_RECEPTOR_SIGNALING_PATHWAY) was enriched in the male LUSC (p = 0.002) and male LUAD (p = 0.021; Fig. 4F) groups. For female patients, the outcomes were not significant (p = 0.079 for LUAD and p = 0.088 for LUSC). Based on miRDB and Starbase databases, we performed Venn diagram and found hsa-miR-873-5p might become the important bridge to connect ceRNA network of TBX5-AS1 and AR (Fig. 4G). Meanwhile catRAPID was used to predict the combination region between TBX5-AS1 and AR, under the consideration of classic interaction of eRNAs and hormone receptors (Table S3). Thus, the schematic diagram was showed in Fig. 4H to exhibited that TBX5-AS1 could stabilized the expression of AR via ceRNA network and interact with AR to influence downstream pathways.

Predicting prognosis in male and female pan-cancer cohorts

Figure 1D shows the prognostic significance of TBX5-AS1 in ACC and UCEC cohorts, and the consistency in predicting prognosis between TBX5-AS1 and TBX5 existed only in female UCEC samples. Thus, given the connection between TBX5-AS1 and gender, the pan-cancer survival analysis was performed after gender grouping. In various cancers, the prognostic significance of TBX5-AS1 and its homologous target became more consistent, which were the risk factors for liver hepatocellular carcinoma (LIHC) in male patients and glioblastoma multiforme (GBM) and UCEC in female patients (Figs. 5A and 5B). Additionally, a significant prediction was explored for both TBX5-AS1 and TBX5.

Functional analyses of TBX5-AS1 in male LUSC samples

The KEGG and GO analyses were performed based on the correlated genes of TBX5-AS1 (R > 0.4). The immune-related terms, “leukocyte migration,” “adaptive immune response based on somatic recombination,” and “humoral immune response mediated by circulating immunoglobin,” were significant in GO analyses (Fig. 5C). Similarly, the result of KEGG analysis also indicated that the impact of TBX5-AS1 on cell biology might focus on immunological pathways, for example, “Th17 cell differentiation,” “leukocyte transendothelial migration,” and “Th1 and Th2 cell differentiation” (Fig. 5C).

Target prediction at the protein level based on TCPA databases indicated that CKIT, PTEN, and PAXILLIN were positively related to TBX5-AS1 and TBX5, and the correlated relationship to the targets ACC1, G6PD, and KEAP1 was negative (Fig. 5D). Besides, the potential targets individually to eRNA were listed as follows and (−)/(+) meant negative/positive correlation: CLAUDIN7 (−), ACC_pS79 (−), CD31 (+), MYOSINIA_pS1943 (+), TFRC (−), ASNS (−), and MAPK_pT202Y204 (+).

Evaluation of infiltrated immune cells in male patients with LUSC

The infiltrated fraction of plasma cell and resting memory CD4 T cells increased in the high-expression group of TBX5 (Figs. 5E and 5G). Meanwhile, the situation in follicular helper T cells was the opposite (Fig. 5F). For TBX5-AS1, the connection to resting memory CD4 T cells and follicular helper T cells was similar to that for TBX5 (Figs. 5H and 5I).

The correlated heatmap exhibited that the infiltrated fraction of resting memory CD4 T cells was connected with the immunological infiltration of CD8 T cells (−0.56), activated CD4 T cells (−0.37), and follicular helper T cells (−0.36; Fig. S3), which supported the aforementioned results.

Prediction of response to ICI treatment according to the expression of TBX5-AS1

Based on the consideration that TBX5-AS1 influenced the T lymphocyte function, its role in response to lung squamous cancer immunotherapy was further explored (Waldman, Fritz & Lenardo, 2020). The IPS files of SCC patients were downloaded from the TCIA database, and the high value suggested that patients were sensitive to immunotherapy. The IPS of programmed cell death protein 1 (PD1) inhibitors was higher in the high-expression group of TBX5-AS1 (mean = 7.146) than in the low-expression group (mean = 6.704; Fig. 6A). The IPS of the Cytotoxic Lymphocyte Antigen 4 (CTLA-4) blocker was higher in the high-expression group of TBX5-AS1 (mean = 7.354) than in the low-expression group (mean = 7.022; Fig. 6B). As for the combination therapy of CTLA-4 blockers and PD1 inhibitors, the IPS was higher in the high-expression group of TBX5-AS1 (mean = 6.917) than in the low-expression group (mean = 6.404; Fig. 6C). All results were statistically significant.

The TIDE algorithm focused on the malignant degree of the tumor microenvironment (TME) (Fu et al., 2020b). When infiltrated CTLs were rare, the TIDE score was equal to the T-cell exclusion score and also to T-cell dysfunction in the high-CTL situation. Figure 6D shows that patients with high expression of TBX5-AS1 might possess worse TME, which reflected in the T-cell dysfunction (mean of high expression = 1.735 and mean of low expression = −0.331; p < 0.0001).

Because TIDE algorithm was developed for lung cancer and melanoma, we only used IPS to predict the immunotherapy responses in male patients with LIHC and female patients with GBM and UCEC. In addition to meaningless in GBM cohort, the outcome in UCEC and LIHC cohorts was significant but opponent to result in SCC patients (Table S4). In this situation, patients with low TBX5-AS1 expression might are sensitive to immunotherapy including PD1 inhibitors, CTLA-4 blockers and combination therapy of both them.