Integrated Whole Transcriptome Profiling and Bioinformatics Analysis for Revealing Regulatory Pathways Associated With Quercetin-Induced Apoptosis in HCT-116 Cells

Quercetin (QUE) is a bioactive component that belongs to the natural flavonoids group, and recent researchers found that it could prevent colorectal cancer (CRC). However, the exact mechanism by which QUE exerts its anti-tumor effects in CRC remains unclear. In this study, MTS assay and flow cytometry were used to detect the anti-tumor effects of QUE on HCT-116 cells. The results showed that QUE could inhibit the proliferation and induce apoptosis of HCT-116 cells. Furthermore, whole transcriptome sequencing was employed to establish the microRNA (miRNA), long non-coding RNA (lncRNA), circular RNA (circRNA), and mRNA profiles. A total of 240 differentially expressed lncRNAs (DElncRNAs), 131 circRNAs (DEcircRNAs), 83 miRNAs (DEmiRNAs), and 1415 mRNAs (DEmRNAs) were identified in the QUE-treated HCT-116 cells compared to the untreated HCT-116 cells. Then, quantitative real-time polymerase chain reaction (qRT-PCR) was used to validate the expression of selected circRNAs, miRNAs, lncRNAs, and mRNAs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed to further investigate RNAs’ biological functions and potential mechanisms. Based on the theory of competing endogenous RNA (ceRNA), the circRNA–miRNA–mRNA and lncRNA–miRNA–mRNA regulatory networks were constructed to illustrate the regulatory relationship between non-coding RNA (ncRNA) and mRNA. Our results provided novel information about the molecular basis of QUE in treating CRC. Our findings indicated that deep RNA sequencing analysis of mRNA and ncRNAs was a promising approach to research anticancer mechanisms.


