A Pan-Cancer Landscape of HOX-Related lncRNAs and Their Association With Prognosis and Tumor Microenvironment

The highly conserved homology cassette family (HOX) as well as 18 referenced long non-coding antisense transcripts (HOXATs) play vital roles in the development of some cancers. Nevertheless, their expression patterns as well as their association with cancer prognosis and the tumor microenvironment (TME) in pan-cancers are still unclear. Here, based on public databases, the expression levels of HOXATs, their prognostic potentials, and correlation with tumor mutation burden (TMB), immune cell infiltration, immune subtype, immune response-related genes, and stemness scores corresponding to 33 tumor types were analyzed systematically using R language. The results of the analysis indicated that different cancer tissues show different HOXAT expression profiles. Further, HOXAT expression showed association with cancer prognosis and immune and stemness regulation. Gene set enrichment analysis also demonstrated that HOXATs participate in cancer- and immune-related pathways, and based on their expression levels, HOTAIRM1 and HOXB-AS1 showed potential involvement in oncogenesis as well as possible involvement in immune regulation across a variety of cancer types. Further investigation also confirmed a significantly higher expression of HOXB-AS1 in GBM than in lower grade glioma tissues. Importantly, in vitro cell function experiments indicated that HOXB-AS1 supports cancer stem cell and plays a fundamental role in glioma metastasis. In conclusion, our results provide valuable resources that can guide the investigation of the mechanisms related to the role of HOXATs in cancers as well as therapeutic analysis in this regard.


INTRODUCTION
For Long non-coding RNAs (lncRNAs), as one of the key regulators of cancer, function not only as important molecular markers for cancer prognosis, but also as markers of molecular characteristics at the pan-cancer level. Increasing studies have shown that several lncRNAs affect cancer tumorigenesis and contribute to cancer pathogenesis (Akbari et al., 2020;Bi et al., 2020).
The tumor microenvironment (TME), which contains multiple types of stromal cells, tumor cells, cancer stem cells (CSCs), tumor-associated immune cells and others, plays a major role in tumor growth and metastasis, and particularly, the immune microenvironment constitutes a vital element of tumor biology. Immune system dysregulation can be a major cause of cancer development (Chen et al., 2019a;Carlsten and Järås, 2019). The most recent studies suggest that CSCs represent another important component of the TME that have differentiation abilities to generate the original lineage cells for regulating tumor occurrence, development, recurrence, metastasis and drug resistance. It has also been suggested that the activities of CSCs and immune cells represent the potential mechanisms that link the TME and cancer. Recently, increasing evidence has indicated that lncRNAs play fundamental roles in the regulation of the immune microenvironment and are closely associated with the development and progression of various CSCs. For example, the lncRNA, NeST has been shown to be related to T-cell activation, and is critical for immune response regulation (Chen et al., 2019b), while the lncRNA, NRON has been shown to maintain a resting state of T cells by sequestering phosphorylated NFAT in the cytoplasm . Moreover, it has also been demonstrated that the lncRNA, lnc-sox5, modulates CRC tumorigenesis by unbalancing TME (Chen et al., 2016), while the lncRNA, LINK-A functions as an oncogene that promotes tumor immune evasion by sensitizing T cells (Dong et al., 2018). Additionally, it has been observed that the depletion of oncogenic HotairM1 promotes the self-renewal of CSCs via the HOXA1-Nanog signaling loop. It has also been observed that lncTCF7 is required for liver CSC self-renewal and tumor propagation (Chi et al., 2019). However, further studies on lncRNAs and their roles in the TME of different cancer types still require investigation.
Transcription factors are encoded by homeobox (HOX) genes, which are assigned to a group of 13 orthologs and paralogs in four clusters and share similar positions in the homozygous frames, with 9-11 genes per cluster (Cruz et al., 2017). Multiple reports have shown that the dysregulation of 39 HOX genes is implicated in various aspects of cancer. Similar to protein-coding transcripts, accumulating experimental evidence indicates that lncRNA transcripts in the Hox gene locus, with 18 referenced non-coding antisense transcripts (HOXATs) are implicated in cancer progression. These 18 HOXATs are located in four HOXATs clusters, with the HOXA region consisting of HOTAIRM1, HOXA-AS2, HOXA-AS3, HOXA10-AS, HOXA11-AS, and HOTTIP, the HOXB region consisting of HOXB-AS1, HOXB-AS2, HOXB-AS3, HOXB-AS4, and PRAC2, the HOXC region consisting of HOXC13-AS, HOXC-AS1, HOXC-AS2, HOXC-AS3, and HOTAIR, and the HOXD region consisting of HAGLR and HOXD-AS2. For example, it has been reported that the expression levels of HOXA10-AS, HOXA11-AS, and HOXA-AS are significantly upregulated in lung adenocarcinoma (LUAD), and their high expression levels in such a cancer type are associated with shorter survival periods in patients (Hao and Yu, 2018;Dasti et al., 2020). Further, it has also been observed that the expression levels of HOTAIRM1, HOXA10-AS, HOXA10-AS, HOXA11-AS, HOXA-AS2, HOXA-AS3, HOXB-AS1, HOXC13-AS, HOXC-AS3, and HOXD-AS2 are upregulated in glioma, and are positively correlated with malignancy Li et al., 2020a), and importantly, some studies have shown that HOXATs can play fundamental roles in immune system regulation (Li et al., 2020b;Liang et al., 2020). However, the different expression levels and the specific functions of the individual HOXATs in pan-cancer as well as their association with the TME still remain unclear. Therefore, an in-depth investigation and analysis of their potential roles in pan-cancer is warranted.
Hence, in this study, we attempted to explore and analyze the expression profile and prognostic significance of HOXATs across multiple cancer types and investigated the potential role of the HOXATs in the modulation of the TME as well as stemness score. Additionally, we also explored and verified the potential biological function and pathways of the HOXATs with respect to the different aspects of cancer.

