Impact Factor 4.137 | CiteScore 4.28
More on impact ›

Original Research ARTICLE

Front. Oncol., 17 December 2019 | https://doi.org/10.3389/fonc.2019.01325

Novel Breast-Specific Long Non-coding RNA LINC00993 Acts as a Tumor Suppressor in Triple-Negative Breast Cancer

Shipeng Guo, Lei Jian, Kai Tao, Chen Chen, Haochen Yu and Shengchun Liu*
  • Department of Endocrine Breast Surgery, The First Affiliated Hospital of Chongqing Medical University, Chongqing, China

Background: Triple-negative breast cancer (TNBC) was characterized by breast cancers that do not express estrogen receptor (ER), progesterone receptor (PR), or human epidermal growth factor receptor (HER)-2 genes. TNBC patients are associated with a shorter median time to relapse and death for the lack of available treatment targets. Long non-coding RNAs (LncRNAs) have been reported to play an important role in the development of TNBC. We identified a novel breast-specific long non-coding RNA LINC00993, but less was known about its expression pattern and functional role in TNBC.

Methods: LINC00993 RNA expression was detected across different types of clinical breast cancer samples by using qRT-PCR. Bioinformatic methods “guilt by association” and gene set enrichment analysis (GSEA) were used to predict LINC00993 functions. Subcellular localization of LINC00993 in cells was detected by RNA fluorescence in situ hybridization (FISH). Effect of LINC00993 on cell growth was measured by plate colony formation assays, typical growth curve, and an in vivo tumor model. Cell cycle analysis was done by flow cytometry analysis. Key cell cycle regulators were detected by Western blot.

Results: LINC00993 was largely downregulated in TNBC, and higher expression indicated better outcome. LINC00993 located mainly in the nucleus. LINC00993 suppressed TNBC growth both in vitro and in vivo. LINC00993 was predicted to be involved in cell cycle pathways by using “guilt by association” and GSEA methods. Key cell cycle regulators like p16INK4A, p14ARF, p53, and p21 were affected by LINC00993 overexpression.

Conclusions: A new breast-specific lincRNA LINC00993 was identified with a tumor-suppressive feature and with prognostic value. This is the first research on LINC00993 function. Our results suggest that controlling LINC00993 level may be beneficial for breast cancer treatment.

Introduction

Triple-negative breast cancers (TNBCs) are breast cancers that do not express three receptor genes: estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor (HER)-2 genes (1). These cancers do not respond to endocrine therapy or other available targeted agents, thus posing a greater clinical treatment challenge (2). Patients diagnosed with TNBC showed a shorter median time to relapse and death (3). However, the mechanism underlying TNBC pathogenesis is not completely understood.

The majority of the human genome are transcribed as non-coding RNAs (ncRNAs) with <2% of the human genome encoding proteins (4). Among these ncRNAs, transcripts of more than 200 nucleotides that are not able to translate into proteins are defined as long non-coding RNAs (lncRNAs). These lncRNAs can be intergenic (long intergenic non-coding RNA, lincRNA), intronic, or natural antisense transcripts. Their transcription are controlled by promoters with divergent enhancers (4). The specific functions of lncRNAs are largely unknown, and they are less conserved across multiple species. LncRNAs can regulate gene expression in cis-acting or in trans-acting manner (5). LncRNAs are reported to regulate gene expression in diverse biological and pathological contexts, including cancer (6).

By using RNA sequencing (RNA-Seq) and microarray technology, a large number of lncRNAs have been identified to be aberrantly expressed in breast cancers, indicating that those lncRNAs might play a pivotal role in breast cancer (7). Some of these identified lncRNAs were related to breast cancer development (8, 9). Metastasis associated with lung adenocarcinoma transcript 1 (MALAT1) was highly expressed in TNBC tissues, which promoted the proliferation of TNBC cells in vitro and promoted TNBC growth and metastasis in vivo (10). Antisense oligonucleotides (ASOs) strategy was able to inhibit MALAT1 in a well-characterized mouse model of luminal B breast cancer (11). Those methods were adopted by several researchers to inhibit cancer (12).

Most lncRNAs found in breast cancer are reported to be elevated in tumor tissues compared to normal tissues (13).Functions of lncRNAs are largely unknown, especially for those whose expressions were suppressed in cancers. We previously reported that intergenic lncRNA LINC00993 was greatly downregulated in TNBC (14). A genome-wide transcriptional survey to examine lncRNAs in 995 breast tissue samples found that LINC00993 was the most credible downregulated lncRNA in ER-positive breast cancer compared with ER-negative breast cancer (7). However, the role of LINC00993 in breast cancer is totally unknown.

In this study, we explored the expression pattern of LINC00993 in clinical breast cancer samples and in The Cancer Genome Atlas (TCGA). LINC00993 was found to be breast-specific and was downregulated in breast cancer. We also studied the function of LINC00993 in breast cancer cells both in vitro and in vivo. Our results suggest that LINC00993 may serve as a tumor suppressor in breast cancer.

Results

Guilt by Association Analysis and GSEA Methods Indicated LINC00993's Role in Cell Cycle Regulation

LINC00993 is a lincRNA, which maps on chromosome 10p11.21. LINC00993 was not conserved across multiple species, with only chimpanzee sharing the same sequence with human (Figure 1A). The closest gene to LINC00993 is ANKRD30A, which encodes a breast-specific protein (15). LINC00993 is 666 nucleotides (nt) long consisting of five exons (Figure 1B). By analyzing data from the GTEx project (16), LINC00993 was expressed in a tissue-specific pattern, mainly expressed in breast, which is similar to ANKRD30A (Figure 1C). We have previously shown that LINC00993 was downregulated in TNBC tissues, but its function in TNBC is largely unknown.

FIGURE 1
www.frontiersin.org