INTRODUCTION
Colorectal cancer (CRC) is one of the most prevalent malignancies, ranking as the second leading cause of cancerinduced death worldwide (Bray et al., 2018). According to global statistics, it is estimated that CRC disease will affect more than 2.2 million new patients and will cause 1.1 million cancer deaths by 2030 (Arnold et al., 2017). In the United States, the number of new cases of CRC reached 145,600, and the number of deaths reached 51,020 in 2019 (Siegel et al., 2019). In recent years, despite advances in diagnosis and treatment of CRC, any of the current therapies cannot offer an effective therapeutic outcome due to the poor prognosis and high recurrence rate. The mortality rates are still high and unacceptable (Hamilton, 2018;Peng et al., 2018). Therefore, there is an urgent need for new CRC therapies.
Quercetin (QUE), a bioactive component belonging to the natural flavonoids group, is distributed in fresh onions, fruits, and vegetables (Russo et al., 2014). Previous studies suggested a positive association between QUE intake and improved outcomes of several types of cancer treatment (Hashemzaei et al., 2017). Recently, intense attention has been paid to the chemopreventive and anti-tumor functions of QUE on colon cancer and increasing experimental evidence verified these effects in vivo and in vitro. Moreover, extensive studies had reported that QUE has anti-cell proliferation, pro-apoptosis, anti-angiogenesis, and anti-metastasis effects on colon cancer (Darband et al., 2018). In addition, Maria et al. confirmed that the CB1 receptor mediated the anti-proliferative and pro-apoptotic of QUE effects in human colon cancer cell lines (Refolo et al., 2015). A transcriptomics and proteomics research was conducted to characterize gene and protein changes occurring in the distal colon mucosa of rats supplemented with QUE; this further enriches the anticancer mechanism of QUE in vitro (Dihal et al., 2008). Although multiple lines of studies demonstrated the efficacy of QUE, the exact mechanism of its anti-tumor effects in CRC remains unclear. Hence, more advanced tools with large-scale data are required for the deep interpretation of the problems.
Non-coding RNAs (ncRNAs) are an abundant class of RNAs that typically do not encode proteins but functionally regulate protein expression (Mattick and Makunin, 2006). As part of RNA-protein complexes in regulating gene expression, ncRNAs can be subdivided into several families based on their size and biogenesis pathways, including microRNAs (miRNAs, with <200 nucleotides), long ncRNAs (lncRNAs, with a length >200 bp), and circular RNAs (circRNAs, with a closed continuous loop) (Cech and Steitz, 2014;Thum, 2014;Thum and Condorelli, 2015). These different ncRNAs are further classified based on sequence or structure conservation, subcellular localization and function, and association with annotated protein-coding genes and other DNA elements of known function (St Laurent, 2015). Accumulating evidence indicated that ncRNAs are involved in a remarkable variety of biological functions. Moreover, multiple lines of evidence have linked ncRNA mutations and dysregulation with various human diseases, especially cancers. Colon cancer is characterized by genetic and epigenetic modifications; thus, ncRNAs were emerging key regulators of gene expression under colon cancer (Hu et al., 2014;Zhou et al., 2016;Li et al., 2018). Hamfjord et al. (2012) profiled at least 37 differentially expressed miRNAs between CRC tissues and normal colon mucosa. Hofsli et al. (2013) detected that 21 miRNAs differentially expressed in serum of CRC patients may be applied to identify colon cancer patients at an early stage of the disease. The expression profile of circRNA in CRC tissues has also been identified, which could be used as new biomarkers for prognosis and diagnosis of CRC (Bachmayr-Heyda et al., 2015). A recent report has summarized that the lncRNA regulating cancer cell proliferation, migration, and apoptosis might also serve as biomarkers for CRC diagnosis and prognosis, including LOC285194, RP11-462C24.1, BANCR, NR_034119, NR_029373, lncRNA-AFAP1-AS1, and NR_026817.79 (Deng et al., 2017).
However, few of transcriptomes concerned about the transcriptomic effect of QUE intervention against CRC were published, and various questions remain unanswered with respect to the ncRNAs in response to QUE treatment. Thus, the whole transcriptome sequencing technique was used to profile the coding transcriptome and ncRNAs changes occurred in CRC cell response to QUE treatment. Additionally, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted to explore the biological roles and potential signaling pathways of these differentially expressed mRNA and ncRNAs. Moreover, co-expression networks were constructed based on these sequencing data and bioinformatics analysis. These results could provide novel insights into the molecular basis of QUE in treating CRC.

Cell Lines and Cell Culture
The human colon cancer cell line (HCT-116) was purchased from Shanghai Institutes for Biological Sciences (SIBS, Shanghai, China). The cells were cultured in Dulbecco's modified Eagle's medium contained in tissue culture flasks. The medium was supplemented with 10% fetal bovine serum (FBS), antibiotics, and 1% glutamine. The cells were cultured at 37°C under a humidified atmosphere containing 5% CO 2 and 95% air.

Cell Apoptosis Detection
In brief, HCT-116 cells were seeded in six-well plates at a density of 4 × 10 4 /well and incubated with QUE at different concentrations (100, 150, and 200 μΜ) for 72 h. Then, cells were collected by centrifugation, washed with phosphate-buffered saline, and resuspended with 400 µl of buffer. After that, cell suspensions were stained with Annexin V-FITC and propidium iodide (PI) according to the manufacturer's instruction. At last, cell apoptosis detection was performed by a flow cytometer (Shanghai Pudi Biotechnology Co., Ltd., BDFACS Calibur) and analyzed by FlowJo 7.6.

RNA Extraction and Quality Monitoring
Transcriptomic analysis was performed on untreated HCT-116 cells and HCT-116 cells treated with 150 μΜ QUE for 48 h. At first, total RNA of HCT-116 cells was extracted from three independent experiments by Trizol reagent (Invitrogen, Carlsbad, CA, USA) in accordance with the manufacturer's protocol. RNA quantity and quality were measured by a spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA) and the Agilent 2100 Bioanalyzer (Agilent Technologies).