Datasets Used
In this study, datasets from the public TCGA pan-cancer data portal (https://xenabrowser.net/datapages/), containing 10,635 samples spanning 33 cancer types with clinically annotated multi-omic data (clinical data, stemness scores based on mRNA (RNAss) and DNA methylation (DNAss), and immune subtypes) were used for analyses.

Correlations Between HOXATs and Tumor Microenvironment
ESTIMATE-based immune and stromal scores were used to analyze the infiltration levels of immune and stromal cells in the different tumor types. Six immune subtypes were defined to measure the levels of immune infiltrate types in the TME based on the TCGA pan-cancer database. Further, tumor stemness features extracted from the transcriptomics and epigenetics of TCGA tumor samples were used to measure the stem-cell-like features of the tumor cells. CIBERSORT was then used for the calculation of the TMB scores of the tumor cells and the correlation between the TMB and HOXATs levels was analyzed . Further, the Spearman or Pearson correlation coefficients obtained were then used to evaluate the relationships between the expression levels of the HOXATs and RNAss, DNAss, immune score, estimate score, TMB, immune infiltration, and the expression of immune checkpoint genes as well as tumor stemness markers.

Cells Culture
Two human glioma cell line samples (U251 and HS683) were obtained from Cell Culture Center, Chinese Academy of Medical Sciences (Shanghai, China). The cells were maintained in a medium containing Dulbecco's Modified Eagles Medium (DMEM), 10% fetal bovine serum (Invitrogen, Carlsbad, CA, United States), and 1% penicillin-streptomycin solution. Thereafter, they were cultured under standard cell culture conditions consisting of 5% CO 2 and 95% relative humidity at 37 C.
lncRNA Knockdown U251 and HS683 cells in the logarithmic growth phase were incoculated into 6-well plates, after the cells grew to 50-60% fusion. The lncRNA-HOXB-AS1 smart silencer (a mixture containing three siRNAs and three antisense oligonucleotides) and its negative control, lncRNA smart silencer, NC, synthesized by RiboBio were used (Guangzhou, China). The cells were transfected using Lipofectamine RNAiMAX (Invitrogen, Carlsbad, CA, United States) following the manufacturer's protocol.

Sphere Formation Assay
Single human glioma cell line suspensions were washed twice with serum-free phosphate buffer solution (PBS) and transferred into a serum-free stem cell complete medium (containing 20 μg/L bFGF and 20 μg/L EGF) followed by the addition of B27. Thereafter, the cells were cultured at 37 C in 5% CO 2 , and sphere formation was analyzed.

Invasion Assay and Wound Healing Assay
To realize invasion assay, transwell chambers (Corning, Tewksbury, MA, United States) were coated with Matrigel (dilution 1:8; BD Biosciences, Bedford, MA, United States). Thereafter, a 120-μL serum-free medium containing 5 × 10 4 cells was added to the upper chamber, while a 600-μL complete medium containing 10% FBS was added to the lower chamber. After transfaction, the cells on the lower surface were fixed with methanol and stained with 0.1% crystal violet and then captured using an inverted microscope (AMG, United States). To realize the wound healing assay, U251 and HS683 cells were incoculated into 6-well plates. When the cell cultures reached 80% confluence, the cells were scraped using a 10-μL micropipette tip. Thereafter, images were captured using an inverted microscope with a camera (AMG, United States). The distance travelled by the migrating cells was quantified to assess cell migration capacity using ImageJ software version 1.46 (National Institutes of Health, Bethesda, MD, United States).

MTT Assay
Cell viability was monitored using an MTT cell viability assay kit (Sigma-Aldrich, St. Louis, MO, United States). The transfected cells were seeded in 96-well plates. Thereafter, 10 μL MTT (5 mg/ml) was added to each well and incubation was continued for 2 h, followed by the dissolution of the precipitate formed in 100 μL of dimethyl sulfoxide (DMSO). Absorbance measurements were then performed at 570 nm under a microplate spectrophotometer.

Bioinformatics and Statistical Analyses
All computational and statistical analyses were performed using R software v4.0.3 and Prism 8 GraphPad (https:// www.graphpad.com/scientific-software/prism/). Correlations between stemness markers and HOXAT expression levels were analyzed using the GEPIA online dataset (http://gepia. cancer-pku.cn/), which includes the gene expression profiles of different cancers and normal tissues in the TCGA database as well as genotype-tissue expression (GTEx). p-values < 0.05 were considered statistically significant. Further, the results were adjusted for multiple comparisons.

Brief Literature Review
We analyzed the connection between the expression of HOXATs and prognosis in patients with cancer by exploring all published related literature. A total of 1,802 potential studies were identified, and after title and abstract screening, 65 studies were found to be eligible for inclusion. Next, we examined the full text of these 65 articles, in which most of the relationship between HOXATs and their related genes, as well as the mechanisms of regulation were investigated. Conversely, we observed that studies in which the relationship between the HOXATs and the immune microenvironment in pan-cancers were investigated were scarce. Further, no study aimed at investigating the immune response of HOXATs in tumors was identified ( Table 1).

HOXATs Expression and Correlation With Pan-Cancers
To clarify the intrinsic expression pattern of HOXATs, we described the differential expression of HOXATs considering the 33 available cancer types in the TCGA pan-cancer database. Notably, most of the HOXATs showed a striking intra-and inter-tumor heterogeneity in these different cancer types. Among them, HOTAIRM1 and HAGLR showed relatively higher expression levels in all the cancer types than the other HOXATs ( Figure 1A; Supplementary Figures S1-S18). Further, we systematically investigated the expression levels of all 18 genes in primary patient tumor samples corresponding 16 cancer types. This investigation also included at least five normal samples. For instance, the expression levels of HOTAIRM1 showed significant inter-tumor heterogeneity in esophageal carcinoma (ESCA), GBM, and kidney papillary cell carcinoma (KIRP), which are associated with very high HOTAIRM1 expression levels. Conversely, breast cancer (BRCA), colon adenocarcinoma (COAD), kidney  chromophobe (KICH), and rectum adenocarcinoma (READ) were found to be characterized by lower HOTAIRM1 expression levels (p < 0.05; Figure 1B). Therefore, these results indicated the existence of intrinsic differences between the different tumor types with respect to the expression of HOXATs. Next, we explored the potential relationship between the HOXATs. In this regard, Pearson correlation analysis revealed the following significant correlations between the HOXATs: HOTAIRM1 with HOXA-AS2 and HOXA-AS3; HOXA-AS2 with HOXA-AS3 and HOXA10-AS; HOXA-AS3 with HOXA10-AS and HOXB-AS2; HOXA10-AS with HOXA11-AS and HOTTIP; HOXA11-AS with PRAC2; HOTTIP with HOXB-AS4 and PRAC2; HOXB-AS1 with HOXB-AS2 and HOXC-AS1; HOXB-AS2 with HOXB-AS3; HOXB-AS3 with HOXB-AS4; HOXC13-AS with HOTAIR; HOXC-AS1 with HOXC-AS2, HOXC-AS3, and HOTAIR; HOXC-AS2 with HOXC-AS3 and HOTAIR; HOXC-AS3 with HOTAIR; and HAGLR with HOXD-AS2 (R > 0.4; Figure 1C). These significant correlations indicated that several HOXATs interact with each other during the pathogenesis of pan-cancer. Interestingly, HOTAIRM1 and HOXA-AS2 showed positive correlation with the other 16 HOXATs, suggesting that they possibly exert a regulatory effect on all the other 16 HOXATs and may be critical factors in cancer development. Notably, the correlation coefficients between HOXA11-AS and HOTTIP, HOXC-AS1 and HOXC-AS2 as well as that between HOXC-AS2 and HOXC-AS3 were all ≥ 0.7, suggesting that they possibly share common features and functions.

Prognostic Analysis of HOXATs in Pan-Cancer
The prognostic significance of HOXATs in the different cancer types was analyzed using Cox regression (Figures 2A,B). In cancers, including adenoid cystic carcinoma (ACC), KIRC, LGG, and thymoma (THYM), approximately half or more than half of the HOXATs showed association with worse survival of patients. Additionally, it was also observed that certain HOXATs possibly exerted obviously different prognostic effects across the various cancer types. For instance, HOXA11-AS predicted poor prognosis in patients with ACC, KICH, KIRC, LGG, and lung adenocarcinoma (LUAD), HOTTIP predicted poor prognosis in patients with ACC, KICH, KIRC, LGG, LUAD, and uveal melanoma (UVM), and HOTAIRM1 predicted poor prognosis in patients with ACC, GBM, KICH, KIRC, LGG, THYM, and UVM patients, but was found to be associated with a survival advantage for bladder urothelial carcinoma (BLCA) and pheochromocytoma/paraganglioma (PCPG). Further, HOXC-AS1 predicted poor prognosis in patients with ACC, COAD, ESCA, acute myeloid leukemia (LAML), LGG, and THYM, but showed association with a survival advantage for KIRP. Furthermore, HOXC-AS2 predicted poor prognosis for ACC, BLCA, COAD, LGG, liver hepatocellular carcinoma (LIHC), THYM, but was associated with a survival advantage for KIRP. It was also observed that HOTAIR predicted poor prognosis for ACC, BRCA, cholangiocarcinoma (CHOL), COAD, GBM, KIRC, LGG, mesothelioma (MESO), and THYM, but was associated with a survival advantage for LAML. Hence, these results indicated that the different expression levels of the HOXATs show different prognostic associations in the different cancer types.
Moreover, this brief review of the literature indicated the existence of an association between HOXATs and OS in different cancer types. To further validate this aforementioned results, meta-analysis involving HOTTIP and HAGLR was also performed. The results obtained indicated that elevated HOTTIP  Figure S19). Moreover, given that several HOXATs influenced LGG prognosis, K-M survival analysis for OS and DSS (Disease Free Survival) were performed to further evaluate the prognostic value of the HOXATs in LGG, based on the TCGA database (Supplementary Figures S20,21). The results obtained showed that the overexpression of all 18 HOXATs could predict poor prognosis in patients with LGG, suggesting that HOXATs might serve as novel prognostic biomarkers for LGG.

Association Between HOXATs and the Immune Microenvironment in Pan-Cancer
We further examined the association between the expression levels of HOXATs and the infiltration levels of stromal and immune cells and tumor purity, which was defined based on stromal scores, immune scores, and estimate scores using the ESTIMATE algorithm ( Figures 3A,B). The higher stromal and immune scores indicated lower tumor purity, and as expected, the degree of association between the HOXATs and stromal scores for different cancer types varied over a wide range. In cancers, including BRCA, LGG, tenosynovial giant cell tumor (TGCT), thyroid carcinoma (THCA), and THYM, approximately half or more than half of the HOXATs showed positive association with the stromal score of the patients. However, some of the HOXATs showed negative association with stromal scores in patients with ACC, BLCA, and KIRP. Additionally, certain HOXATs displayed different correlations across various cancer types. For instance, HOXB-AS1 (p < 0.001) showed the highest correlation with stromal scores across the different cancer types followed by HOTAIRM1 (p < 0.001), and HOXA-AS2 (p < 0.05). More specifically, we observed a significant positive correlation between HOXB-AS1 and the stromal scores of CHOL, ESCA, LIHC, LGG, and THCA; HOTAIRM1 and the stromal scores of CHOL, LGG, LIHC, and THCA; HOXA-AS2 and the stromal scores of BRCA, TGCT, and THCA, as well as HOXD-AS2 and the stromal score of TGCT (p < 0.0001). Conversely, there was a significant negative correlation between HOXA11-AS and the stromal scores of ACC and pancreatic adenocarcinoma (PAAD), and HOTAIR and the stromal scores of ACC and sarcoma (SARC) (p < 0.001). Further, we also analyzed the correlation between HOXATs and the immune score. Thus, we observed a significant positive correlation between HOXB-AS1 and the immune scores of CHOL, diffuse large B-cell lymphoma (DLBC), LGG, LIHC, prostate adenocarcinoma (PRAD), and THCA, and between HOTAIRM1 and the immune scores of CHOL, LGG, LIHC, and THCA ( Figure 3B).
To further clarify the relationship between each of the HOXATs and the immune components, we determined the correlation between HOXATs and the immune infiltrates in the cancer tumors. Six classes of immune infiltrates corresponding either to tumor promotion or suppression, namely, C1 (wound healing), C2 (INF-r dominant), C3 (inflammatory), C4 (lymphocyte depleted), C5 (immunologically quiet), and C6 (TGFβ dominant), were identified. Patients characterized by C3 and C5 immune subtypes showed significantly better survival rates than those characterized by the other immune subtypes, whereas patients with C4 and C6 subtypes showed the least favorable survival, given that C1 and C2 as well as C4 and C6 showed associated with more aggressive immune subtypes . A comparison of the expression levels of HOXATs corresponding to the different immune infiltrate subtypes showed a strong relationship between the expression of HOXATs genes and the immunity of the patients (p < 0.001; Figure 3C). Most of the HOXATs showed a lower expression level in the C5 immune subtype, indicating that they possibly play tumor promoter roles. Conversely, HAGLR and HOXD-AS2 were highly expressed in all the six immune subtypes, showing higher expression in C5, which in future, can be used as a new direction for immunosuppression. Additionally, higher levels of HOXA11-AS, HOTTIP, HOXB-AS3, HOXB-AS4, PRAC2, and HOXC13-AS were found to be associated with types 1, 2, and 6 infiltrates (i.e., C1, C2, and C6), indicating the tumor promoter role of these gene members, given that patients belonging to these categories, who showed a higher cancer cell proliferation rate as well as TGFβ enrichment, had lower survival rates . Conversely, HOTAIRM1, HOXA-AS2, HOXA10-AS, HOXB-AS1, HOXC-AS1, HOXC-AS2, and HOXD-AS2 showed association with C3 and C6 immune subtypes, and the decreased expression in the C5 subtype indicated that the differential expression of HOXATs in the different immune subtypes may partially explain why HOXATs played contrasting roles in the prognosis of various cancers.

Correlation Between HOXAT Expression and Immune Marker Sets
The application of cancer immunotherapy has received significant attention in recent years. Here, we analyzed the relationship between the expression of HOXATs and 47 immune checkpoint genes in different immune infiltrated cells. Notably, the expression levels of most of the HOXATs showed a close relationship with the expression of immune markers in the various immune cells (Supplementary Figures  S22-S39). Specifically, the expression levels of the majority of resting T cells, NK cells, eosinophils, and mast cells exhibited a significant correlation with the expression of HOXATs. This observation suggested that in some tumors, HOXATs possibly play a role in modulating the pattern of tumor immunity by regulating the expression levels of these immune checkpoint genes. Interestingly, HOXATs, including HOXB-AS1, HOXC-AS2, HAGLR, HOTAIR, and HOTAIRM1 showed significant positive correlation with the expression levels of immune checkpoint genes, including CTLA-4 and PD-L1, suggesting that these HOXATs are involved in the regulation of immune responses in different cancer types ( Figures 4A,B). Therefore, these results confirmed the specific correlation between HOXAT expression levels and immune marker sets.

Correlation Between HOXAT Expression and Tumor Mutation Burden in Pan-Cancer
TMB, which is considered an essential factor that affects the occurrence and progression of tumors, has been shown to affect response to cancer immunotherapy (Mao et al., 2018). Thus, we examined its association with HOXAT expression levels ( Figure 4C). Interestingly, we observed that approximately half or more than half of the HOXATs were positively associated with TMB in cancers, including ACC, LGG, ESCA, HNSC, and MESO, but negatively correlated with TMB in SKCM. Notably, the higher expression levels of HOXATs in LGG showed their greater likelihood to bring about exposure of immune checkpoints, inducing a stronger anti-tumor immune response.

Association Between HOXATs and Tumor Stemness in Pan-Cancers
To conduct a preliminary study on the association between HOXATs and cancer stemness, tumor stemness based on RNA stemness scores (RNAss) and DNA stemness score (DNAss) were utilized ( Figures  5A,B); lower RNAss or higher DNAss are associated with higher immune infiltration. Thus, we observed that HOXATs were negatively associated with RNAss in BRCA, KIRC, KIRP, LGG, READ, SKCM, TGCT, and THCA, and positively associated with DNAss in CHOL, GBM, LGG, and STAD. Notably, with the exception of HOTAIRM1 and HAGLR, all the other HOXATs showed negative association with RNAss and DNAss in TGCT.

Potential Roles of HOXATs in Glioma Tumor
Immunity and the TME In this study, we also compared the expression levels of HOXAT genes across the immune subtypes in LGG and GBM based on the TCGA RNAseq database. First, the genes showed contrasting intra-tumor heterogeneities in the different immune subtypes of the different grades of glioma. Specifically, HOXB-AS3, HOXC13-AS, and HAGLR showed significant differences in expression level across C1 and C4 immune subtypes, indicating that HOXAT immune subtypes are minimally correlated with GBM (p < 0.05; Figure 6A). High HOXAT expression levels in GBM enhance immunosuppression and facilitate the escape of GBM from immune supervision. However, no significant association was observed in GBM in this regard. Similarly, our results were consistent with those of other studies, indicating that elevated expression levels of HOXATs, such as HOXA10-AS and HOXA-AS2, are associated with the malignancy of gliomas. Conversely, except for HOTTIP, HOXB-AS4, and HAGLR, the expression levels of the other 15 HOXATs were significantly different across the C1 and C4 immune subtypes in LGG, suggesting that higher HOXAT expression levels in LGG elicits a higher immune response, while a lower expression level induces immune escape (p < 0.001; Figure 6B). Next, an investigation of the correlation between immune checkpoint genes and the HOXATs in GBM and LGG showed the existence of a significant correlation between the majority of the immune checkpoint genes and the HOXATs in LGG compared with GBM (p < 0.001; Figures 7A,B), suggesting that HOXATs possibly regulate malignancy degree of gliomas by influencing the immune system. The correlation analysis between HOXATs and cancer stem cell-associated genes was investigated. It was found that HOXATs and CD133, SOX2 were differentially association FIGURE 4 | Association of the homeobox antisense transcripts (HOXATs) with immune cell marker expression and tumor mutation burden (TMB). Heatmaps represent the association of the HOXATs expression with PD-L1 and CTLA4 in the pan-cancer (A, B). For each pair, the top-right triangle represents the p-value, while the bottom-left triangle represents the correlation coefficient. *p < 0.05, **p < 0.01, ***p < 0.001. Heat map showing the correlation between the HOXAT expression and TMB (C). The numbers on the heatmap represent the Spearman correlation values, where red and blue represent positive and negative correlations, respectively. *p < 0.05, **p < 0.01, ***p < 0.001.
Frontiers in Molecular Biosciences | www.frontiersin.org November 2021 | Volume 8 | Article 767856 within glioma tissues of different grades ( Figure 7C). The advance of tumor grade, was related to the increasing of SOX2 and CD133 cells. To further verify the role of signaling activation in the aberrant expression of HOXATs in pan-cancer, functional annotation was performed using GSEA. Interestingly, based on Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, except for HOXA-AS3 and HAGLR, 16 out of the 18 HOXATs were found to be enriched in several immune response pathways, including primary immunodeficiency, antigen processing and presentation, and the intestinal immune network corresponding to IGA.

Critical Role of HOXB-AS1 in the Maintenance of the Properties of GSCs
To evaluate the association between HOXB-AS1 and the properties of GSCs in glioma, the expression level of HOXB-AS1 and the stemness markers that are associated with stemness in GBM and LGG were analyzed using the GEPIA online dataset ( Figures 8A-C). As expected, HOXB-AS1 expression was significantly higher in GBM than in lower grade glioma tissues. Additionally, the stemness marker, CD133, which was significantly upregulated in GBM compared with normal tissues, could well distinguish between different grades of gliomas, while the other upregulated stemness marker, SOX2, was not different for the different grades of glioma. Moreover, correlation analysis demonstrated a strong correlation at transcriptional level between HOXB-AS1 and the stemness markers, CD133 and SOX2, especially in GBM (p < 0.001; Figures 8D-G). Further in vitro cell function experiments, which primarily involved cultured GBM and LGG cells in the normal culture medium with adherent culture conditions, confirmed that the established GSCs have a significantly higher neurosphere formation capacity than non-GSCs (p < 0.05; Figures 8H,I). Furthermore, sphere formation assay revealed that HOXB-AS1 knockdown inhibited spheroid formation in GSCs, indicating that HOXB-AS1 supports cancer stem cell growth. In conclusion, HOXB-AS1 was identified as being critical for the maintenance of the properties of GSCs ( Figures 8J,K).

Inhibition of Proliferation, Migration, and Invasion in Glioma Cells Following HOXB-AS1 Knockdown
To further investigate the biological role of HOXB-AS1 in glioma cells, we depleted the expression of HOXB-AS1 using lncRNA Smart Silencer. As shown in Figures 9A,B, MTT assays were performed to detect cell viability. Thus, it was observed that U251 and HS683 cells transfected with HOXB-AS1 smart silencer exhibited a markedly increased rate of apoptosis compared with the control cells (p < 0.05). Further, transwell migration assay showed that HOXB-AS1 downregulation decreased U251 and HS683 migration (p < 0.05; Figures 9C,D), and similarly, wound healing assay indicated that HOXB-AS1 downregulation could significantly decreased invasion (p < 0.01; Figures 9E,F). These results indicated that HOXB-AS1 knockdown inhibited glioma cell proliferation, migration, and invasion.

DISCUSSION
Presently, it is well known that over 75% of the human genome is functional and encodes a large number of ncRNAs. A recent study provided strong evidence that lncRNAs play important roles in numerous biological processes, including chromatin modification, RNA processing, structural scaffold invasion modulation, and apoptosis (Zhou et al., 2017;Miyazato and Hayakawa, 2020). However, lncRNAs continue to be discovered and are yet to be functionally annotated. Therefore, studies with a focus on the elucidation of the biological functions of these lncRNAs and the identification of those that are disease-related are still required.
HOXATs play pivotal roles in tumor growth and metastatic spread (Meti et al., 2018). For instance, HOXB-AS1, which might function as a carcinogen and a promising therapeutic biomarker for the clinical treatment of patients with multiple myeloma, GBM, and endometrial carcinoma patients, has been shown to regulate cancer cell proliferation, migration, and invasion (Klebanov et al., 2019;Obaid et al., 2021). In this study, the expression levels of HOXATs in different cancer types were analyzed based on the TCGA database. The results of Cox regression analysis suggested that the high expression levels of half or more than half of the HOXATs was primarily associated with unfavorable outcomes in ACC, KIRC, LGG, and THYM. Additionally, certain HOXATs, including HOXA-AS2, HOXB-AS1, and HOTAIR possibly exerted distinct prognostic effects across the various cancer types (Quinn and Chang, 2016;Petermann et al., 2019).
Further, accumulating evidence also suggests that lncRNAs play important roles in immune regulation (Sharma et al., 2011;Sheng et al., 2018). For example, Li et al. identified immune lncRNAs as potential oncogenic biomarkers for pan-cancer characterization, based on the TCGA database (Sun et al., 2016). It has also been reported that the HOTAIR, influences the development of cancer invasion via immune responses (Wang et al., 2019a;Wang et al., 2019b). However, only a few examples have been identified thus far, this study represents the first comprehensive genome-wide analysis of HOXAT expression in pan-cancer, wherein the potential relationships between HOXATs and the TME in pan-cancers are explored via a bioinformatics analysis. Notably, HOXATs, such as HOTAIRM1, HOXA11-AS, HOXB-AS1, HOXB-AS4, HAGLR, and HOXC-AS1, showed significant correlation with immune cell infiltration and were likely to exhibit higher expression levels in immune cells, while the correlation between their expression levels and tumor dryness scores suggested that they possibly act primarily as tumor promoters in the TME. Additionally, as described above, we observed that HOXAT expression levels varied across the molecular subtypes within the same tumor type, and that the tumor suppression or promotion effect of the HOXATs may be subtype-specific. The mechanisms of immune evasion by cancer cells are the targets of immunotherapies, and immune checkpoint inhibitors, including antibodies that block the programmed cell death of protein-1 (PD-1)/PD-L1 and the cytotoxic T lymphocyte associated protein 4 (CTLA-4), are receptors that attenuate T cell response. These have brought about changes in the strategies for the management of several cancers (Wang et al., 2015). Meanwhile, the co-expression of HOXATs in immune checkpoints and immune cells across 33 different cancer types indicated that the expression levels of HOXATs, such as HOXB-AS1, HOXC-AS2, HAGLR, HOTAIR, and HOTAIRM1 are significantly positively correlated with the expression of CTLA-4, PD-L1, and other immunocyte markers as well as immune cells, such as resting T cells, NK cells, eosinophils, and mast cells. Additionally, the identified signaling pathways showed association with immunity. Specifically, the results of this study indicated that these HOXATs may be important factors in the evolution of cancer development and might be involved in the regulation of immune response.
During cancer progression, tumor cells can progressively lose their differentiated phenotypes and acquire stem and progenitor-like features, which are believed to be responsible for tumor recurrence after standard chemotherapy and radiotherapy (Zhou et al., 2017). Recently, it has been hypothesized in several studies that cancer stem-like-cells, also known as tumor-initiating cells, mediate tumor initiation, progression, and metastasis (Wu et al., 2017). Interestingly, most of the HOXATs showed significant correlation with cancer stem cell-like features (DNAss and RNAss), suggesting that they likely play a role in tumor-initiating cells and are related to tumor resistance to Silencing of HOXB-AS1 inhibits the proliferation and tumorigenicity of GSCs in U251and HS683. The effects of siRNA targeting of HOXB-AS1 on spheroid-forming ability in GBM and LGG cells were detected by sphere formation assay (K, J). Bar: 100 μm *p < 0.05, **p < 0.01, ***p < 0.001.
Frontiers in Molecular Biosciences | www.frontiersin.org November 2021 | Volume 8 | Article 767856 FIGURE 9 | the homeobox antisense transcripts (HOXATs) enhanced the progression of glioma cells. Knockdown of HOXB-AS1 inhibited glioblastoma proliferation, metastasis, and invasion. Cell viability was measured by MTT assay in U251 and HS683 cell lines, and cell viability of cells with HOXB-AS1 knockdown was significantly inhibited (A, B). Transwell migration assay: HOXB-AS1 knockdown significantly decreased the migration distance in U251and HS683 cells (C, D). Invasion capacity was measured after HOXB-AS1 downregulation by transwell assay in U251 and HS683 cells (E, F). Bar: 100 μm *p < 0.05, **p < 0.01, ***p < 0.001. drug treatment (Zhao et al., 2018a). HOXB-AS1 was identified as being critical for the maintenance of the properties of GSCs. The immune microenvironment of GBM consists of anti-tumor immune response components (Yao et al., 2019;Yuan et al., 2020). Even though little is known regarding immunosuppression in this regard, a better understanding of the involvement and regulation of the immune system is essential for the treatment of glioblastoma. In this study, to observe the effects of HOXATs on the tumor immune microenvironment in gliomas, GBM and LGG were also compared in multiple ways, including immune subtype expression, immune markers, stromal scores, immune scores, and estimate scores using the ESTIMATE algorithm. The results obtained revealed that the expression levels of HOXATs are associated with diverse immune infiltration levels as well as tumor purity levels in gliomas, with high levels of intra-tumor heterogeneity in the immune infiltration corresponding to the different grades of the gliomas. Moreover, the HOXATs showed overexpression in the C4 immune subtype of GBM, i.e., the overexpression of HOXATs largely matched the immune profile of the C4 subtype, suggesting that HOXATs may play a negative role in the microenvironment, leading to an immunosuppressive TME, with poor prognosis, when compared with LGG. Further, the HOXATs also showed lower correlation with GBM than LGG in terms of both the immune and stromal scores, indicating that they are closely associated with the tumor immune microenvironment of LGG, but do not elicit a high immune response. Furthermore, the tumor immune microenvironment of GBM was altered when the glioma developed into GBM. Even though there was a high HOXAT expression that elicited an immune response, the HOXATs did not show any association with the TME. Therefore, we hypothesized that the early knockdown of HOXATs could inhibit the development of gliomas. HOXB-AS1 knockdown inhibited glioma cell proliferation, migration, and invasion.
With the continued exploration of lncRNAs (Zhang et al., 2018a) and the development of RNA-centric approaches (Zhang et al., 2018b), we will be able to identify the interaction targets of lncRNAs in a native context . The identification of these targets will further improve our understanding of lncRNA-mediated pathway regulation and advance our understanding of lncRNA functions (Zhao et al., 2019). Simultaneously, the relationship between lncRNAs and TME represented an additional layer of tumor microenvironment complexity. Hence, it is important to clarify the roles of lncRNAs in TME. Additionally, the TME has a complex tissue-like structure with a rich phenotype as well as functional heterogeneity. Different concentrations of chemokines and cytokines, as well as interactions between the cells that make up the TME can activate epithelial-mesenchymal transition-related signaling pathways and control the generation of CSCs (Zhao et al., 2018b). In this study, we described the contribution of HOXATs in pan-cancer prognosis and TME modulation, demonstrating that HOXATs might serve as biomarkers for prognosis and immune response across cancers.
This study had some limitations. First, one of the major limitation of this study is that all of the experiments in this study were done in cultured cells. Second, the effects of the different HOXATs in the different cancer subtypes were not analyzed in detail, implying that further studies are needed in this regard to obtain more accurate results. These limitations notwithstanding, our results can significantly facilitate the clarification of the role of HOXATs in tumorigenesis, specifically with respect to immune response, the TME, and prognosis. This may be beneficial to researchers in this field with respect to the study of the roles of HOXATs and the development of personalized medicine for cancer treatment.

CONCLUSION
This first comprehensive analysis of HOXATs expression in 33 different cancer types, analysed the role of the tumour microenvironment, and compared results in different grades of glioma. Therefore, we hypothesized that the early knockdown of HOXATs could stop the development of gliomas. We believe that the continued investigation of the HOXATs identified here will aid in the development of better immunotherapies for human cancer and other diseases.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
QX, YGG, and YG conceived and designed the study. WS and QD collected the data. YGG and WS analyzed and interpreted the data. JX, ZW, and QD verified the data. YG and JX performed all the statistical analysis. QX, WS, YGG, and ZH prepared figures and prepared the manuscript. All authors have read and approved the final version.