Figure 1. Information about LINC00993 from public websites. (A) Location of LINC00993 in human genome. Conservation of LINC00993 in multiple species. (B) Construction and sequence of LINC00993. (C) Detailed LINC00993 expression in different tissues shown by GTEx data.

It has been a challenging task to define lncRNA's functions. A “guilt by association” method has been used to predict ncRNA functions (17, 18). To examine the function of LINC00993 in a tissue-specific way, breast cancer data containing 1,208 breast cancer samples from TCGA database (1,096 are from cancer tissues and 112 are from peritumor tissues) were downloaded. We performed correlation analysis between LINC00993 and all 19,547 coding genes in those 1,096 samples by Spearman method. Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed on the top 5% genes (978) according to correlation index. The GO analysis results contain three parts: cellular component, biological process, and molecular function (19). The results showed that LINC00993-related genes took part in the formation of spindle, chromosomal region, and other cellular components, which are mainly located in the nucleus (p < 0.001; Figure 2A). Those genes were involved in chromosome segregation, nuclear division, and other biological processes, which also happened in the nucleus (p < 0.001; Figure 2B). The molecular function of GO analysis showed that LINC00993-related genes were involved in microtubule motor activity (p < 0.001; Figure 2C). The KEGG pathway analysis on those LINC00993-related genes showed that those genes mainly enriched in cell cycle pathway (p < 0.001; Figure 2D). As shown in Figure 2B, some of the genes are involved in the regulation of cell cycle-phase transition (p < 0.001).

FIGURE 2
www.frontiersin.org

Figure 2. Annotation of LINC00993 via guilt-by-association correlation analysis. Gene Ontology (GO) analysis of LINC00993-related genes in cellular component (A), biological process (B), molecular function (C). (D) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.

Gene set enrichment analysis (GSEA) is a statistical method to determine if predefined sets of genes are differentially expressed in different phenotypes (20). To analyze which pathway LINC00993 might be involved with, 186 KEGG gene sets from Molecular Signatures Database (MSigDB) were used in GSEA. Most enriched pathways were inhibited by LINC00993 (Figure S1A). Cell cycle was the most significant gene set enriched with a normalized enrichment score as −2.6668 (Supplementary Supplementary Figure S1B). The enrichment score plot showed that LINC00993 negatively regulates cell cycle pathway (Supplementary Figure S1C). For each gene set, a 1,000 times permutation according to gene labels was done to produce a gene set null distribution. As shown in Supplementary Figure S1D, the normal Enrichment Score (ES) for the gene set ranges from −0.4 to 0.4, the ES for gene set cell cycle is −0.674 (p < 0.001, FWER < 0.001, FDR < 0.001) (Supplementary Figure S1D). A heat map was drawn based on the expression of genes in the mentioned gene set (Supplementary Figure S1E).

Hallmark gene sets summarize and represent specific well-defined biological states or processes and display coherent expression (21). To examine what potential target genes LINC00993 regulates, 50 hallmark gene sets were analyzed. There were 17 of 50 hallmark gene sets enriched (Supplementary Figure S2A). LINC00993 positively regulates estrogen response (p < 0.001, FWER < 0.01, FDR < 0.001), which is consistent with our data that LINC00993 and estrogen receptor (ESR)1 had a strong correlation. The most enriched hallmark data set was E2F targets with a NES at −3.12 (p < 0.001, FWER < 0.01, FDR < 0.001) (Supplementary Figure S2B). The running score plot showed that LINC00993 negatively regulated E2F targets (Supplementary Figure S2C). The null distribution plot showed that the result was reliable (Supplementary Figure S2D). A heat map was drawn based on the expression of genes in the E2F targets gene set (Supplementary Figure S2E).

Taken together, those data suggested that LINC00993 might mainly function in the nucleus and was involved in the regulation of the cell cycle.

LINC00993 Suppressed the Growth of TNBC Cells in vitro

Next, we intended to explore the function of LINC00993 in breast cancer. To choose appropriate cell lines to manipulate the expression of LINC00993, a couple of breast cancer cell lines were detected. LINC00993 was relatively higher in MCF-7 cells when compared to those in MDA-MB-231 cells and BT-549 cells (Supplementary Figure S3A). Lentivirus shRNA method was adopted to knock down LINC00993 in MCF-7, but no significant fold change of LINC00993 was detected in shRNA-stable MCF-7 (Supplementary Figures S3B–F). We thus decided to overexpress LINC00993 in MDA-MB-231 and BT-549, which are two TNBC cell lines that showed a lower LINC00993 expression level. To confirm the location of LINC00993 in TNBC cells, total RNA, cytoplasmic RNA, and nuclear RNA were isolated from MDA-MB-231 cells. Expression of LINC00993 was evaluated in these fractions by qRT-PCR assays. Our results showed that LINC00993 expression was located mainly in the nucleus compared to that in the cytoplasm (16.82 ± 1.69 vs. 1.08 ± 0.03, p < 0.01; Figure 3A). Results from RNA fluorescence in situ hybridization (FISH) assay also confirmed that most of the LINC00993 located in the nucleus in MDA-MB-231 cells (Figure 3B, red arrow), and some LINC00993 located around nuclear membrane (Figure 3B, white arrow). Those results indicated that LINC00993 mainly resided in the nucleus.

FIGURE 3
www.frontiersin.org

Figure 3. Location of LINC00993 in breast cancer cells. (A). Location of LINC00993 detected by qRT-PCR. Nuclear and cytosolic compartments of MDA-MB-231 cells were extracted. Nu means nuclear, CP means cytoplasm. ***P < 0.001, based on Student's t-test. Data were presented as mean ± SEM. (B). MDA-MB-231 cells were labeled with a LINC00993 probe via RNA fluorescence in situ hybridization (FISH) and counterstained with DAPI: 4',6-diamidino-2-phenylindole (DNA to visualize the nucleus). A representative image is shown.