RNA Library Construction, Quality Control, and Sequencing
The small RNA library (sRNA library) was constructed according to the manufacturer's instruction of TrueSeq small RNA library prep kit (Illumina San Diego CA, USA), and 2.5 ng RNA per sample was used as the initial amount. The T4 RNA Ligase 1 was ligated to the 3′ end of the RNA, followed by the ligation of T4 RNA Ligase 2 (truncated) to the 5′ adapter. After that, the RNA was reverse transcribed to synthesize cDNA. Finally, small RNA libraries were obtained by screening these recovered fragments using gel separation technology. In the construction of lncRNA library (chain-specific library for removal of ribosomal RNA), the epicenter Ribo-ZeroTM kit was used to remove the ribosomal RNA. Then, rRNA-depleted RNA was fragmented and used as a template to construct the cDNA library. The qualities of the libraries were further tested following these steps: 1) Initial quantification was performed using Qubit 2.0, and the insert size of the library was tested using Agilent 2100. 2) The Q-PCR method was used to accurately quantify the effective concentration of the library (effective library concentration > 2 nM). After the libraries' quality tests were passed, different libraries performed pooling according to the amount of target data, followed by sequencing on the Illumina HiSeq platform.

Data Processing
Based on Sequencing-By-Synthesis (SBS) technology, the Illumina HiSeq high-throughput sequencing platform sequenced cDNA libraries to produce a large amount of high-quality data, i.e., raw data. In order to ensure the accuracy of information analysis, quality control of the original data is required to obtain a highquality sequence, which is clean reads. The exclusion criteria for quality control of the original sequence measured by the deribosomal library included the following: 1) low-quality data, 2) reads containing jointed-sequence, and 3) reads containing N (undetermined base information) ratio greater than 5%. The exclusion criteria for quality control of the original sequence measured by the sRNA library included the following: 1) lowquality data, 2) reads containing N ratio greater than 10%, 3) reads without 3′ linker sequence, and 4) sequences shorter than 15 or longer than 35 nucleotides.
The targets for miRNAs were predicted using miRanda (Betel et al., 2008) and targetscan (Lewis et al., 2003). The targets for lncRNAs were predicted based on position relationship (within 100 kb of lncRNA) and Pearson correlation coefficient (when sample size >5, the absolute value of correlation >0.9, and p value < 0.01) between lncRNA and mRNA. The source genes for circRNAs were acquired according to the position of the circRNA sequence; namely, the gene on which the circRNA sequence was located was deemed as its source gene.

Quantitative Real-Time PCR Validation
To validate the expression levels of the selected lncRNAs, miRNAs, circRNAs, and mRNAs by quantitative real-time polymerase chain reaction (qRT-PCR), RNA samples from the untreated HCT-116 cells and HCT-116 cells treated with 150 μΜ QUE for 48 h were collected. Total RNA was isolated by using the Trizol reagent (Invitrogen, USA), and then reverse-transcribed into cDNA using SuperScript III Reverse Transcriptase (Invitrogen, USA) according to the manufacturer's instruction. A Gene Amp PCR System 9700 (Applied Biosystems, USA) and 2× PCR Master Mix (Arraystar, USA) were used to perform qRT-PCR in accordance with the manufacturer's instructions. The circRNA, lncRNA, and mRNA expression levels were normalized to β-actin (Human), an endogenous reference transcript, and calculated using the standard curve-based method (Larionov et al., 2005). The expression levels of the selected miRNAs, being normalized against the U6 (human), were calculated by the 2 −ΔΔCt method using untreated control as the calibrator (Livak and Schmittgen, 2001). The sequences of all primers are shown in Table 1. The data represent the means of three experiments.

