- 1Guangzhou Panyu Central Hospital, Guangzhou, China
- 2Guangdong Provincial Hospital of Chinese Medicine, Guangzhou, China
Over the past few decades, researchers have become aware of the importance of non-coding RNA, which makes up the vast majority of the transcriptome. Long non-coding RNAs (lncRNAs) in turn constitute the largest fraction of non-coding transcripts. Increasing evidence has been found for the crucial roles of lncRNAs in both tissue homeostasis and development, and for their functional contributions to and regulation of the development and progression of various human diseases such as cancers. However, so far, only few findings with regards to functional lncRNAs in cancers have been translated into clinical applications. Based on multiple factors such as binding affinity of miRNAs to their lncRNA sponges, we analyzed the competitive endogenous RNA (ceRNA) network for the colorectal cancer RNA-seq datasets from The Cancer Genome Atlas (TCGA). After performing the ceRNA network construction and survival analysis, the lncRNA KCNQ1OT1 was found to be significantly upregulated in colorectal cancer tissues and associated with the survival of patients. A KCNQ1OT1-related lncRNA-miRNA-mRNA ceRNA network was constructed. A gene set variation analysis (GSVA) indicated that the expression of the KCNQ1OT1 ceRNA network in colorectal cancer tissues and normal tissues were significantly different, not only in the TCGA-COAD dataset but also in three other GEO datasets used as validation. By predicting comprehensive immune cell subsets from gene expression data, in samples grouped by differential expression levels of the KCNQ1OT1 ceRNA network in a cohort of patients, we found that CD4+, CD8+, and cytotoxic T cells and 14 other immune cell subsets were at different levels in the high- and low-KCNQ1OT1 ceRNA network score groups. These results indicated that the KCNQ1OT1 ceRNA network could be involved in the regulation of the tumor microenvironment, which would provide the rationale to further exploit KCNQ1OT1 as a possible functional contributor to and therapeutic target for colorectal cancer.
Introduction
Colorectal cancer (CRC) is a common malignant cancer and is the second-highest contributor to the worldwide incidence of cancer-related deaths (Miller et al., 2019). CRC develops sporadically from some inflammatory bowel diseases or hereditary cancer syndromes. The development of colorectal cancer is based on the adenoma-carcinoma sequence. So far, the molecular mechanism of the adenoma-carcinoma sequence has been only partly identified. CRC prognosis depends on factors related to the patient and treatment. The expertise of the treatment team is one of the most important determinants of the outcome. Early detection of CRC cells and cancer precursor cells significantly reduce morbidity and improve patient prognosis (Sung et al., 2021).
The mortality of colorectal cancer can be effectively reduced by screening for the cancer. The most common screening procedures include flexible sigmoidoscopy, double-contrast barium enema, fecal occult blood tests, and colonoscopy (Bibbins-Domingo et al., 2016). There is no consensus regarding which screening method is the best, and it appears that no one test is better than the other. Risk, cost, and effectiveness are the main factors to be considered when discussing different options (Issa and Noureddine, 2017). Undoubtedly, a complete colonoscopy has the advantages of allowing the entire colon to be assessed, material for a biopsy to be collected, and a polypectomy to be carried out all within the same examination time; however, it also has disadvantages of higher costs as well as risks, discomfort and inconvenience for the patient being examined. Therefore, it is for all practical purposes necessary to develop effective biomarkers of CRC for applications in screening, diagnosis and prognosis.
Most of the RNA transcribed from the human genome does not encode for proteins. Some of these non-coding RNAs (ncRNAs) have been found to dysregulate the normal expression of genes, including tumor suppressor genes and oncogenes. Therefore, ncRNAs are considered to be new promising targets for studying tumorigenesis. ncRNAs include long non-coding RNAs (lncRNAs), microRNAs (miRNAs), circular RNAs, and small interfering RNAs (siRNAs). An miRNA is a highly conserved non-coding RNA approximately 21–24 nucleotides in length, and interacts with target mRNAs to regulate gene expression. Some miRNAs have been reported to be involved in the occurrence and development of cancer (Hayes et al., 2014), and miRNA-based therapeutics have in fact reached the stage of clinical development (Toden et al., 2021). In recent years, the importance of lncRNA has gradually become recognized. An lncRNA is essentially an ncRNA with more than 200 base pairs. So far, some lncRNAs have been shown to play a major regulatory role in genetic regulation. Recent studies have shown multiple roles for lncRNAs in tumorigenesis (de Oliveira et al., 2019). The competitive endogenous RNA (ceRNA) hypothesis is related to lncRNAs and miRNAs, proposed by Salmena and others (Salmena et al., 2011), has been described as the “Rosetta Stele”, used to decode the RNA language in order to regulate RNA crosstalk and regulate biological functions. Many studies have shown that miRNA-mediated ceRNA regulation plays a crucial role in the occurrence and development of cancer (de Oliveira et al., 2019). Long non-coding RNA KCNQ1 opposite strand/antisense transcript one gene (KCNQ1OT1) were markedly upregulated in gastric cancer tissues and cells (Zhong et al., 2021). High expression of lncRNA KCNQ1OT1 was significantly related to poor survival in patients with CRC in a pan-cancer meta-analysis. The upregulation of KCNQ1OT1 in CRC tissues and cell lines was also confirm the important role of KCNQ1OT1 in CRC (Lin et al., 2021). In our current work, we investigated the KCNQ1OT1 by mining the CRC RNA-Seq dataset in The Cancer Genome Atlas (TCGA), tested its prognostic potential in a CRC cohort and constructed the KCNQ1OT1-related lncRNA-miRNA-mRNA ceRNA network. We aimed to provide the rationales to the further exploitation of KCNQ1OT1 as a possible functional contributor to and therapeutic target for CRC.
Materials and Methods
Differentially Expressed Genes From TCGA RNA-Seq Data
RNA-seq data of a colorectal cancer cohort (TCGA-COAD) was downloaded with gdc-client (version 1.3.0) from the data portal of TCGA (the dbGaP accession: phs000178. v11. p8. Release date: December 18, 2019). The RNA sequencing read counts of sample were obtained from TCGA. All these samples were sequencing by Illumina Genome Analyzer IIX. According to the bioinformatics pipeline for mRNA analysis in TCGA (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/), quality assessment was performed with FASTQC (version 0.11.8), the alignment was performed using a two-pass method with STAR (version 2.4.2a). The reads were aligned to the GRCh38 reference genome and then were quantified with HTSeq (version 0.6.1). CPM normalization was performed to correct library size differences between samples. As an initial filter, we retained only genes with log2CPM >1 in more than half of the samples. To compare expression levels between samples, for these genes alone, we then re-normalised raw count data by TMM (implemented in edgeR version 3.32.1) and transformed these by voom in limma (version 3.46.0) (Robinson and Oshlack, 2010; McCarthy et al., 2012; Ritchie et al., 2015). After carrying out this normalization, mRNAs, microRNA and lncRNAs differentially expressed in tumor group and solid normal tissues were identified. Thresholds for differential expression were set each as the adjusted p value less than 0.01 and |log2fold change|>1.
Construction of a lncRNA-miRNA-mRNA ceRNA Network
All of the differentially expressed genes were considered as candidates in the lncRNA-miRNA-mRNA ceRNA network construction, which was based on multiple factors such as binding affinity of miRNAs to their lncRNA sponges, RNA secondary structures and RNA-binding proteins, and the abundance and subcellular localization of ceRNA components (Qi et al., 2015). To construct the ceRNA network, the R package GDCRNATools (version 1.10.1) was used (Li et al., 2018). with five databases on miRNA-mRNA interactions including STarMir (version 2.2) (Kanoria et al., 2016), StarBase (version 2.0) (Li et al., 2014), miRcode (version 11) (Jeggari et al., 2012), spongeScan (version 1.0) (Furió-Tarí et al., 2016), and mirTarBase (version 7.0) (Chou et al., 2018) incorporated for the interactions analysis. The criteria for identifying competing lncRNA-mRNA interactions are: a) the strength of positive association between expression of lncRNA and its target mRNAs, b) the hypergeometric probability of shared miRNAs on the lncRNA-mRNA pair, c) the strength of regulation similarity of all shared miRNAs on the lncRNA-mRNA pair. Based on above criteria, Pearson’s correlation was used to measure the association between expression of lncRNA and mRNA. Hypergeometric distribution was measured by Fisher’s exact test. Regulation similarity was calculated based on the total number of the lncRNA-mRNA shared miRNAs, the Pearson’s correlation between the miRNA with lncRNA, as well as miRNA with mRNA (Li et al., 2018). After the construction, the ceRNA network was plotted using Cytoscape (version 3.7.0) software (Smoot et al., 2011).
Survival Analysis
We investigated the prognostic values of the main differentially expressed lncRNAs in the ceRNA network for CRC patients. Kaplan-Meier curves for survival analysis were depicted to present the survival-related lncRNAs. A log-rank test was performed to compare the survival distributions of samples grouped by lncRNA differential expression levels of the patient cohort.
Functional Annotation
To identify the biological function of the ceRNA network possibly contributing to tumor development, we performed a functional annotation analysis with multiple pathway databases, including Gene Ontology (GO) (Ashburner et al., 2000), Reactome (Haw et al., 2011), and Speed2 (Signalling Pathway Enrichment using Experimental Datasets 2) (Rydenfelt et al., 2020). A p value of less than 0.05 was considered to indicate statistically significant enrichment. The R package clusterProfiler V3.11 (Yu et al., 2012) was used for the GO and Reactome pathway enrichment analysis and visualization of significant modules. The R package SPEED2 was used for checking the upstream pathway activity from the genes in the ceRNA network, and a Bates test was used to calculate the test statistics for pathway enrichment.
Optimization of the ceRNA Network
There were more than 100 genes in the ceRNA network (Figure 2), so for practical applicability, we confirmed the target structural accessibility and selected the critical lncRNAs associated with the survival of CRC patients as the hub genes. Only the mRNAs both correlated with the hub gene and in the original ceRNA network were used to reconstruct the optimized network. Target structural accessibility for miRNA target recognition was calculated and visualized using STarMirDB and Sfold (version 2.2) (Kanoria et al., 2016; Rennie et al., 2019). The optimized network was visualized using Cytoscape.
Network Signature Analysis
Gene set variation analysis (GSVA) was used to determine the network expression level of each single sample, analogously to a competitive gene set test (Hänzelmann et al., 2013). By performing GSVA scoring, we were able to estimate the variation in the gene enrichments of the networks of the samples, and to do so independently for tumor and normal tissues. GSVA scores are designed to range from −1 to 1, with negative scores indicating relative decreases in network expression while positive scores indicated elevations. For practical applicability, only lncRNAs and mRNAs in networks were listed as the gene set signature for the GSVA assessment. The R package GSVA (v3.11) was used to calculate the GSVA score of networks over the samples in the CRC transcriptome dataset.
Estimation of Immune Cell Infiltration
Immune cells infiltrated in the tumor microenvironment play crucial roles in tumor invasion and metastasis. To estimate immune cell infiltration in samples with different network expression levels, the web-based tool ImmuCellAI was applied to calculate the abundance of 24 immune cell subsets via their gene expression profiles (Miao et al., 2020). The immune cell subsets estimated in this study included 18 T cells subsets: CD4+T cells, CD8+T cells, naïve CD4+T cells, naïve CD8+T cells, natural regulatory T (nTreg) cells, induced regulatory T (iTreg) cells, gamma delta T (γδ T) cells, central memory T (Tcm) cells, effector memory T (Tem) cells, natural killer T (NKT) cells, T helper 1 (Th1) cells, T helper 2 (Th2) cells, and T helper 17 (Th17) cells, cytotoxic T (Tc) cells, exhausted T (Tex) cells, type 1 regulatory T (Tr1) cells, follicular T helper (Tfh) cells, mucosal-associated invariant T (MAIT) cells, and other six immune cell subests: B cells, natural killer (NK) cells, monocytes, macrophages, neutrophils and dendritic cells (DCs).
The pipeline of this study is depicted in the flowchart shown in Figure 1.
Results
CRC Tumor Samples and Normal Samples Were Significantly Distinguished on the Basis of Differentially Expressed lncRNAs
We obtained the aligned read counts (GRCh38 (hg38) version) of tissue samples from 478 primary tumor, one metastatic tumor, one recurrent tumor and 41 normal solid tissue from TCGA. Samples from primary tumor, metastatic tumor and recurrent tumor were combined as tumor group. The age and sex of patients in normal group was matched to those in tumor group (Supplementary Table 1). There were total number of 60483 genes in this dataset. In prefiltering step, 14768 genes were kept for the downstream analyses. After TMM normalization and voom transformation, we explored the DGEs (differentially expressed genes) based on GLM (generalized linear model) likelihood ratio test, 2935 CRC tissue-specific mRNAs and 213 lncRNAs via the differential expression analysis of the TCGA-COAD dataset (Supplementary Figures S1A,B, Supplementary Table S2). The heatmaps of the differential expression of lncRNAs between CRC tissues and solid normal tissues are shown in Supplementary Figure S1C. The expression of 213 differentially expressed lncRNAs separated the tumor group and normal group clearly.
Construction of the lncRNA-miRNA-mRNA ceRNA Network
In the current work, certain lncRNAs and mRNAs were shown to co-express in the ceRNA networks. We built the ceRNA network based on the co-expression patterns of the lncRNAs, miRNAs and mRNAs. In total, 133 nodes and 288 edges constituted the ceRNA network. The preliminary network is shown in Figure 2.
Knowledge-Driven Pathway Analysis of ceRNA Network
According to the Reactome analysis, biological processes of the genes in the ceRNA network from tumor tissue were mainly related to MET pathways (Figure 3A). MET is a receptor tyrosine kinase (RTK), and like other related RTKs such as EGFR, MET can be activated by binding to its ligand, namely hepatocyte growth factor/scatter factor (HGF/SF), resulting in MET dimerization and trans-autophosphorylation.
FIGURE 3. Functional annotation for the genes in ceRNA network of CRC. (A). Reactome pathway enrichment of the genes in ceRNA network of CRC. (B). The top 20 enriched molecular in GO terms. (C). Speed2 pathway activity ranking. (D). CRC related lncRNAs similarity.
A total of 106 GO terms were extracted using the GO analysis (Supplementary Table S3). Of the 106 GO terms, 97 were Biological Process terms, eight were Cellular Component terms and one was a Molecular Function term. The BP GO terms have all been confirmed to be related to the regulation of neurogenesis and mesenchymal cell proliferation (Figure 3B) with, for instance, GO:0050768 (negative regulation of neurogenesis), GO:0010975 (regulation of neuron projection development) and GO:0051961 (negative regulation of nervous system development) involving neurogenesis, and GO:0010464 (regulation of mesenchymal cell proliferation) and GO:0051591 (mesenchymal cell proliferation) involving proliferation of human mesenchymal cells.
SPEED2 analyses allow one to infer upstream pathways. In the current work, the VEGF pathway was indicated from this analysis to be the signaling pathway that likely caused the genes in the ceRNA network to be deregulated (Figure 3C, Supplementary Table S4). We used another web-based tool LnCeVar-Cluster which provides cluster profiles between differential expressed lncRNAs and ceRNAs, especially the similarity profiles between differential expressed lncRNAs (Wang et al., 2020). The differential expressed lncRNAs clustering profile of TACG-COAD indicated that NEAT1, MALAT1, KCNQ1OT1 and LINC00969 had the similar expression pattern in CRC (Figure 3D). Functional annotations showed the ceRNA network genes to be closely related to the tumor pathogenesis in general.
KCNQ1OT1 Determined to Be the Critical Hub Gene for the Network Optimization
We chose the hub lncRNAs (degree >5) and their related miRNAs and mRNAs in the ceRNA network. We identified 22 lncRNAs involved in preliminary ceRNA networks (MIER3, MET, ANKRD13B, ATP11A, TTYH3, EPHB4, FOSL1, MIR17HG, TARBP1, COL1A1, SCML1, ATP2B1, H19, CPEB2, PMEPA1, KCNQ1OT1, SOX12, CDC7, HECW2, MALAT1, CSF1, PI15) (Figure 2). Then we analyzed the association of the hub lncRNAs with the clinically obtained survival data to identify the lncRNAs crucial to CRC prognosis. Only the lncRNA KCNQ1OT1 was significantly differentially expressed according to the log-rank test of survival analyses (Table 1). A Kaplan-Meier estimate showed poorer prognoses for patients with higher levels of KCNQ1OT1 (Figure 4).
FIGURE 4. Survival analysis of KCNQ1OT1 ceRNA network and gene co-expressions. Kaplan-Meier and ROC curves of CRC cohort in high and low KCNQ1OT1 expression level.
Of the lncRNAs present in the preliminary ceRNA network, only one gene, namely KCNQ1OT1, was indicated to be associated with survival in the CRC cohort. Thus, a correlation analysis was performed based on this critical KCNQ1OT1 hub gene. Based on a hypergeometric test and correlation analysis (Table 2), the mRNAs and lncRNAs not significantly correlated with KCNQ1OT1 or without a direct link with KCNQ1OT1 were removed from the network. The KCNQ1OT1-miRNA-mRNA subnetwork consisted of one centroid lncRNA node, nine miRNA nodes and nine mRNA nodes. This KCNQ1OT1 ceRNA network was reconstructed and visualized using Cytoscape (Figure 5). The expression level of the lncRNA KCNQ1OT1 was positively correlated with the expression levels of the nine mRNAs (Figure 6). This result was consistent with the ceRNA theory that the lncRNA regulated other RNA transcripts by competing for shared microRNAs. We checked the expression profile of KCNQ1OT1 and other hub network genes in tumor group and normal group in TCGA data set (Supplementary Table S2, Supplementary Figure S2). An additional searching in The Genotype-Tissue Expression (GTEx) portal (dbGaP accession number phs000424. vN.pN) which is a comprehensive public dataset for tissue-specific gene expression and regulation in non-diseased cohort for verification of the expression of KCNQ1OT1 in colon was conducted (Supplementary Figure S3). The results showed that the gene expression of KCNQ1OT1 in both transverse colon and sigmoid colon was higher than in blood. And the gene expression of KCNQ1OT1 across all tissues showed that, although KCNQ1OT1 did not represent the highest expression in colon, its expression level was not low either, consistent with what we found in TCGA-COAD: the KCNQ1OT1 expression in normal tissue was low (CPM: median: 2.41, mean: 3.46). But its expression level in tumor tissue in TCGA-COAD was higher (CPM: median: 4.73, mean: 10.67). All this data showed that KCNQ1OT1 had enough expression to physiologically function as a miRNA sponge.
Target structural accessibility for miRNA target recognition to lncRNA KCNQ1OT1 or mRNAs was calculated and visualized using Sfold and STarMirDB (Supplementary Figures S4, S5).
The KCNQ1OT1 ceRNA Network Signature Was Highly Expressed in CRC Transcriptomic Profiles
A GSVA assessment of network signatures in tumor and normal samples and using the TCGA-COAD dataset showed higher network signature GSVA scores for the tumor tissues than for the normal samples. This result was expected because we constructed this network from DEGs in the TCGA-COAD dataset. For validation, we applied the GSVA assessment of network signatures in three GEO CRC transcriptome datasets (GSE21510, GSE113513 and GSE32323). In these GEO datasets, KCNQ1OT1 ceRNA network signature GSVA scores were significantly higher for tumor samples than for normal samples, which indicated the disease specificity of this KCNQ1OT1 ceRNA network (Figure 7).
FIGURE 7. GSVA scores of KCNQ1OT1 ceRNA network for cancer and normal samples in four different CRC datasets.
Different Immune Cell Infiltration Levels for High- and Low-GSVA-Score Tumor Samples
The immune cell infiltration levels for 24 immune cell subsets in tumor samples from the TCGA-COAD dataset are shown in Figure 8. Overall, the profiles of immune infiltration varied significantly between the tumor samples with high GSVA scores and those with low scores (Figure 8A). Also, lower CD4+/CD8+ T-cell ratios were found in the tumor samples with high network GSVA scores than in the tumor samples with low scores. In the high-score group, along with the increase in the levels of CD8+ T cells in general was observed an increase in those of naïve CD8+ T cells. Notably, in contrast, lower levels of cytotoxic T cells were found for the high-score group than for the low-score group. Also, along with the relatively low levels of CD4+ T cells in the tumor samples with high scores, were relatively low levels of the T helper subsets Th2 and Th17. Other T cell subsets, namely Treg (nTreg and iTreg), Tex, Tem and MAIT cells were also downregulated. However, the B cells were upregulated. Regarding another lymphoid cell line, natural killer cell levels showed relatively low levels for the high-score group. Further, regarding the myeloid cell line, lower levels of monocytes, macrophages and neutrophils were also found for the high-score group.
FIGURE 8. The immune infiltration of 24 immune cells subsets in the tumor samples with high GSVA score and low scores. (A). The comparisons with statistically significant difference. (B). The comparisons without statistically significant difference.
Patient sex, tumor stage and age were also investigated. However, the network GSVA score was neither associated with sex (p = 0.8) nor with age (p = 0.77) nor with tumor stage (p = 0.78), as shown in Supplementary Figures S6A–C, respectively.
Discussion
There has been increasing evidence for a link between dysregulation of lncRNAs and cancers (Gutschner and Diederichs, 2012). For instance, studies have shown the lncRNA MALAT1 to be associated with the development and metastasis of cancer cells (Gutschner et al., 2013; Tripathi et al., 2013). HOTAIR to be implicated in cancer metastasis regulation by targeting the chromatin repressor polycomb protein (Gupta et al., 2010), and linc00673 to activate WNT/β-catenin signaling and aggravate lung adenocarcinoma by binding between casein kinase 1ε (CK1ε) and DEAD box RNA helicase DDX3 (Guan et al., 2019).
Recent studies have also demonstrated important roles played by lncRNAs in the tumorigenesis of CRC. Colorectal cancer associated transcript 1 (CCAT1) was reported to be a specific biomarker for CRC (Xiang et al., 2014), and to be expressed at high levels not only in pre-malignant conditions but also throughout the various disease stages of CRC (Ozawa et al., 2017). Further, CCAT2, a lncRNA derived from the human chromosome MYC-335 region, has been shown to enhance metastasis and invasion through the MYC pathway-regulating miRNAs miR-20a29 and miR-17–5p. These two CRC-specific lncRNAs were shown to be transcribed from the 8q24 region where previous studies have found common genetic variants related to the risk of CRC (Yang et al., 2019).
In the current work, we showed that the lncRNA KCNQ1OT1, an antisense lncRNA transcribed from the human chromosome 11p15.5 KCNQ1 locus, is related to the pathogenesis and progression of CRC. It has been shown in previous work to function as a cis-silencer of the imprinted KCNQ1 cluster and to be involved in the metastasis and proliferation of various tumors, such as hepatocellular carcinoma, cholangiocarcinoma, ovarian and breast cancer tumors (Feng et al., 2018; Luo and Jin, 2019).
In our research, there were 213 lncRNAs differentially expressed between tumor and normal tissue. Among these 213 lncRNAs, we identified 22 lncRNAs involved in ceRNA networks. We analyzed the survival association with these 22 lncRNAs. Only KCNQ1OT1 significantly associated with the clinically obtained survival data. We compared 228 patients expressing KCNQ1OT1 at high levels with 228 patients expressing it at low levels. According to our Kaplan-Meier survival analysis, the patients with overexpressed KCNQ1OT1 did not survive on average for as long as did those displaying low KCNQ1OT1 expression. This result was consistent with the latest research (Lin et al., 2021). KCNQ1OT1 was reported to be upregulated in CRC tissue according to a study by Li and others (Li F. et al., 2019). KCNQ1OT1 knockdown in HCT116 and SW480 CRC cells downregulated the expression of Atg4B, which has been shown to cleave LC3 and promote the formation of autophagosomes (Hemelaar et al., 2003). The viability of these cells decreased after being treated with oxaliplatin, which implicated KCNQ1OT1 in inducing autophagy protection and chemo-resistance. Moreover, the relationship between the upregulation of KCNQ1OT1 and poor prognosis of CRC patients also suggested that higher KCNQ1OT1 levels in patients make them resistant to chemotherapy or other anti-cancer treatments (Li Y. et al., 2019). These results suggested that KCNQ1OT1 might become a promising therapeutic target for use in CRC patients.
The development of an understanding of RNA-RNA interactions is expected to provide deep insights into gene regulatory networks potentially implicated in various cancer diseases. Our optimized ceRNA network showed that the critical hub gene KCNQ1OT1 indirectly regulated 9 other mRNAs, namely those encoding ATP-binding cassette subfamily B member 6 (ABCB6), Rho guanine nucleotide exchange factor 19 (ARHGEF19), helicase lymphoid-specific (HELLS), gamma-glutamylcyclotransferase (GGCT), cell division cycle 7 (CDC7), sex comb on midleg-like 1 (SCML1), collagen type VII alpha 1 chain (COL7A1), TNF-receptor-associated factor 5 (TRAF5) and ankyrin repeat domain 13B (ANKRD13B) (Figure 5). According to our correlation analysis, these 9 mRNAs showed a co-expression relationship with KCNQ1OT1. A GeneCards (www.genecards.org) (Stelzer et al., 2016) analysis showed these gene products to be involved in some cancer-related pathways, for instance TRAF5 in the RANKL/RANK (receptor activator of NFKB (ligand)) signaling pathway, COL7A1 in ERK signaling, HELLS in the AMPK enzyme complex pathway and BRCA1 pathway, ARHGEF19 in RET signaling, and GGCT in a cancer metabolism pathway. These mRNAs and KCNQ1OT1 influence each other’s level by competing for the same pool of microRNAs: miR-29c-3p, miR-29b-3p, miR-326, miR-330-5p, miR-29a-3p, miR-7-5p, miR-140-5p, miR-24-3p and miR-335-5p. Among these microRNAs, the binding sited of miR-7-5p, miR-29a-3p, miR-29c-3p, miR-140-5p, miR-326 and hsa-miR-335-5p with KCNQ1OT1 have been confirm by experimental studies (Hu et al., 2018; Sun et al., 2018; Cheng et al., 2020; Mu et al., 2020; Yao et al., 2020; Zhou et al., 2021).
As the GSVA method can be used to score a gene set signature and depends neither on the composition nor size of a dataset, we applied this method to measure the optimal network signature across different datasets. Samples with a high network score across the independent datasets were found, based on the results, to be particularly enriched in the CRC tumor group.
The identified subnetwork was reported previously to be significantly more reproducible between different disease cohorts then were individual marker genes (Chuang et al., 2007), and in our study, the results of all these datasets indicated dramatically higher network scores for the tumor samples than for the normal samples (Figure 7). These results suggested that the KCNQ1OT1 ceRNA network might play a role in the mechanism of CRC pathogenesis. This network signature could be a potentially way to distinguish CRC tissue form normal tissue.
To gain more knowledge about the function of the KCNQ1OT1 ceRNA network related to immune infiltration, we used ImmuCellAI to estimate the abundance of immune cells from individual samples. The estimation implied that the overall profiles of immune infiltration differed significantly between the tumor samples with high GSVA scores and those with low scores (Figure 8A, Supplementary Table S7, S8). Lower CD4+/CD8+ T-cell ratio s were indicated for the higher-score groups. Both the reduction of the CD4+ T cell population and increase of the CD8+ T cell population contributed to the lower ratio of CD4+/CD8+ T-cell. The CD4+/CD8+ T-cell ratio was considered as an immunostimulatory marker in the general population (Gojak et al., 2019). A low CD4+/CD8+ T-cell ratio has been shown to be an immune risk phenotype related to chronic inflammation, persistent immune dysfunction and immune senescence (Wikby et al., 2005). In some investigations, the CD4+/CD8+ T-cell ratio was described as a significant marker for prognostic prediction in HIV/AIDS patients (Castilho et al., 2019; Gojak et al., 2019; F; Li F. et al., 2019). Notably, in our study, the trend found for CD8+ T cells was opposite that found for cytotoxic T cells (cytokine-produced CD8+ T cells) in the network score comparisons. The generation of a functional cytotoxic T-cell response in general depends on the activation of Th1 cells. However, in the current work, the proportions of Th1 cells in the high- and low-GSVA-score groups showed no significant difference. Another possibility was that a recruitment of CD8+ T cells in the tumor tissue accompanied the reduction of cytotoxic T cells, and that this recruitment reflected a suppression of the cytotoxic machinery of the infiltrates, suggesting that the dysfunctional status of the effector cells was due to the microenvironment in the samples with high network scores.
Beside the immune infiltration, the network score differences in sex, tumor stage and age were also investigated. Neither any sex bias nor age association with the network scoring was found according to the statistical analysis (Supplementary Figure S6).
Our findings investigated a network signature for the lncRNA KCNQ1OT1, which was shown in the current work to be a representative of the ceRNA network transcribed in CRC tissues. This ceRNA network could be a potential regulator in colorectal cancer immunity. These findings indicated an oncogenic role for KCNQ1OT1 and its downstream target mRNAs in the pathogenesis and progression of CRC.
Conclusion
The results of our study indicated that the KCNQ1OT1 ceRNA network could be involved in the regulation of the CRC tumor microenvironment, providing a rationale for the further exploitation of KCNQ1OT1 as a possible functional contributor to and therapeutic target for CRC.
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 author.
Ethics Statement
Ethical review and approval was not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.
Author Contributions
JL wrote the manuscript. JD contributed the design of the study. WL and SL organized the references and databases. JD and JL explored the RNAseq data. JL and WL revised manuscript. All authors contributed to write and approve the final manuscript.
Funding
This study was supported by the Medical and Health Science and Technology Project of Panyu District, Guangzhou (No.2019-Z04-88; 2018-Z04-59), the Guangdong Provincial Bureau of traditional Chinese Medicine Research Project (No. 20201317), Special Funding for TCM Science and Technology Research of Guangdong Provincial Hospital of Chinese Medicine (no. YN2019QL12) and the Guangzhou Health and Family Planning Commission Program (No. 20181A011118).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
Appreciation to National Cancer Institute (NCI), the National Human Genome Research Institute (NHGRI) due to the joint efforts on The Cancer Genome Atlas (TCGA). We would like to express our gratitude to the contributors of GEO transcriptome datasets (GSE21510, GSE113513 and GSE32323). The Genotype-Tissue Expression (GTEx) Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.684002/full#supplementary-material
References
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25, 25–29. doi:10.1038/75556
Bibbins-Domingo, K., Grossman, D. C., Curry, S. J., Davidson, K. W., Epling, J. W., García, F. A. R., et al. (2016). Screening for Colorectal Cancer. JAMA 315, 2564–2575. doi:10.1001/jama.2016.5989
Castilho, J. L., Turner, M., Shepherd, B. E., Koethe, J. R., Furukawa, S. S., Bofill, C. E., et al. (2019). CD4/CD8 Ratio and CD4 Nadir Predict Mortality Following Noncommunicable Disease Diagnosis in Adults Living with HIV. AIDS Res. Hum. Retroviruses. 35, 960–967. doi:10.1089/AID.2019.0064
Cheng, P., Lu, P., Guan, J., Zhou, Y., Zou, L., Yi, X., et al. (2020). LncRNA KCNQ1OT1 Controls Cell Proliferation, Differentiation and Apoptosis by Sponging miR-326 to Regulate C-Myc Expression in Acute Myeloid Leukemia. neo 67, 238–248. doi:10.4149/neo_2018_181215N972
Chou, C.-H., Shrestha, S., Yang, C.-D., Chang, N.-W., Lin, Y.-L., Liao, K.-W., et al. (2018). miRTarBase Update 2018: a Resource for Experimentally Validated microRNA-Target Interactions. Nucleic Acids Res. 46, D296–D302. doi:10.1093/nar/gkx1067
Chuang, H. Y., Lee, E., Liu, Y. T., Lee, D., and Ideker, T. (2007). Network‐based Classification of Breast Cancer Metastasis. Mol. Syst. Biol. 3, 140. doi:10.1038/msb4100180
Feng, W., Wang, C., Liang, C., Yang, H., Chen, D., Yu, X., et al. (2018). The Dysregulated Expression of KCNQ1OT1 and its Interaction with Downstream Factors miR-145/CCNE2 in Breast Cancer Cells. Cell Physiol Biochem. 49, 432–446. doi:10.1159/000492978
Furió-Tarí, P., Tarazona, S., Gabaldón, T., Enright, A. J., and Conesa, A. (2016). spongeScan: A Web for Detecting microRNA Binding Elements in lncRNA Sequences. Nucleic Acids Res. 44, W176–W180. doi:10.1093/nar/gkw443
Gojak, R., Hadžiosmanović, V., Baljić, R., Zečević, L., Ćorić, J., and Mijailović, Ž. (2019). CD4/CD8 Ratio as a Predictor for the Occurrence of Metabolic Syndrome in HIV/AIDS Patients during 6 Months of cART Therapy. J. Med. Biochem. 38, 489–495. doi:10.2478/jomb-2018-0049
Guan, H., Zhu, T., Wu, S., Liu, S., Liu, B., Wu, J., et al. (2019). Long Noncoding RNA LINC00673-V4 Promotes Aggressiveness of Lung Adenocarcinoma via Activating WNT/β-catenin Signaling. Proc. Natl. Acad. Sci. USA. 116, 14019–14028. doi:10.1073/pnas.1900997116
Gupta, R. A., Shah, N., Wang, K. C., Kim, J., Horlings, H. M., Wong, D. J., et al. (2010). Long Non-coding RNA HOTAIR Reprograms Chromatin State to Promote Cancer Metastasis. Nature 464, 1071–1076. doi:10.1038/nature08975
Gutschner, T., and Diederichs, S. (2012). The Hallmarks of Cancer. RNA Biol. 9, 703–719. doi:10.4161/rna.20481
Gutschner, T., Hämmerle, M., Eissmann, M., Hsu, J., Kim, Y., Hung, G., et al. (2013). The Noncoding RNA MALAT1 Is a Critical Regulator of the Metastasis Phenotype of Lung Cancer Cells. Cancer Res. 73, 1180–1189. doi:10.1158/0008-5472.CAN-12-2850
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics. 14, 7. doi:10.1186/1471-2105-14-7
Haw, R. A., Croft, D., Yung, C. K., Ndegwa, N., D'Eustachio, P., Hermjakob, H., et al. (2011). The Reactome BioMart. Database(Oxford) 2011, bar031. doi:10.1093/database/bar031
Hayes, J., Peruzzi, P. P., and Lawler, S. (2014). MicroRNAs in Cancer: Biomarkers, Functions and Therapy. Trends Mol. Med. 20, 460–469. doi:10.1016/j.molmed.2014.06.005
Hemelaar, J., Lelyveld, V. S., Kessler, B. M., and Ploegh, H. L. (2003). A Single Protease, Apg4B, Is Specific for the Autophagy-Related Ubiquitin-like Proteins GATE-16, MAP1-LC3, GABARAP, and Apg8L. J. Biol. Chem. 278, 51841–51850. doi:10.1074/jbc.M308762200
Hu, H., Yang, L., Li, L., and Zeng, C. (2018). Long Non-coding RNA KCNQ1OT1 Modulates Oxaliplatin Resistance in Hepatocellular Carcinoma through miR-7-5p/ABCC1 axis. Biochem. Biophysical Res. Commun. 503, 2400–2406. doi:10.1016/j.bbrc.2018.06.168
Issa, I. A., and Noureddine, M. (2017). Colorectal Cancer Screening: An Updated Review of the Available Options. Wjg 23, 5086–5096. doi:10.3748/wjg.v23.i28.5086
Jeggari, A., Marks, D. S., and Larsson, E. (2012). miRcode: a Map of Putative microRNA Target Sites in the Long Non-coding Transcriptome. Bioinformatics 28, 2062–2063. doi:10.1093/bioinformatics/bts344
Kanoria, S., Rennie, W., Liu, C., Carmack, C. S., Lu, J., and Ding, Y. (2016). “STarMir Tools for Prediction of microRNA Binding Sites,” in RNA Structure Determination: Methods and Protocols, Methods in Molecular Biology. Editors D.H. Turner, and D.H. Mathews (New York, NY: Springer), 73–82. doi:10.1007/978-1-4939-6433-8_6
Li, F., Sun, Y., Huang, J., Xu, W., Liu, J., and Yuan, Z. (2019a). CD4/CD8 + T Cells, DC Subsets, Foxp3, and Ido Expression Are Predictive Indictors of Gastric Cancer Prognosis. Cancer Med. 8, 7330–7344. doi:10.1002/cam4.2596
Li, J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein-RNA Interaction Networks from Large-Scale CLIP-Seq Data. Nucl. Acids Res. 42, D92–D97. doi:10.1093/nar/gkt1248
Li, R., Qu, H., Wang, S., Wei, J., Zhang, L., Ma, R., et al. (2018). GDCRNATools: an R/Bioconductor Package for Integrative Analysis of lncRNA, miRNA and mRNA Data in GDC. Bioinformatics 34, 2515–2517. doi:10.1093/bioinformatics/bty124
Li, Y., Li, C., Li, D., Yang, L., Jin, J., and Zhang, B. (2019b). lncRNA KCNQ1OT1 Enhances the Chemoresistance of Oxaliplatin in colon Cancer by Targeting the miR-34a/ATG4B Pathway. Ott 12, 2649–2660. doi:10.2147/OTT.S188054
Lin, Z.-B., Long, P., Zhao, Z., Zhang, Y.-R., Chu, X.-D., Zhao, X.-X., et al. (2021). Long Noncoding RNA KCNQ1OT1 Is a Prognostic Biomarker and Mediates CD8+ T Cell Exhaustion by Regulating CD155 Expression in Colorectal Cancer. Int. J. Biol. Sci. 17, 1757–1768. doi:10.7150/ijbs.59001
Luo, Z. P., and Jin, H. (2019). Effects of LncRNA KCNQ1OT1 on Proliferation and Migration of Ovarian Cancer Cells by Wnt/β-Catenin. Eur. Rev. Med. Pharmacol. Sci. 23, 8788–8794. doi:10.26355/eurrev_201910_19273
McCarthy, D. J., Chen, Y., and Smyth, G. K. (2012). Differential Expression Analysis of Multifactor RNA-Seq Experiments with Respect to Biological Variation. Nucleic Acids Res. 40, 4288–4297. doi:10.1093/nar/gks042
Miao, Y. R., Zhang, Q., Lei, Q., Luo, M., Xie, G. Y., Wang, H., et al. (2020). ImmuCellAI: A Unique Method for Comprehensive T‐Cell Subsets Abundance Prediction and its Application in Cancer Immunotherapy. Adv. Sci. 7, 1902880. doi:10.1002/advs.201902880
Miller, K. D., Nogueira, L., Mariotto, A. B., Rowland, J. H., Yabroff, K. R., Alfano, C. M., et al. (2019). Cancer Treatment and Survivorship Statistics, 2019. CA A. Cancer J. Clin. 69, 363–385. doi:10.3322/caac.21565
Mu, L., Deng, X., Wu, K., and Yang, C. (2020). LncRNA KCNQ1OT1 Promotes Cisplatin Resistance of Osteosarcoma Cancer Cells through the miR-335-5p/β-Catenin axis. Genes Dis. doi:10.1016/j.gendis.2020.06.002
Oliveira, J. C., Oliveira, L. C., Mathias, C., Pedroso, G. A., Lemos, D. S., Salviano‐Silva, A., et al. (2019). LncRNAs in Cancer: Another Layer of Complexity. J. Gene Med. 21, e3065. doi:10.1002/jgm.3065
Ozawa, T., Matsuyama, T., Toiyama, Y., Takahashi, N., Ishikawa, T., Uetake, H., et al. (2017). CCAT1 and CCAT2 Long Noncoding RNAs, Located within the 8q.24.21 'gene Desert', Serve as Important Prognostic Biomarkers in Colorectal Cancer. Ann. Oncol. 28, 1882–1888. doi:10.1093/annonc/mdx248
Qi, X., Zhang, D.-H., Wu, N., Xiao, J.-H., Wang, X., and Ma, W. (2015). ceRNA in Cancer: Possible Functions and Clinical Implications. J. Med. Genet. 52, 710–718. doi:10.1136/jmedgenet-2015-103334
Rennie, W., Kanoria, S., Liu, C., Carmack, C. S., Lu, J., and Ding, Y. (2019). Sfold Tools for MicroRNA Target Prediction. Methods Mol. Biol. 1970, 31–42. doi:10.1007/978-1-4939-9207-2_3
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007
Robinson, M. D., and Oshlack, A. (2010). A Scaling Normalization Method for Differential Expression Analysis of RNA-Seq Data. Genome Biol. 11, R25. doi:10.1186/gb-2010-11-3-r25
Rydenfelt, M., Klinger, B., Klünemann, M., and Blüthgen, N. (2020). SPEED2: Inferring Upstream Pathway Activity from Differential Gene Expression. Nucleic Acids Res. 48, W307–W312. doi:10.1093/nar/gkaa236
Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA Hypothesis: the Rosetta Stone of a Hidden RNA Language? Cell 146, 353–358. doi:10.1016/j.cell.2011.07.014
Smoot, M. E., Ono, K., Ruscheinski, J., Wang, P.-L., and Ideker, T. (2011). Cytoscape 2.8: New Features for Data Integration and Network Visualization. Bioinformatics 27, 431–432. doi:10.1093/bioinformatics/btq675
Stelzer, G., Rosen, N., Plaschkes, I., Zimmerman, S., Twik, M., Fishilevich, S., et al. (2016). The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr. Protoc. Bioinformatics. 54, 1.30.1–1.30.33. doi:10.1002/cpbi.5
Sun, H., Li, Y., Kong, H., Dai, S., and Qian, H. (2018). Dysregulation of KCNQ1OT1 Promotes Cholangiocarcinoma Progression via miR-140-5p/SOX4 axis. Arch. Biochem. Biophys. 658, 7–15. doi:10.1016/j.abb.2018.09.019
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A. Cancer J. Clin. 71, 209–249. doi:10.3322/caac.21660
Toden, S., Zumwalt, T. J., and Goel, A. (2021). Non-coding RNAs and Potential Therapeutic Targeting in Cancer. Biochim. Biophys. Acta (Bba) - Rev. Cancer 1875, 188491. doi:10.1016/j.bbcan.2020.188491
Tripathi, V., Shen, Z., Chakraborty, A., Giri, S., Freier, S. M., Wu, X., et al. (2013). Long Noncoding RNA MALAT1 Controls Cell Cycle Progression by Regulating the Expression of Oncogenic Transcription Factor B-MYB. Plos Genet. 9, e1003368. doi:10.1371/journal.pgen.1003368
Wang, P., Li, X., Gao, Y., Guo, Q., Ning, S., Zhang, Y., et al. (2020). LnCeVar: A Comprehensive Database of Genomic Variations that Disturb ceRNA Network Regulation. Nucleic Acids Res. 48, D111–D117. doi:10.1093/nar/gkz887
Wikby, A., Ferguson, F., Forsey, R., Thompson, J., Strindhall, J., Löfgren, S., et al. (2005). An Immune Risk Phenotype, Cognitive Impairment, and Survival in Very Late Life: Impact of Allostatic Load in Swedish Octogenarian and Nonagenarian Humans. Journals Gerontol. Ser. A: Biol. Sci. Med. Sci. 60, 556–565. doi:10.1093/gerona/60.5.556
Xiang, J.-F., Yin, Q.-F., Chen, T., Zhang, Y., Zhang, X.-O., Wu, Z., et al. (2014). Human Colorectal Cancer-specific CCAT1-L lncRNA Regulates Long-Range Chromatin Interactions at the MYC Locus. Cell Res. 24, 513–531. doi:10.1038/cr.2014.35
Yang, T., Li, X., Montazeri, Z., Little, J., Farrington, S. M., Ioannidis, J. P. A., et al. (2019). Gene-environment Interactions and Colorectal Cancer Risk: An Umbrella Review of Systematic Reviews and Meta‐analyses of Observational Studies. Int. J. Cancer. 145, 2315–2329. doi:10.1002/ijc.32057
Yao, L., Yang, L., Song, H., Liu, T., and Yan, H. (2020). MicroRNA miR-29c-3p Modulates FOS Expression to Repress EMT and Cell Proliferation while Induces Apoptosis in TGF-Β2-Treated Lens Epithelial Cells Regulated by lncRNA KCNQ1OT1. Biomed. Pharmacother. 129, 110290. doi:10.1016/j.biopha.2020.110290
Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J. Integr. Biol. 16, 284–287. doi:10.1089/omi.2011.0118
Zhong, X., Wen, X., Chen, L., Gu, N., Yu, X., and Sui, K. (2021). Long Non‐coding RNA KCNQ1OT1 Promotes the Progression of Gastric Cancer via the miR‐145‐5p/ARF6 axis. J. Gene Med. 23, e3330. doi:10.1002/jgm.3330
Zhou, W., Li, H., Shang, S., and Liu, F. (2021). lncRNA KCNQ1OT1 Reverses the Effect of Sevoflurane on Hepatocellular Carcinoma Progression via Regulating the miR-29a-3p/CBX3 axis. Braz. J. Med. Biol. Res. 54, e10213. doi:10.1590/1414-431X2020e10213
Glossary
ABCB6 ATP binding cassette subfamily B member 6
ANKRD13B ankyrin repeat domain 13B
ARHGEF19 rho guanine nucleotide exchange factor 19
BP biological process
CC cellular component
CCAT1 colorectal cancer associated transcript 1
CDC7 cell division cycle 7
COL7A1 Collagen Type VII Alpha 1 Chain
ceRNA competitive endogenous RNA
CK1ε casein kinase 1ε
CRC colorectal cancer
DC dendritic cell
DEGs differentially expressed genes
γδ T gamma delta T cell
GGCT gamma-glutamylcyclotransferase
GO gene ontology
GSVA gene set variation analysis
GTEx The Genotype-Tissue Expression
HELLS helicase, lymphoid-specific
HGF hepatocyte growth factor
iTreg induced regulatory T cell
KM Kaplan-Meier
lncRNAs long non-coding RNAs
MAIT mucosal-associated invariant T cell
MF molecular function; miRNA, microRNA
ncRNAs non-coding RNAs
NK natural killer cell
NKT natural killer T cell
nTreg natural regulatory T cell
RANKL/RANK receptor activator of NFKB (ligand)
RTK receptor tyrosine kinase
SCML1 sex comb on midleg-like 1 (Drosophila)
SF scatter factor
siRNAs small interfering RNA
Speed2 Signalling Pathway Enrichment using Experimental Datasets 2
Tc cytotoxic T cell
TCGA The Cancer Genome Atlas
Tcm central memory T cell
Tex exhausted T cell
Tem effector memory T cell
Tfh follicular T helper cell
Th1 T helper 1 cell
Th2 T helper 2 cell
Th17 T helper 17 cell
TMM trimmed mean of M values
TRAF5 TNF receptor associated factor 5
Tr1 type 1 regulatory T cell.
Keywords: long non-coding RNA, KCNQ1OT1, ceRNA network, colorectal cancer, immunity
Citation: Liu J, Lv W, Li S and Deng J (2021) Regulation of Long Non-coding RNA KCNQ1OT1 Network in Colorectal Cancer Immunity. Front. Genet. 12:684002. doi: 10.3389/fgene.2021.684002
Received: 22 March 2021; Accepted: 09 September 2021;
Published: 22 September 2021.
Edited by:
Shaveta Kanoria, Wadsworth Center, United StatesCopyright © 2021 Liu, Lv, Li and Deng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jingwen Deng, dengjingwen@gzucm.edu.cn