We then aimed to elucidate the function of LINC00993 in TNBC. Recombinant adenovirus vector expressing LINC00993 under control of a Murine cytomegalovirus (MCMV) promotor and RFP under a Cytomegalovirus (CMV) promotor were constructed (Figure 4A). Fluorescent microscopic analysis of the expression levels of Red fluorescent protein (RFP) confirmed that the infection efficiency was high (Figure 4B), and qRT-PCR assay showed high LINC00993 expression with a 6,000-fold change in RNA expression in adenovirus-infected TNBC cell line MDA-MB-231 (Figure 4C).

FIGURE 4
www.frontiersin.org

Figure 4. LINC00993 suppresses the growth of triple-negative breast cancer (TNBC) cells in vitro. (A) Structure of LINC00993 expression plasmid for adenovirus. (B) LINC00993 expression adenovirus infection efficiency showed by fluorescence microscope in MDA-MB-231 cells. Original magnification, ×400. Scale bars, 50 μm. (C) Expression of LINC00993 detected by qRT-PCR. MDA-MB-231 cells were infected by adenovirus for 24 h, and RNA was extracted. (D) Image of clone formation assay. (E) Number of clones were counted 2 weeks after plantation. (F) MDA-MB-231 cells were planted into 24-well plates. Twenty-four hours later, adenovirus was added to each well. Three wells of cells were digested and counted every 24 h. (G) LINC00993 expression caused apoptosis shown by TUNEL assay. Green points reflected apoptosis, and we used DAPI to stain DNA. Positive control cells were treated with DNase I, negative control cells were collected without adding TUNEL reaction buffer. Original magnification, ×100. Scale bars, 50 μm. (H) Effect of LINC00993 on cell cycle detected by flow cytometry. (I) Flow cytometry cell cycle results shown in a bar plot. (J) Invasive ability tested by Transwell assay. Twenty thousand cells were plated in each well, cells were observed after 24 h of incubation. (K) Bar plot for the number of cells that migrated across the membrane. ***P < 0.001, based on Student's t-test. Data were presented as mean ± SEM.

To study the role of LINC00993 in cell proliferation, plate colony formation assays were performed with LINC00993 overexpressed MDA-MB-231 cells. Overexpression of LINC00993 decreased the colony formation ability with the number of cell clones dropped from 245 ± 41.58 to 115 ± 5 (p < 0.01; Figures 4D,E). LINC00993 decreased the growth rate of MDA-MB-231 (p < 0.01; Figure 4F). TUNEL assay showed that LINC00993 overexpression induced MDA-MB-231 cell apoptosis (18 ± 1.0 vs. 3.3 ± 0.5, p < 0.01; Figure 4G). Next, flow cytometry analysis of PI-stained MDA-MB-231 cells was used to examine LINC00993's effect on cell cycle. Overexpression of LINC00993 significantly increased the proportion of cells in the G0/G1 phases (84.26 ± 1.86 vs. 56.89 ± 2.06, p < 0.001), while decreasing cells in S phases (15.96 ± 8.7 vs. 32.08 ± 1.60, p < 0.01; Figures 4H,I), suggesting a potential role of LINC00993 in cell cycle regulation, which was in line with the results from KEGG pathway analysis. As metastasis is an important feature of TNBC, Transwell migration assay was performed to detect the migrating ability of MDA-MB-231 cells. The results showed that LINC00993 strongly suppressed cell migration with less cell transport through membrane (100 ± 8.7 vs. 602 ± 20.56, p < 0.001; Figures 4J,K). Similar results were found in BT-549 cells (Supplementary Figures S4A–G). Those results indicated that LINC00993 played a tumor suppressor role in TNBC cells by inhibiting cell proliferation.

LINC00993 Suppressed the Growth of TNBC Cells in vivo

To confirm the in vitro results, nude mice xenograft models of MDA-MB-231 cells were established. About 2 weeks after plantation with tumor size at ~100 mm3, the mice were intratumorally injected with Ad-LINC00993, Ad-RFP, or PBS. The volume of the tumor was measured twice a week till 5 weeks. The average tumor volume in Ad-LINC00993 group was much smaller than that in the Ad-RFP group (287.14 ± 68.73 vs. 1,628.57 ± 475.09, p < 0.01; Figure 5A). Five weeks after cell injection, the mice were sacrificed and the weight of tumors were examined (Figure 5B). The average weight of Ad-LINC00993-treated tumors was markedly lower than control (0.34 ± 0.05 vs. 1.39 ± 0.27, p < 0.01; Figure 5C). The results showed that ectopic expression of LINC00993 in MBA-MD-231 xenografts inhibited tumor growth, thus LINC00993 played a tumor suppressor role in TNBC in vivo.

FIGURE 5
www.frontiersin.org

Figure 5. LINC00993 suppresses the growth of triple-negative breast cancer (TNBC) in vivo. (A) The average tumor volume in Ad-LINC00993 group was much smaller than that in Ad-RFP or PBS group. (B) Image of nude mouse burden with breast cancer. Five weeks after cell injection, the mice were sacrificed, and the weight of tumors were examined. (C) The average weight of Ad-LINC00993-treated tumors was markedly lower than control. *P < 0.05 and **P < 0.01 based on Student's t-test. Data were presented as mean ± SEM.

LINC00993 Functions Through Cell Cycle-Related Genes

Results from KEGG pathway analysis indicated that LINC00993 regulates the cell cycle pathway. Our previous data showed that LINC00993 induced G0/G1 arrest. We were interested to check whether key G0/G1 arrest regulators were changed after LINC00993 overexpression.