GO and KEGG Pathway Analysis
To better understand the biological functions and potential mechanisms of ncRNAs and mRNAs in the mechanism of QUE acting on CRC, we applied GO enrichment and KEGG pathway analyses on these differentially expressed mRNAs (DEmRNAs), source genes of differentially expressed circRNAs (DEcircRNAs), predicted target genes of differentially expressed miRNAs (DEmiRNAs), and predicted target genes of differentially expressed lncRNAs (DElncRNAs). Briefly, GO analyses (www. geneontology.org) consisted of three components: biological process (BP), cellular component (CC), and molecular function (MF). KEGG analyses were carried to investigate the potential significant pathways (http://www.genome.jp/kegg/).

Construction of Co-expression Network
CircRNA, lncRNA, and mRNA could inhibit target gene regulation of miRNA to indirectly regulate gene expression, serving as a miRNA sponge to bind miRNA competitively with its binding sites, which was called competing endogenous RNA (ceRNA). Based on the theory of ceRNA, the circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA network were constructed and visually displayed using the Cystoscope software V3.5.0 (San Diego, CA, USA), as previous described. Different shapes and colors represent different RNA types and regulated relationships, respectively.

Statistical Analysis
Statistical differences were determined by using SPSS 22.0 software. A p value less than 0.05 was considered to be significant. For sequencing data, we analyzed DE ncRNAs and DEmRNAs using DEseq (Anders and Huber, 2010) software, with p < 0.05 and |log2 (fold change)| > 2 as screening criteria.

QUE Inhibits the Proliferation of HCT116 Cells
As shown in Figure 1A, QUE at concentrations of 100, 150, 200, and 250 μΜ was used to treat HCT-116 cells for 24, 48, and 72 h, respectively, and the results showed that the viability of HCT-116 cells was markedly reduced by QUE, in both a dose-and timedependent manner. As shown in Figure 1B, the IC 50 values for QUE were 220.1, 179.7, and 143.23 μΜ.

QUE Induces Apoptosis in HCT-116 Cells
Cells were treated with 100, 150, and 200 μM QUE for 72 h and stained with Annexin V-FITC and PI for the analysis of apoptosis. As shown in Figure 2, the percentage of apoptotic cells was significantly increased in a dose-dependent manner after QUE treatment. Taken together, these results clearly demonstrated QUE-induced apoptosis of HCT-116 cells.

Sequencing and Mapping of the Transcriptome
The cDNA and sRNA libraries of cell samples from three groups of untreated HCT-116 cells and three groups of QUE-treated HCT-116 cells were sequenced. Moreover, counts of clean reads and mapped ratio of sequencing data are shown in Table 2.

Differentially Expressed ncRNAs and mRNAs
Information of the top 20 up-regulated and 20 down-regulated lncRNAs, circRNAs, miRNAs, and mRNAs in the QUE-treated cells compared with the untreated cells are listed in Tables 3-6 In total, our project detected 15,755 lncRNAs, 25,843 mRNAs, 2,735 miRNAs, and 6,911 new circRNAs. As shown in Figure  4A, there were 240 DElncRNAs (89 down-regulation and 151 up-regulation), 131 DEcircRNAs (37 down-regulation and 94 up-regulation), 83 DEmiRNAs (19 down-regulation and 64 up-regulation), and 1,415 DEmRNAs (901 down-regulation and 514 up-regulation) in the QUE-treated cells compared to the untreated cells, respectively. In an effort to further pinpoint the genes involved in the mechanism of QUE anti-CRC, we performed a global overlapping gene analysis. The intersections of DEmRNA, DEmiRNA-Target mRNA, DELncRNA-Target Trans.mRNA, and DEcircRNA-Hostgene are showed in Figure 4B. This left us with six reliable core mRNAs validated by different platforms and across different institutions, including ENSG00000101187, ENSG00000111344, ENSG00000116183, ENSG00000125319, ENSG00000150637, and ENSG00000169035.

GO Enrichment Analysis
Based on the GO enrichment analysis of the cis targeted genes of lncRNA (Figure 6A), the most significantly enriched BP, CC, and MF were protein ADP-ribosylation, actin cytoskeleton, and RNA-directed DNA polymerase activity, respectively. Based on the GO enrichment analysis of the trans targeted genes of lncRNA ( Figure 6B), the most significantly enriched BP, CC, and MF were stimulatory C-type lectin receptor signaling pathway, spindle pole, and protein N-terminus binding, respectively. Based on the GO enrichment analysis of DEmRNAs (Figure 6C), the most significantly enriched BP, CC, and MF were defense response to virus, apical plasma membrane, and structural molecule activity, respectively. Based on the GO enrichment analysis of source genes of DEcircRNAs (Figure 6D), the most significantly enriched BP, CC, and MF were protein phosphorylation, nuclear matrix, and phospholipid binding, respectively. Based on the GO enrichment analysis of targeted genes of DEmiRNAs (Figure 6E), the most significantly enriched BP, CC, and MF were positive regulation of apoptotic process, cell junction, and phospholipid binding, respectively.

KEGG Pathway Enrichment Analysis
The most significantly enriched KEGG pathways are shown in Figure 7. For the cis targeted genes of DElncRNA ( Figure 7A), alcoholism, systemic lupus erythematosus, and antigen processing and presentation were the most significant pathways for enrichment. For the trans targeted genes of DElncRNA ( Figure 7B), cell cycle, pyrimidine metabolism, and glycolysis/gluconeogenesis were the most significant enriched pathway. For the DEmRNAs (Figure 7C), Alzheimer's disease, carbon metabolism, and biosynthesis of amino acids were the most significant enriched pathway. For the hostgenes of DEcircRNA (Figure 7D), Fanconi anemia pathway, glycosaminoglycan biosynthesis-keratan sulfate, and DNA replication were the most significant enriched pathway. For the DE miRNAs ( Figure 7E), focal adhesion, cAMP signaling pathway, and miRNAs in cancer were the most significant enriched pathway. The main biochemical pathways and signal transduction pathways determined by KEGG analysis will provide further insight into future research directions of ncRNAs and mRNA.

Regulatory Network of ncRNAs and mRNA
To explore the molecular mechanism of ncRNAs, a circRNA-miRNA-mRNA regulatory network and lncRNA-miRNA-mRNA   network were constructed based on the theory of ceRNA. By using circRNA as a decoy, miRNA as the center, and mRNA as the target, the circRNA-miRNA-mRNA regulation network containing 437 mRNAs, 23 circRNAs, and 17 miRNAs was generated (Figure 8). By using lncRNA as a decoy, miRNA as center, and mRNA as target, the lncRNA-miRNA-mRNA regulation network that contained 331 mRNAs, 24 lncRNAs, and 13 miRNAs was built (Figure 9). In these two networks, different shapes represent different RNA types; red and green represent up-and down-regulation, respectively. Interestingly, circRNA (8:93786223|93822563), ENST00000313807, and ENST00000449307 were all able to regulate LRG1 expression, through competitively binding with miR-5096. These results suggested that circRNAs and lncRNAs harbor miRNA response elements and play pivotal regulatory roles in the mechanisms of QUE in anti-CRC.

DISCUSSION
CRC is one of the most common types of malignant tumors, ranking as the second leading cause of cancer-induced death worldwide. Multiple lines of studies have demonstrated the efficacy of QUE in CRC (Dihal et al., 2008;Refolo et al., 2015;Darband et al., 2018), but the exact mechanism of its anti-tumor effects in CRC remains unclear. To the best of our knowledge, this is the first comprehensive report of lncRNA, mRNA, circRNA, and miRNA to reveal regulator pathways with regard to QUEinduced apoptosis in HCT-116 cells.
In general, our data confirmed that QUE could inhibit the proliferation and induce apoptosis of HCT-116 cells. With FC ≥ 2.0 and p value < 0.05 thresholds, a total of 240 lncRNAs, 131 circRNAs, 83 miRNAs, and 1415 mRNAs with significant differential expression were identified in the QUE-treated cells compared with the untreated cells. We found that several DEmRNAs and DEmiRNAs may be connected with CRC. However, the great majority of DElncRNAs and DEcircRNAs were not known, mainly due to rarely related research. In addition, 12 dysregulated ncRNAs and mRNAs identified were selected for qRT-PCR validation, and the results confirmed the sequencing analysis findings to some extent. On the basis of the KEGG analysis, two significantly enriched pathways, PI3K-Akt and Ras signaling pathway, were participated by all four RNAs. Previous studies have reported that these two cancer-related pathways may be the possible molecular mechanisms of QUE in anti-CRC (Kim et al., 2005;Psahoulia et al., 2007;Xavier et al., 2009;Yang et al., 2016). Our research results further confirmed this possibility, more deeply, comprehensively, and systematically. Firstly, we focused on the differentially expressed coding genes. The top 20 DEmRNAs were regarded as the most important ones involved in mechanism of QUE acting on CRC. Among them, pioneering studies demonstrated that their dysregulation may result in the progression of CRC, such as AZGP1 (Ji et al., 2013;Chang et al., 2014;Xue et al., 2014b), APOBEC3G (Ding et al., 2011;Lan et al., 2014), BST2 (Mukai et al., 2017), TRIM29 (Jiang et al., 2013;Xu et al., 2016), and S100A4 (Dahlmann et al., 2014;Fei, 2017). Our results first suggested that MAGEA6, TNFSF12-TNFSF13, ANK1, PRR15L, BAG6, CD79B, KRT20, CTD-2410N18.5, RP11-385D13.1, HLA-DRB1, TTYH1, TNFRSF6B, PPP1R1B, SMIM5, and PTGDS may also exert their functions. The KEGG pathway analysis indicated that QUE-responsive gene alterations in HCT-116 cells were significantly enriched in the MAPK signaling pathway, a widely known cancer-related pathway. LRG1, which enriched in the MAPK pathway, has been strongly associated with worse overall survival for CRC, and may be considered as an independent prognostic indicator for CRC (Ladd et al., 2012;Zhang et al., 2016;Zhou et al., 2017;Zhang et al., 2018). In this study, LRG1 was proved to be down-regulated in QUE-treated HCT-116 cells compared to the untreated HCT-116 cells. These investigations indicated that LRG1 might be a potential drug target of QUE acting on CRC by regulating the MAPK signaling pathway. Subsequently, the effects of QUE on ncRNAss, including miRNAs, lncRNAs, and circRNAs, were evaluated in the current study. Previous studies indicated that miR-338-3p dysregulation may contribute to the progression of CRC (Xue et al., 2014a). In our study, miR-338-3p was significantly up-regulated in HCT-116 cells treated with QUE, and was the top regulated ones. We suggested that miR-338-3p may be an important one participating in QUE's anti-CRC mechanism. Additionally, four DEmiRNAs were referred as the most likely candidate miRNA associated with the mechanism of QUE. Among these, the expression of miR-320b, miR-320c, and miR-320d was significantly lower in CRC tissues than in normal tissues (Li et al., 2012). In our research, these three miRNAs were up-regulated in HCT-116 cells treated with QUE. These observations suggested that miR-320b, miR-320c, and miR-320d may play vital roles in anti-CRC mechanism of QUE. MiR-125b-2-3p, significantly down-regulated miRNA in CRC, was also pointed out as a novel diagnostic and prognostic biomarker in human CRC (Zhou et al., 2018). Our data showed that miR-125b-2-3p was up-regulated in QUEtreated HCT-116 cells. These studies suggested that miR-125b-2-3p may play an important role in the anti-CRC mechanism of QUE. Furthermore, some DEmiRNAs, such as novel_miR_873, novel_miR_710, novel_miR_20, and novel_miR_885, were also found to be significantly different between CRC cells with and without treatment of QUE, which suggested that these miRNAs play an important role in the anti-CRC mechanism of QUE. The relationships between these genes and CRC were firstly reported.
We noticed that a significant GO term of DElncRNAs and their target gene were related with DNA repair. This phenomenon Increasing evidence indicate that circRNAs can influence miRNA activity as endogenous sponges and affect mRNA splicing and transcription by interacting with the Pol II complex in the nucleus (Tay et al., 2014;Xie et al., 2017). As many circRNAs failed to be allocated to functional modules, little public data about these circRNAs could be found. In this study, based on the constructed circRNA-miRNA-mRNA co-expression network, we observed that many circRNAs contained one or more miRNA binding sites. Thus, circRNA (8:93786223|93822563) was able to interact with LRG1, through competitively binding with miR-5096. This competitively binding mode was similar with ENST00000313807 and ENST00000449307. Therefore, further study was deserved to reveal the interaction relationships of circRNA (8:93786223|93822563)-miR-5096-LRG1 in QUE's action mechanism. (I-L) Unsupervised clustering analysis showing expression profiles of top 40 differentially expressed miRNAs, lncRNAs, circRNAs, and mRNAs between HCT116 cells and QUE-treated HCT116 cells. QUE-1, QUE-2, and QUE-3 represent three groups of HCT-116 cells treated with QUE. CON-1, CON-2, and CON-3 represent three groups of untreated HCT-116 cells. Abbreviation: QUE, quercetin. CON, control.
Although altered ncRNAs and mRNAs were identified and their possible roles in anti-CRC mechanisms of QUE were investigated, several limitations should be considered in interpreting our findings. Firstly, a previous study reported that maximum plasma concentrations of QUE after the ingestion of 100 mg were lower than 10 µM (Graefe et al., 2001), and only one-tenth of the lower QUE concentration is used in this in vitro assays. Recently, the daily dose of QUE in a clinical study was up to 1250 mg for three consecutive days/week (Justice et al., 2019). Merging these results, one may speculate that those effective tissue levels could be attained. Thus, we selected relatively high concentrations of QUE to investigate its effects on the CRC cell line. In future studies, we could perform an analysis by exposing cells to lower QUE concentrations, aiming to investigate whether lower concentrations of QUE have anti-CRC efficacy and different impacts on gene expression profiles. Secondly, the analysis was only performed on HCT-116 cancer cells. Global ncRNA and mRNA changes in other CRC cell lines and CRC model animals treated with QUE should be also determined in further studies to more accurately reflect the anti-CRC mechanisms of QUE. Thirdly, RNA-sequence technology should be applied to unravel previously inaccessible transcriptome complexities, and because the functions of ncRNAs remain largely unknown, the comprehension of our data was not straightforward (Wang et al., 2018). Finally, further research should select more ncRNAs and mRNAs for sequencing data validation and observe whether the apoptotic process of HCT-116 cells will change by manipulating the expression of top nRNA and mRNA candidates. To solve this problem, further studies are now highly warranted.

CONCLUSION
In summary, changes in the coding and non-coding transcriptomes of HCT-116 cells without and with QUE intervention were identified, and combining it with bioinformatics analysis may provide a better understanding of the potential roles of miRNA, lncRNA, circRNA, and mRNA in CRC and QUE treatment. To better reveal the mechanism of QUE-mediated apoptosis in HCT-116 cells, we now need to conduct further experiments through manipulating the action of the top potential candidates in this study. Moreover, the corresponding roles and molecular mechanisms of these ncRNA and mRNA need to be further explored.

DATA AVAILABILITY
All the relevant data are contained within the manuscript.

AUTHOR CONTRIBUTIONS
BY contributed to the conception and the design of this study. ZZ, BL, and PX conducted the experiments. ZZ and BY drafted the main manuscript text. All the authors participated in the interpretation of results. All the authors have read and approved the final manuscript.

FUNDING
This work was financially supported by the Hunan Science and Technology Department Project (No. 2014SK3039) and the implementation plan issued by the Office of the State Administration of Traditional Chinese Medicine for the fifth batch of national academic succession work for veteran Chinese medicine experts (2012, No. 40).
FIGURE 7 | The KEGG enrichment analysis of the cis targeted genes of DElncRNA (A), the trans targeted genes of DElncRNA (B), DEmRNA (C), hostgenes of DEcirRNA (D), and targeted genes of DEmiRNA (E). The abscissa GeneRatio represents the proportion of genes of interest in the pathway, and the ordinate represents each pathway. The size of the dots represents the number of genes annotated in the pathway, and the color of the dots represents the corrected p value of the hypergeometric test. Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.