E2F was pivotal in the control of G1 to S phase transition. p16INK4A was reported to inhibit Cyclin D-CDK4/6 complex, which could promote translocation of E2F from cytoplasm to the nucleus. On the other hand, translocation of E2F was also controlled by Cyclin E-CDK2 complex (22, 23). p21 was reported to take part in this process (24). A couple of studies had proved that p14ARF-p53-p21 signal pathway was involved in the control of G1/S transition (25, 26). We then decided to detect the expression of those key factors after overexpression of LINC00993. Our results showed that p16INK4A was significantly increased 8-fold (0.18 ± 0.007 vs. 1.48 ± 0.10, p < 0.001; Figures 6A,B). p21 expression was raised after LINC00993 expression (0.4 ± 0.02 vs. 1.04 ± 0.02, p < 0.001; Figures 6A,D). LINC00993 significantly increased p14ARF to 1.6-fold change (0.62 ± 0.02 vs. 0.93 ± 0.01, p < 0.001; Figures 6A,C). At the same time, p53 was also elevated (0.59 ± 0.05 vs. 1.02 ± 0.01, p < 0.001; Figures 6A,E). Taken together, our results showed that LINC00993 might inhibit Rb/E2F pathway by upregulating p16INK4A, p14ARF.

FIGURE 6
www.frontiersin.org

Figure 6. LINC00993 regulates cell cycle-related genes. (A) Western blot bands show cell cycle-related proteins infected by LINC00993 overexpression. Integrated optical density (IOD) was measured by Photoshop. The experiment was performed at least three times independently. (B–E) Bar graph for Western blot results of p16INK4A, p21, p14ARF, p53 from (A). ***P < 0.001 based on Student's t-test. Data were presented as mean ± SEM. (F) Schematic model of LINC00993 in breast cancer. LINC00993 upregulates p16INK4A, p14AR, inducing G0/G1 arrest via E2F pathway.

RNA-protein binding assays was performed to investigate potential binding partner interacting with LINC00993. In vitro-transcribed biotinylated LINC00993 or its antisense control RNA was incubated with cell extracts from MDA-MB-231, and potential RNA-protein complexes were captured using streptavidin magnetic beads and subjected to mass spectrometry (Supplementary Figure S5A). The most abundant RNA binding protein identified in these pull-downs was heterogeneous nuclear ribonucleoprotein A1 (hnRNPA1) (Supplementary Figure S5B). The binding between hnRNPA1 and LINC00993 was confirmed by Western blotting analysis, showing that sense LINC00993 could pull down hnRNPA1 while its antisense RNA showed less ability to bind hnRNPA1 (Supplementary Figure S5C). Our results clearly demonstrated that LINC00993 interacted with hnRNPA1. hnRNPA1 belongs to the hnRNP family. Some of the hnRNP family members were reported to regulate lncRNAs (27, 28). hnRNPA1 could bind and result in instability of mature p16INK4A and p14ARF mRNAs (29). The lncRNA UCA1 was reported to bind and sequester hnRNPA1, preventing its effect on stability of p16INK4A and p14ARF mRNAs (29). LINC00993 thus might interact with hnRNPA1 this way, regulating the cell cycle-related gene expression. A possible schematic model of LINC00993 function in breast cancer was shown in Figure 6F.

Validation of LINC00993 Expression in Breast Cancer Samples

To study the expression pattern of LINC00993 in clinical breast cancer samples, LINC00993 was examined in 98 tumor and peritumor paired clinical tissues by qRT-PCR, details of those patients were shown in Table 1. Compared to each paired peritumor tissue, the expression of LINC00993 was suppressed in breast cancer tissues in about 76.5% of the 98 patients (Figure 7A). Subgroup analysis showed that LINC00993 was downregulated in breast cancer tissues when compared to its paired peritumor tissues in TNBC group (n = 25, p < 0.001; Figure 7B). Similar results were observed in HER2 group (n = 31, p < 0.001; Figure 7C). No significant changes of LINC00993 levels were found in Luminal A group (n = 19, p > 0.05) and Luminal B group (n = 23, p > 0.05; Figures 7D,E).

TABLE 1
www.frontiersin.org

Table 1. Clinicopathological features of the study cohort (n = 98).

FIGURE 7
www.frontiersin.org

Figure 7. LINC00993 expression in clinical samples. (A) LINC00993 expression in 98 tumor and peritumor paired tissues by qRT-PCR. (B) LINC00993 expression in triple-negative breast cancer (TNBC) group. (C) LINC00993 expression in HER2 group. (D) LINC00993 expression in Luminal A group. (E) LINC00993 expression in Luminal B group. ns: not significant, ***P < 0.001 based on Student's t-test.

To confirm these findings, Breast Invasive Carcinoma (BRCA) data set from TCGA project was analyzed. Breast cancer subtype data were generated by using PAM50 predictor bioclassifier (30). Paired tumor and peritumor samples were extracted, and LINC00993 was found markedly downregulated in basal-like group compared to that in Normal group (n = 20, p = 0.007; Figure 8A). No significant difference was found in Her2 group (n = 17, p = 0.818; Figure 8B) and in Luminal B group (n = 29, p < 0.0001; Figure 8D). However, in Luminal A group (n = 330), LINC00993 was slightly elevated (n = 40, p = 0.031; Figure 8C).

FIGURE 8
www.frontiersin.org

Figure 8. LINC00993 expression in The Cancer Genome Atlas (TCGA) database. (A) Paired boxplot of LINC00993 expression in triple-negative breast cancer (TNBC) group. (B) Paired boxplot of LINC00993 expression in HER2 group. (C) Paired boxplot of LINC00993 expression in Luminal A group. (D) Paired boxplot of LINC00993 expression in Luminal B group. (E) Correlation analysis between LINC00993 and estrogen receptor (ESR)1 expression in GTEx mammary tissues using the Spearman method. (F) Correlation analysis between LINC00993 and ESR1 expression in TCGA BRCA tumor tissues using the Spearman method. (G) Overall survival Kaplan-Meier analysis was performed using the R2 database. (H) Relapse-free survival Kaplan-Meier analysis was performed using the R2 database.

As LINC00993 was downregulated in TNBC group but upregulated slightly in Luminal A group, we tested whether LINC00993 expression correlates with ER status or not. ER was encoded by ESR1 gene, a correlation analysis was carried out by using GEPIA. LINC00993 and ESR1 showed a high positive correlation index at 0.82 in GTEx mammary tissues (p < 0.0001; Figure 8E). In TCGA breast cancer tissues, they still showed a high positive correlation index at 0.45 (p < 0.0001; Figure 8F).

Due to the low expression level of LINC00993 in breast cancer, a Kaplan-Meier survival analysis was performed based on LINC00993 gene expression by using R2 database (http://r2.amc.nl). The overall survival curves suggested that breast cancer patients benefited from high LINC00993 expression (p < 0.0001; Figure 8G). The relapse-free survival analysis showed a similar result that a higher LINC00993 expression had a better outcome (p < 0.001; Figure 8H). Our results indicated that LINC00993 was strongly downregulated in TNBC, and its expression in breast cancer might serve as a prognostic biomarker. Taken all these results together, LINC00993 was breast-specific and was identified with a tumor suppressive feature and with prognostic value, suggesting that controlling LINC00993 level may be helpful for breast cancer treatment.

Discussion

Here we identified a novel breast-specific long non-coding RNA LINC00993 with tumor suppressive feature. We proved that LINC00993 was largely downregulated in breast cancer, especially in TNBC. LINC00993 suppressed breast cancer growth by causing G0/G1 arrest. We further showed that LINC00993 regulated key cell cycle-related genes such as p16INK4A, p14ARF, p21, and p53. As E2F played a significant role in regulation of G1 to S phase transition, LINC00993 could regulate translocation of E2F through two possible pathways, one is p16INK4A-CDK4/6- Rb/E2, the other is p14ARF- p53- p21-CDK2- Rb/E2 (Figure 6F). One interesting thing is that p16INK4A is a tumor suppressor protein encoded by the CDKN2A gene, while p14ARF is an alternate reading frame protein product of the CDKN2A locus. Because they were upregulated by LINC00993, future studies are needed to find out how LINC00993 affects p16INK4A, p14ARF expression. One possible mechanism is that LINC00993 might regulate CDKN2A gene transcription. Based on our findings, LINC00993 showed the potential to be a novel cancer therapeutic target.

Although LINC00993 was greatly downregulated in TNBC, how it was regulated remains unknown. By analyzing 862 pieces of breast-related ChIP-seq data sets from Cistrome database, FOXA1 binding peaks were found around LINC00993 gene in several data sets, indicating that FOXA1 might regulate LINC00993 expression (data not shown). FOXA1 could maintain luminal related genes and suppress basal-like genes (31). Expression of FOXA1 in luminal cancer was much higher than that in basal-like breast cancer. FOXA1 is the primary determinant of ER binding and transcriptional activity in breast cancer cells (32). This might explain why LINC00993 correlated with ER but did not affect each other (data not shown). They might be downstream elements of FOXA1. FOXA1 might serve as a key regulator of LINC00993, but how FOXA1 regulates LINC00993 remains to be investigated.

A large proportion of the human and mouse genome is transcribed as ncRNAs. Although some lncRNAs are verified to be functional, most of them remain largely unknown. Researchers have adopted guilt-by-association method to annotate lncRNAs (8). Although there are some database annotate ncRNAs this way, none of them offered information for LINC00993 (33). In this article, we combined two different bioinformatical method to annotate breast-specific LINC00993. First, we got a gene list containing genes mostly correlated with LINC00993 in expression. We did GO and KEGG analysis based on the gene list. Our results showed a strong relationship between LINC00993 and cell cycle process. We then used another method to annotate LINC00993, which we called single gene GSEA. Usually, GESA analysis needs at least two phenotypes to rank gene list. Here, we created the needed two groups by separating 1,096 breast cancer samples according to LINC00993 expression. To make the result more accurate, we set a threshold at 10%, the top 10% as LINC00993 high group and bottom 10% as low group. Our results showed that LINC00993 was involved in cell cycle pathway, and it might regulate G1/S transition (34). Those results were confirmed by wet lab experiments. Meanwhile, our GSEA analysis offered some other information about LINC00993, like chromosomal segregation, cilium pathway, and sperm formation. Considering that most of the lncRNAs are tissue-specific and expressed with a low abundance, annotation of them is a hard task. Based on RNA-Seq data from TCGA, “guilt by association” and GSEA method might be potential tools in exploring functions of cancer-related lncRNAs. However, more experiments were needed to evaluate the feasibility.

In conclusion, LINC00993 is a breast-specific lincRNA, which was downregulated in breast cancer, especially in TNBC. LINC00993 suppresses TNBC growth both in vitro and in vivo. Our findings here help with the understanding of TNBC development. LINC00993 showed the potential to be a novel cancer therapeutic target.

Materials and Methods

Cell Culture

The human breast cancer cell line MDA-MB-231 and BT-549 were purchased from the cell bank of Shanghai Institute of Biological Sciences, Chinese Academy of Science. These cells were maintained with DMEM basal medium containing 10% fetal bovine serum in a humidified incubator with 5% CO2 at 37°C.

RNA Isolation and Quantitative Real-Time PCR (qRT-PCR)

Total RNA was isolated using TRIzol reagent (Ambion, 15596-026, USA) according to manufacturer instructions. PrimeScript™ RT reagent Kit (Takara, RR037A) was used to synthesize the first-strand cDNA from an equal amount of the RNA sample. qRT-PCR was done by using SYBR® Premix Ex Taq® II (Takara, RR820A). The relative expression of mRNA was calculated with the 2−ΔΔT method. Each sample was repeated in triplicate. qRT-PCR data were analyzed and converted to relative fold changes. The following qRT-PCR assay primer sequences were used: LINC00993: 5′-AGTGCGGGGCTCATCTATCC (forward) and 5′TTACTGGGAGGTGGCAAGAG (reverse); Glyceraldehyde 3-phosphate dehydrogenase (GAPDH), 5′-CTCTGCTCCTCCTGTTCGAC (forward) and 5′ GCGCCCAATACGACCAAATC (reverse).

Plasmid Construction and Adenovirus Production and Infection

The shuttle plasmid of adenovirus was constructed by inserting synthesized LINC00993 sequence into pHBAd-MCMV-RFP vector (Hanbio, Shanghai). The construction was done by using Gibson Assembly method with the restriction enzyme sites Not I and Nsi I. The shuttle plasmid was confirmed by sequencing. The shuttle plasmids were transfected into HEK293 cells together with the adenovirus packaging plasmid pBHGloxdelE13cre using the Lipofectamine® 2000 DNA Transfection Reagent (Thermo Fisher Scientific). Six days later, first-generation adenoviruses were harvested. All adenoviruses were amplified in HEK293 cells and purified by ultracentrifugation on cesium chloride (CsCl) gradients. The titer of adenovirus was determined by TCID50 method. For MDA-MB-231 cell infection, the multiplicity of infection (MOI) equals to 40 was adopted after multiple tests.

Flow Cytometry Analysis

Flow cytometry analysis was used to analyze cell cycle. In brief, after transfection with adenovirus vectors for 48 h, cells were trypsinized and washed twice with ice-cold PBS, and then fixed with 75% ethanol at 4°C overnight. The fixed cells were suspended in propridium iodide (PI, 100 mg/ml) or RNase (10 mg/ml) or PBS, and incubated at 37°C for 30 min in the dark. All samples were assessed with FACScan system (BD Biosciences, USA). Data were analyzed by Cell Quest software (BD Biosciences, USA). Measurements were repeated independently three times.

Transwell Assay

Cell invasion capability was measured using Corning BioCoat™ 24-well-cell permeable supports (8 μm pore size, Corning, USA). Briefly, the upper chamber of the Transwell was coated with 150 μl diluted BD Matrigel™ Basement membrane Matrix (BD, USA) gel solution (200 μg/ml) at 37°C for 3 h. Medium (800 μl) with 10% FBS was added to the lower chamber, and 2 × 104 cells were then plated in the upper chamber with 100 μl serum-free medium. After incubation at 37°C with 5% CO2 for 24 h, cells were fixed with methanol and acetic acid (3:1) and stained with Giemsa staining (SolarBio, Beijing, China). Cells remaining on the upper membrane were gently removed by wiping with cotton swabs. Cell invasion was assessed by counting stained cells under a microscope on five random fields. Three replicates were performed.

TUNEL Assays

To analyze cell apoptosis, terminal deoxynucleotidyl transferase-mediated dUTP nick-end labeling (TUNEL) assays were performed with One Step TUNEL Apoptosis Assay Kit (FITC, Beyotime, China) according to the manufacturer's instructions. DNase I-treated cells were used as positive control. FITC-labeled apoptotic cells were observed under a fluorescence microscope.

CCK-8 Assay, Cell Growth Curve, and Colony Formation Assays

Cell viability was measured using Cell Counting Kit-8 (CCK-8) (Dojindo, Japan). For cell growth assay, cells were plated in a 24-cell plate, and after adenovirus infection, three wells of cells were randomly picked up, digested, and counted. For the colony formation assay, 1–3 × 103 cells were plated in each well of a six-well plate and incubated at 37°C for 1–2 weeks. The cells were fixed with 4% paraformaldehyde and stained with 1% crystal violet (Sigma-Aldrich). Cell colonies were counted and analyzed.

In vivo Tumor Model

The animal experimental protocols were approved by the Committee on Ethical Use of Animal of the First Affiliated Hospital of Chongqing Medical University. Female BALB/c nude mice of 4~6 weeks old were maintained under SPF (specific pathogen-free) conditions in the Experimental Animal Department of the Chongqing Medical University. The mice (n = 21) were randomized and grouped with seven each. MDA-MB-231 cells (2 × 106) in 200 μl serum-free medium were subcutaneously injected to the right dorsal flank of the mice. When ~100 mm3 tumors were observed, the mice were received intratumoral injection with 109 plaque-forming units (PFU) Ad-LINC00993, Ad-RFP, or PBS in a total volume of 50 μl. Delivery of viral vectors or PBS started at the second week and was repeated every 2 weeks. Five weeks after cell injection, the mice were sacrificed. The total volumes of tumors were assessed twice a week using a caliper and was calculated with the longest diameter (X) as well as the shortest diameter (Y) using the formula: X * Y2/2. The weight of tumor was measured after mice were sacrificed.

RNA Pulldown and Mass Spectrometry

In vitro, biotin-labeled RNAs were transcribed with the Biotin RNA Labeling Mix (Roche) and T7 RNA polymerase, treated with 2 μl RNase-free DNase I at 37°C for 15 min to remove DNA and add 2 μl 0.2 M EDTA (pH = 8.0) to stop the reaction. Biotinylated RNA (1 μg) in RNA structure buffer (10 mM Tris pH 7, 0.1 M KCl, 10 mM MgCl2) was heated to 95°C for 2 min, on ice for 3 min, and then left at room temperature (RT) for 30 min to allow proper secondary structure formation. Folded RNA was then mixed with human MDA-MB-231 cell extract (containing 1 mg proteins) in 500 μl RIP buffer and then incubated at RT for 1 h. Washed streptavidin agarose beads (50 μl) (Invitrogen, Cat. No. SA10004) was added to each binding reaction and further incubated at RT for another 1 h. Beads were washed briefly with RIP buffer three times and boiled in SDS buffer. Then, the retrieved proteins were detected by Western blot or resolved in gradient gel electrophoresis followed by mass spectrometry (MS) identification.

Western Blot

Cell lysates were separated by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). Autoradiograms were quantified by densitometry (Quantity One software; Bio-Rad). GAPDH antibody was used as control. p16INK4A, p14ARF, p21, p53, and GAPDH antibody were purchased from Proteintech (Wuhan, China).

Bioinformatics Methods and Data Analysis

Correlation analysis was done by using GEPIA using Spearman method (35). Survival analysis was carried out by using R2 (http://r2.amc.nl), a Genomics Analysis and Visualization Platform. BRCA data were downloaded from TCGA by using GDC Data Transfer Tool (36). Those data were arranged by using R language and were and normalized by “DESeq2” package (37). GSEA analysis requires at least two different phenotypes. One thousand ninety-six BRCA RNA-Seq data were downloaded from TCGA by using GDC tools. Those samples were arranged according to LINC00993 expression. The top 10% (110) of the data that expressed a high level of LINC00993 were chosen as high group; the bottom 10% (110) as low group. GSEA analysis was performed by using the two groups. GSEA analysis was done by using the R script from Board Company. The statistical significance of comparisons between two groups was determined with Student's t-test. P < 0.05 were considered statistically significant.

Patients and Tumor Specimens

The study was approved by the Ethics Committee of Chongqing Medical University, and it was performed in compliance with the Declaration of Helsinki principles. The written informed consent was obtained from all patients. Breast cancer samples from 98 patients consecutively from July 2015 until Jun 2018 were analyzed in this study. Those patients underwent surgical resection of the breast cancer at the First Affiliated Hospital of Chongqing Medical University. Tumor specimens were from surgical resection and frozen in liquid nitrogen for 1 h within 30 min and then stored in −80°C for long time preservation.

Data Analysis

The values represent the mean ± SEM. The statistical significance of comparisons between two groups was determined with Student's t-test. P < 0.05 were considered statistically significant.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/; https://gtexportal.org/; http://gepia.cancer-pku.cn/index.html; https://hgserver1.amc.nl/cgi-bin/r2/main.cgi.

Ethics Statement

The studies involving human participants were reviewed and approved by Human Research Ethics Committee of the First Affiliated Hospital of Chongqing Medical University. The patients/participants provided their written informed consent to participate in this study. This animal study was reviewed and approved by Animal Experimentation Ethics Committee of Chongqing Medical University.

Author Contributions

SG, LJ, KT, and CC performed the research. SL designed the research study. SG contributed essential reagents or tools. SG, LJ, HY, and CC analyzed the data. SG and SL wrote the paper. All authors have read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (NSFC) Grant 81772979 and 81472658 (SL).

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2019.01325/full#supplementary-material

References

1. Perou CM. Molecular stratification of triple-negative breast cancers. Oncologist. (2011) 16(Suppl. 1):61–70. doi: 10.1634/theoncologist.2011-S1-61

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Foulkes WD, Smith IE, Reis-Filho JS. Triple-negative breast cancer. N Engl J Med. (2010) 363:1938–48. doi: 10.1056/NEJMra1001389

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Dent R, Trudeau M, Pritchard KI, Hanna WM, Kahn HK, Sawka CA, et al. Triple-negative breast cancer: clinical features and patterns of recurrence. Clin Cancer Res. (2007) 13(15 Pt 1):4429–34. doi: 10.1158/1078-0432.CCR-06-3045

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Djebali S, Davis CA, Merkel A, Dobin A, Lassmann T, Mortazavi A, et al. Landscape of transcription in human cells. Nature. (2012) 489:101–8. doi: 10.1038/nature11233

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kopp F, Mendell JT. Functional classification and experimental dissection of long noncoding RNAs. Cell. (2018) 172:393–407. doi: 10.1016/j.cell.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lin C, Yang L. Long noncoding RNA in cancer: wiring signaling circuitry. Trends Cell Biol. (2018) 28:287–301. doi: 10.1016/j.tcb.2017.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Van Grembergen O, Bizet M, de Bony EJ, Calonne E. Portraying breast cancers with long noncoding RNAs. Sci Adv. (2016) 2:e1600220. doi: 10.1126/sciadv.1600220

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lin A, Li C, Xing Z, Hu Q, Liang K, Han L, et al. The LINK-A lncRNA activates normoxic HIF1alpha signalling in triple-negative breast cancer. Nat Cell Biol. (2016) 18:213–24. doi: 10.1038/ncb3295

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liu B, Sun L, Liu Q, Gong C, Yao Y, Lv X, et al. A cytoplasmic NF-kappaB interacting long noncoding RNA blocks IkappaB phosphorylation and suppresses breast cancer metastasis. Cancer Cell. (2015) 27:370–81. doi: 10.1016/j.ccell.2015.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Jin C, Yan B, Lu Q, Lin Y, Ma L. Reciprocal regulation of Hsa-miR-1 and long noncoding RNA MALAT1 promotes triple-negative breast cancer development. Tumour Biol. (2016) 37:7383–94. doi: 10.1007/s13277-015-4605-6

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Arun G, Diermeier S, Akerman M, Chang KC, Wilkinson JE, Hearn S, et al. Differentiation of mammary tumors and reduction in metastasis upon Malat1 lncRNA loss. Genes Dev. (2016) 30:34–51. doi: 10.1101/gad.270959.115

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Arun G, Diermeier SD, Spector DL. Therapeutic targeting of long non-coding RNAs in cancer. Trends Mol Med. (2018) 24:257–77. doi: 10.1016/j.molmed.2018.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wu Y, Shao A, Wang L, Hu K, Yu C, Pan C, et al. The role of lncRNAs in the distant metastasis of breast cancer. Front Oncol. (2019) 9:407. doi: 10.3389/fonc.2019.00407

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chen C, Li Z, Yang Y, Xiang T, Song W, Liu S. Microarray expression profiling of dysregulated long non-coding RNAs in triple-negative breast cancer. Cancer Biol Ther. (2015) 16:856–65. doi: 10.1080/15384047.2015.1040957

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Jager D, Filonenko V, Gout I, Frosina D, Eastlake-Wade S, Castelli S, et al. NY-BR-1 is a differentiation antigen of the mammary gland. Appl Immunohistochem Mol Morphol. (2007) 15:77–83. doi: 10.1097/01.pai.0000213111.05108.a0

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Carithers LJ, Ardlie K, Barcus M, Branton PA, Britton A, Buia SA, et al. A novel approach to high-quality postmortem tissue procurement: the GTEx project. Biopreserv Biobank. (2015) 13:311–9. doi: 10.1089/bio.2015.0032

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Rinn JL, Chang HY. Genome regulation by long noncoding RNAs. Ann Rev Biochem. (2012) 81:145–66. doi: 10.1146/annurev-biochem-051410-092902

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Lefever S, Anckaert J, Volders PJ, Luypaert M, Vandesompele J, Mestdagh P. decodeRNA- predicting non-coding RNA functions using guilt-by-association. Database. (2017) 2017:bax042. doi: 10.1093/database/bax042

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. (2000) 25:25–9. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhu D, Xu G, Ghandhi S, Hubbard K. Modulation of the expression of p16INK4a and p14ARF by hnRNP A1 and A2 RNA binding proteins: implications for cellular senescence. J Cell Physiol. (2002) 193:19–25. doi: 10.1002/jcp.10147

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Morris L, Allen KE, La Thangue NB. Regulation of E2F transcription by cyclin E-Cdk2 kinase mediated through p300/CBP co-activators. Nat Cell Biol. (2000) 2:232–9. doi: 10.1038/35008660

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Giono LE, Manfredi JJ. Mdm2 is required for inhibition of Cdk2 activity by p21, thereby contributing to p53-dependent cell cycle arrest. Mol Cell Biol. (2007) 27:4166–78. doi: 10.1128/MCB.01967-06

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Wei W, Hemmer RM, Sedivy JM. Role of p14(ARF) in replicative and induced senescence of human fibroblasts. Mol Cell Biol. (2001) 21:6748–57. doi: 10.1128/MCB.21.20.6748-6757.2001

PubMed Abstract | CrossRef Full Text | Google Scholar

26. DeGregori J. The Rb network. J Cell Sci. (2004) 117(Pt 16):3411–3. doi: 10.1242/jcs.01189

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Huarte M, Guttman M, Feldser D, Garber M, Koziol MJ, Kenzelmann-Broz D, et al. A large intergenic noncoding RNA induced by p53 mediates global gene repression in the p53 response. Cell. (2010) 142:409–19. doi: 10.1016/j.cell.2010.06.040

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Carpenter S, Aiello D, Atianand MK, Ricci EP, Gandhi P, Hall LL, et al. A long noncoding RNA mediates both activation and repression of immune response genes. Science. (2013) 341:789–92. doi: 10.1126/science.1240925

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Kumar PP, Emechebe U, Smith R, Franklin S, Moore B, Yandell M, et al. Coordinated control of senescence by lncRNA and a novel T-box3 co-repressor complex. Elife. (2014) 3:e02805. doi: 10.7554/eLife.02805

CrossRef Full Text | Google Scholar

30. Parker JS, Mullins M, Cheang MC, Leung S, Voduc D, Vickery T, et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J Clin Oncol. (2009) 27:1160–7. doi: 10.1200/JCO.2008.18.1370

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Bernardo GM, Bebek G, Ginther CL, Sizemore ST, Lozada KL, Miedler JD, et al. FOXA1 represses the molecular phenotype of basal breast cancer cells. Oncogene. (2013) 32:554–63. doi: 10.1038/onc.2012.62

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Hurtado A, Holmes KA, Ross-Innes CS, Schmidt D, Carroll JS. FOXA1 is a key determinant of estrogen receptor function and endocrine response. Nat Genet. (2011) 43:27–33. doi: 10.1038/ng.730

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Theurillat JP, Zurrer-Hardi U, Varga Z, Storz M, Probst-Hensch NM, Seifert B, et al. NY-BR-1 protein expression in breast carcinoma: a mammary gland differentiation antigen as target for cancer immunotherapy. Cancer Immunol Immunother. (2007) 56:1723–31. doi: 10.1007/s00262-007-0316-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Lohrum MA, Ashcroft M, Kubbutat MH, Vousden KH. Contribution of two independent MDM2-binding domains in p14(ARF) to p53 stabilization. Curr Biol. (2000) 10:539–42. doi: 10.1016/S0960-9822(00)00472-3

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. (2017) 45:W98–102. doi: 10.1093/nar/gkx247

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Grossman RL, Heath AP, Ferretti V, Varmus HE, Lowy DR, Kibbe WA, et al. Toward a shared vision for cancer genomic data. N Engl J Med. (2016) 375:1109–12. doi: 10.1056/NEJMp1607591

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: LINC00993, LncRNA, TNBC, cell cycle, guilt by association

Citation: Guo S, Jian L, Tao K, Chen C, Yu H and Liu S (2019) Novel Breast-Specific Long Non-coding RNA LINC00993 Acts as a Tumor Suppressor in Triple-Negative Breast Cancer. Front. Oncol. 9:1325. doi: 10.3389/fonc.2019.01325

Received: 11 July 2019; Accepted: 13 November 2019;
Published: 17 December 2019.

Edited by:

Jian-Guo Zhou, University of Erlangen Nuremberg, Germany

Reviewed by:

Jiang-Bo Liu, First Affiliated Hospital of Henan University of Science and Technology, China
Liqun Zou, Sichuan University, China

Copyright © 2019 Guo, Jian, Tao, Chen, Yu and Liu. 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: Shengchun Liu, liushengchun1968@163.com

These authors have contributed equally to this work