Analysis of mRNAs and ncRNAs and Prediction Competing Endogenous RNA Networks in Colorectal Cancer Chemo-Resistance by Whole-Transcriptome Sequencing

Background: Chemo-resistance is a major clinical obstacle to the treatment of colorectal cancer (CRC), mRNAs and non-coding RNAs (ncRNAs) have been reported to modulate the development of chemo-resistance. However, the proles of mRNAs and ncRNAs as well as competing endogenous RNA (ceRNA) networks in CRC chemo-resistance are still unclear, and whether different drug resistance of CRC have the same mechanisms also needs to be explored. This study aims to uncover the expression of mRNAs and ncRNAs in parental cell lines and different chemo-resistant cell lines, and construct ceRNA regulatory networks by whole-transcriptome sequencing. Methods: The expression of mRNAs and ncRNAs in parental cell lines and drug-resistant cell lines were identied by whole-transcriptome sequencing and bioinformatics methods. Results: A total of 1779 mRNAs, 64 miRNAs, 11 circRNAs and 295 lncRNAs were common differentially expressed in two different chemo-resistant cell lines when compared with the control. In addition, 5,767 lncRNA-miRNA-mRNA relationship pairs and 47 circRNA-miRNA-mRNA pathways were constructed according to ceRNA regulatory rules, in which AC109322.2-hsa-miR-371a-5p-BTNL3 and hsacirc_027876-hsa-miR-582-3p-FREM1 were identied as the most potential ceRNA networks involved in drug resistance to CRC. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of two ceRNA regulatory networks showed that the TNF signaling pathway may be crucial in the process of CRC drug resistance. Conclusions: A large number of mRNAs and ncRNAs in chemo-resistant cell lines were different expressed, which may play pivotal roles in development of drug resistance through the ceRNA regulatory network. This study may improve our understanding of the underlying mechanisms and provide a promising therapeutic strategy for CRC chemo-resistance. and sensitivity to tamoxifen. Molecular mechanisms study further elucidated that hsa_circ_0025202 upregulated FOXO3 expression by acting as a ceRNA to affect the expression of miR-182-5p, indicating hsa_circ_0025202 could be a promising biomarker for tamoxifen resistance [28]. At present, the ceRNA regulatory networks in CRC chemo-resistant are still unclear and needs to be further explored. In this study, we carried out whole transcriptome sequencing technology (RNA-sequencing) to analysis proles of mRNAs and ncRNAs in chemo-resistant CRC cell lines (HCT8/5-Fu and HCT8/DDP) and chemo-sensitive cells (HCT8). Differentially expressed (DE) mRNAs and ncRNAs were identied. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were subsequently performed to understand the function of differentially expressed genes. The ceRNA networks were contructed and predicted among the common differently expressed mRNAs and ncRNAs in different CRC cell lines. This study may provide new inspirations of the mechanisms in chemo-resistance and a theoretical basis for the clinical treatment of CRC. to analyze the proles of mRNAs, miRNAs, lncRNAs and circRNAs between two chemo-resistant CRC cell lines and the parental cell lines through high-throughput sequencing technology. CeRNA regulatory networks were predicted, including lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA, GO and KEGG pathway analyses were used to explore potential functions of differentially expressed genes, of them, TNF signaling pathway was both enriched in two ceRNA networks, demonstrating this pathway may be related to CRC 5-Fu and DDP resistance. Our study might provide novel insights into the progress of chemoresistance and facilitate to improve the strategy for the treatment of advanced colorectal cancer patients.


Background
Colorectal cancer (CRC), a common malignant tumor of the digestive system, ranks third in terms of incidence and the second leading cause of cancer-related mortality around the world [1]. At present, surgery is still an effective treatment for CRC, but for advanced-stage CRC patients who are not surgical candidates, chemotherapy has become an optimal method [2]. Cytotoxic drugs commonly used in clinical treatment of CRC include 5-uorouracil (5-Fu) and cisplatin (DDP) [3]. However, the therapeutic effect is limited and the 5-year survival rate of patients is remains low due to the development of drug resistance [4,5]. The underline mechanisms of chemo-resistance for CRC need to be fully identi ed to improve the e cacy of chemotherapy and prognosis of patients.
Non-coding RNAs (ncRNAs) are RNA molecules that lack of coding protein ability, but play a pivotal role in biological processes of many cancers, including cell proliferation, migration, metastasis, apoptosis and chemo-resistance [6][7][8]. In fact, the genomes of eukaryotes are pervasively transcribed, while less than 2% of the transcripts are protein-coding mRNAs, and most of the rest are ncRNAs [9]. Few years ago, ncRNAs, such as microRNAs (miRNAs), siRNAs and PIWIinteracting RNAs (piRNAs), have been researched a lot [10,11]. It was found that miRNAs, around 22 nucleotides in length, exert an important impact on the progression of cancers via binding to the 3' untranslated regions (UTRs) of their target mRNAs to regulate gene expression [12,13]. Recently, with the development of high-throughput sequencing and bioinformatics technologies, a variety of ncRNAs, such as long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs), have been explored and revealed [14,15]. Accumulating evidence shows that lncRNAs and circRNAs are usually differentially expressed in various cancers and function as oncogenes or tumor suppressor genes [16][17][18]. It was also reported that some ncRNAs are related to the regulation of chemoresistance in CRC. For example, lncRNA-XIST was dramatically up-regulated in CRC tissues and cell lines, the chemo-resistance of CRC cells to 5-uorouracil, cisplatin, mitomycin and adriamycin was enhanced when overexpression XIST. Mechanistically, XIST may alter chemosensitivity of CRC cells at least partly through regulation of miR-30a-5p/ROR1 axis [19]. Additionally, Kha Wai Hon et al indicated that exosomal has_circ_0000338 was highly expressed in FOLFOXresistant HCT116 cells compared with parental cells, and knockdown circ_0000338 could improve the chemo-sensitivity, demonstrating circ_0000338 may exhibit regulatory roles in chemo-resistant of CRC [20]. Besides, it was found that miR-148a is down-regulated in cisplatin-resistant CRC cells, and overexpression of miR-148a could inhibit wnt10b expression and the activation of β-catenin signaling pathway, thereby increasing cell invasion and migration ability as well as chemo-sensitivity [21].
Mounting evidence has demonstrated lncRNAs and circRNAs can serve as sponges of miRNAs via competitive endogenous RNA (ceRNAs) network to eliminate the inhibitory effect of miRNAs on their target genes [22][23][24]. For lncRNA, circRNA and mRNA, they may compete with each other for the binding miRNAs to form a ceRNA regulatory network, mainly includes lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks [25,26]. These ceRNA networks might contribute to progression of several cancers, including chemo-resistance. For instance, Zhang et al indicated that lncRNA KCNQ1OT1 is upregulated in chemo-resistant tongue carcinoma tissues compared with chemo-sensitive specimens. KCNQ1OT1 competitively bound with miR-211-5p to increase the expression of target gene Ezrin, resulting in the activation of Ezrin/Fak/Src signaling pathway and enhanced chemo-resistance to cisplatin [27]. In addition, a report by Sang et al found that lowly expressed hsa_circ_0025202 in tamoxifen-resistant breast cancer cell line could inhibit cell proliferation and migration, while increase cell apoptosis and sensitivity to tamoxifen. Molecular mechanisms study further elucidated that hsa_circ_0025202 upregulated FOXO3 expression by acting as a ceRNA to affect the expression of miR-182-5p, indicating hsa_circ_0025202 could be a promising biomarker for tamoxifen resistance [28].
At present, the ceRNA regulatory networks in CRC chemo-resistant are still unclear and needs to be further explored. In this study, we carried out whole transcriptome sequencing technology (RNA-sequencing) to analysis pro les of mRNAs and ncRNAs in chemo-resistant CRC cell lines (HCT8/5-Fu and HCT8/DDP) and chemo-sensitive cells (HCT8). Differentially expressed (DE) mRNAs and ncRNAs were identi ed. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were subsequently performed to understand the function of differentially expressed genes. The ceRNA networks were contructed and predicted among the common differently expressed mRNAs and ncRNAs in different CRC cell lines. This study may provide new inspirations of the mechanisms in chemo-resistance and a theoretical basis for the clinical treatment of CRC.

Cell culture
The human CRC cell lines HCT8 were purchased from China Center for Type Culture Collection (Wuhan, China) and cultured in RPMI-1640 (Gibco, USA) supplemented with 10% fetal bovine serum (Gibco, USA). The chemo-resistant CRC cell lines (HCT8/5-Fu and HCT8/DDP) were separately established via stepwise exposure to 5-uorouracil (5-Fu) and cisplatin (DDP) (Sigma-Aldrich, USA) for about seven months, and cultured in RPMI-1640 medium plus with 10% fetal bovine serum as well as 5 µg/ml 5-Fu or 1 µg/ml DDP to maintain drug resistance at 37 ℃ in a 5% CO 2 incubator.
MTT assay MTT assay was used to detect cell viability. The chemo-sensitive and resistant CRC cell lines at the logarithmic growth phase were seeded into 96-well plates at a density of 5×10 3 cells/well for overnight culture. Then different concentrations (0, 5, 10, 20, 40, 80 μg/ml) of 5-Fu or (0, 3.75, 7.5, 15, 30, 60 μg/ml) DDP were added to the HCT8 and HCT8/5-Fu or HCT8/DDP cells, with ve replicate wells for each concentration. After 72 h treatment, 100 μL 5mg/ ml MTT (sigma, USA) was added into the wells and cultured at 37 °C for 4 h. The medium was removed and 100μL DMSO was added into the wells for 15 minutes, and the absorbance of each well was measured at 490 nm by Microplate reader, the viability was calculated and IC 50 values were assessed.

RNA extraction and quali cation
Total RNA was extracted from the cell samples using TRIzol reagent (Invitrogen, USA) according to the manufacturer's protocol. The concentration and purity of the total RNA were detected by Nanodrop2000 spectrophotomete (Thermo Scienti c, USA). The RNA integrity was measured by agarose gel electrophoresis and Agilent 2100 Bioanalyzer (Agilent, USA). The RNA Integrity Number (RIN) was generated by the 2100 BioAnalyzer to re ect the quantitative value of RNA integrity. Total RNA ≥1ug, concentration≥100 ng/µL, RIN≥8.0, OD260/230≥1.0 and the value of OD260/280 should between 1.8 and 2.2 were required before library construction.

Construction of strand-speci c library
Three pairs of chemo-resistant CRC cell lines (HCT8/5-Fu and HCT8/DDP) and parental cell lines HCT8 were selected for deep Illumina sequencing to identify the differentially expressed mRNAs and ncRNAs. The rRNA in the total RNA was removed by the Ribo-Zero Magnetic kit (EpiCentre, USA) according to the manufacturer's instructions, and the remaining RNAs were interrupted to fragments about 300bp in length. Random hexamer primers and reverse transcriptases were used to synthesize rst-strand complementary DNA (cDNA), dUTPs were used instead of dTTPs when synthesizing second-strand cDNA. 3 'end of the double-stranded cDNA were added an A base connected the Y-shaped linker, and the second strand of cDNA was digested with UNG enzyme, then polymerase chain reaction (PCR) was used to enrich fragments and QuantiFluor® dsDNASystem (Promega, USA) was used to quantify libraries. Paired-end (PE) sequencing (150 bp) of each library was constructed by Illumina HiSeq platform.

Construction of small RNA Library
Total RNA was precipitated with polyethylene glycol and was ligated by 3' adapters followed by hybridization of multiplex SR RT primers, and ligation of multiplex 5' SR adapters at 5' ends of miRNAs. First-strand cDNA was synthesized using Truseq TM Small RNA sample prep Kit (Illumina, USA), after reverse transcription reaction, PCR ampli cation was performed for 11-12 cycles, and the products were puri ed by 6% Novex TBE PAGE gel. QuantiFluor® dsDNASystem (Promega, USA) was used to quantify libraries, and Single End (SE) sequencing (50 bp) was performed to generate at least 15 million reads per sample on an Illumina Hiseq X-ten sequencer (Illumina, USA).

Reads mapping and transcript assembly
The samples were sequenced based on Illumina platform and then generated raw data, cutadapt was used to remove reads containing adaptors and lowquality reads to obtain clean data [29]. Content and quality distribution of bases, as well as average quality distribution of reads were conducted to assess data quality. HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml) was used to map clean data to the reference genome, then reads aligned to the genome were performed comparison regional distribution and gene coverage uniformity analysis [30]. Finally, RNAs expression level was estimated by calculating fragments per kilobase per million reads (Fig. S1).

RNA expression and Functional enrichment analysis
Differentially expressed RNAs were screened with |log 2 fold change (FC) | > 1 and p < 0.05, and then were performed function annotation through Gene Ontology (GO) analysis (http://geneontology.org/). GO enrichment analyses of differentially expressed genes were classi ed into three categories, including biological process (BP), molecular function (MF) and cell component (CC), and the top 10 signi cantly enriched GO terms in each category were displayed. Pathway analysis can further understand the genes involved in metabolic pathways and speci c biological functions. Kyoto Encyclopedia of Genes and Genomes (KEGG) is the main public database (http://www.kegg.jp/) which can determine the biological pathway of the differentially expressed genes signi cantly enriched, and p < 0.05 was considered as signi cant.
Competing endogenous RNA (ceRNA) network analysis As shown in Fig. S2, LncRNA (cirRNA)-miRNA-mRNA regulatory network were constructed based on differently expressed mRNAs and ncRNAs according to the ceRNA hypothesis. First, miRNA-target genes, including mRNA, lncRNA, cirRNA, were predicted by using miRanda and TargetScan software, and the correlation between lncRNA (cirRNA) and mRNA expression was calculated via pearson's correlation coe cient, mRNA and lncRNA (cirRNA) in the ceRNA network generally have a positive correlation. Pearson's correlation coe cient >0.7 and p <0.05 were selected for further analysis. Then, for lncRNA (cirRNA) and mRNA that have a positive correlation and contain common miRNA response elements (MERs), the number of miRNAs they shared was conducted through hypergeometric test. According to the expression correlation between lncRNA (cirRNA)-miRNA-mRNA, two scoring methods were used to evaluate the regulatory function of miRNA on a pair of competitive lncRNA (cirRNA) -mRNA. 1) Regulate similarity score to compare the similarity between miRNAs-lncRNA (cirRNA) expression correlation and miRNAs-mRNA expression correlation; 2) Sensitivity correlation score, and the average value of lncRNA (cirRNA)-miRNA-mRNA correlation were used as the sensitivity correlation between lncRNA (cirRNA) and mRNA. Finally, lncRNA (cirRNA)-mRNA ceRNA networks (sensitivity correlation score > 0.3) were visualized using Cytoscape v3.8 software.

Statistical analysis
The differences of mRNA and non-coding RNAs with |log 2 fold change | > 1 and p < 0.05 between chemo-sensitive and chemo-resistant cells were considered to be statistically based on RNA-sequencing analysis. All experiments were repeated at least three times, and statistical analysis was performed by independent-sample t test using SPSS statistical software version 17.0. All data were expressed as the means ± standard deviation (SD). p <0.05 was regarded statistically signi cant.

Establishment and identi cation of colorectal cancer drug-resistant cell lines
To generate the chemo-resistant CRC cell lines (HCT8/5-Fu and HCT8/DDP), the parental cells HCT8 were exposed to increasing concentrations of 5uorouracil (5-Fu) or cisplatin (DDP) for more than seven months. MTT assay showed that the cell viability was decreased when exposed to the increasing

Screening of DE-miRNAs and its functional enrichment analysis in CRC drug-resistant cell lines
From whole-transcriptome sequencing data, a total of 1542 miRNAs was found in HCT8/5-Fu cells, including 98 (52 upregulation and 46 downregulation) miRNAs with signi cant difference (Fig. 3A). And HCT8/DDP cells were found to have 1361 miRNAs, with 79 DE-miRNAs (50 upregulation and 29 downregulation) (Fig. 3B). A distinguishable miRNAs expression pro le among samples was displayed by hierarchical clustering analysis (Fig. 3C). There are 64 common DE-miRNAs between two comparison groups (Fig. 3D), and the top 10 up-and down-regulated common miRNAs in two chemo-resistant cells were listed in Table S2, hsa-miR-205-5p and hsa-miR-452-5p was the most upregulated and downregulated miRNAs, respectively. Besides, the enriched GO functions for the target genes of miRNAs in HCT8/5-Fu and HCT8/DDP cells were shown in Fig. 3E and 3F, the results showed that the majority of the 10 most enriched CC, MF, BP terms among two chemo-resistant cells are the same, and the BP enrichment of DE-miRNAs mainly related to anatomical structure development and anatomical structure morphogenesis. Pathway analysis results indicated that DE-miRNAs most enriched in inositol phosphate metabolism, endocrine and other factor−regulated calcium reabsorption and TGF-β signaling pathway in both of drug-resistant cells (Fig. 3G and 3H).

Analysis of differentially expressed circRNAs and functional enrichment in CRC drug-resistant cells
The RNA-sequencing results indicated that a total of 7393 circRNAs were screened out in HCT8/5-Fu cells, among which 48 circRNAs were differently expressed, with 16 up-regulated and 32 down-regulated (Fig. 4A). In addition, 90 DE-circRNAs (42 upregulation and 48 downregulation) were found among 7385 circRNAs in HCT8/DDP cells (Fig. 4B). The heat map analysis showed the expression of the DE-circRNAs visually, the three repeats of each group clustered together, while the chemo-resistant groups and control group were clustered separately (Fig. 4C). Besides, only 11 circRNAs were consistently expressed differently in both chemo-resistant groups (Fig. 4D), which were listed in Table S3, hsacirc_023607 and hsacirc_027876 was the most up-regulated and down-regulated DE-circRNAs, respectively. Additionally, GO enrichment (Fig. 4E, 4F) and KEGG pathway analysis (Fig. 4G, 4H)  up-regulated and 342 down-regulated (Fig. 5A). Meanwhile, and 601 (274 upregulation and 327 downregulation) out of all detected 22,001 lncRNAs were differentially expressed in cisplatin resistant cells (Fig. 5B). As shown in Fig. 5C, a heat map showed the general expression pro les of DE-lncRNAs were clearly different in control and drug-resistant groups. 295 DE-lncRNAs were common in two comparisons (Fig. 5D), the top 10 up-and down-regulated common lncRNAs in two drug-resistant groups were listed in Table S4, the most upregulated lncRNAs was MSTRG.30161.48 and the most downregulated lncRNAs was MSTRG.22339.1. The cis-target genes 10 kb upstream and downstream and trans-target genes with pearson's correlation > 0.95 and p value <0.05 of the DE-lncRNAs in HCT8/5-Fu cells were selected for GO and KEGG pathway analyses. The GO enrichment results of cis-target genes of DE-lncRNAs were shown in Fig. 5E and 5F, and KEGG pathway analysis indicated that Hippo signaling pathway was the most signi cant enriched pathway in two drugresistant groups (Fig. 5G and 5H). In addition, the GO functions and KEGG pathways for trans-target genes of DE-lncRNAs in HCT8/5-Fu and HCT8/DDP cells were presented in Fig. S3.

Construction of ceRNAs network of DE-lncRNAs and DE-circRNAs
It is well documented that competing endogenous RNAs (ceRNAs) members can regulate each other through competing for the same miRNA response elements (MREs). To further understand the roles and the molecular mechanisms of the dysregulated mRNAs and ncRNAs in drug-resistant CRC cell lines, lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA networks were constructed based on RNA-sequencing data. Co-expressed genes from differentially expressed RNAs in two chemo-resistant cell lines were screened to construct ceRNA networks according to the following screening conditions: 1) Pearson's correlation coe cient between DE-lncRNAs/DE-circRNAs and DE-mRNAs was > 0.7, and sensitivity correlation was > 0.3. 2) Compared with parental cells, the expression trends of lncRNAs/circRNAs and mRNAs are the same in two drug-resistant cell lines, while the variation trend of miRNA is opposite to that of lncRNAs/circRNAs and mRNAs. At last, there were 5,767 lncRNA-miRNA-mRNA and 47 circRNA-miRNA-mRNA relationship pairs.
In the lncRNA-related ceRNA network, lncRNAs was upregulated in 626 relationship pairs of 5,767 pathways (Fig. 6A). The fold change of lncRNAs upregulation in two chemo-resistant cells were in nity in 307 pathways, in which miRNAs were in nitely downregulated in 29 pathways; Besides, lncRNAs and mRNAs were in nitely upregulated in 6 pathways (Fig. 6B), indicating these networks may act as core regulators of chemo-sensitivity, and the relationship pairs were shown in Table S5. Besides, AC109322.2-hsa-miR-371a-5p-BTNL3 network was identi ed to be involved in drug resistance of CRC, with the fold change of lncRNAs and mRNAs upregulation were both in nity, and the hsa-miR-371a-5p was downregulated with 9.6-fold change that was largest fold changes in 6 pathways, and the heatmap of this ceRNA pathway was shown in Fig. S4. In addition, there were 5,141 relationship pairs with lncRNAs downregulation (Fig. 6C), it was found that lncRNAs and mRNAs were in nitely downregulated, while miRNAs were in nitely upregulated in 29 pathways ( Fig.  6D and Table S6).
In the circRNA-associated ceRNA network, a total of 47 lncRNA-miRNA-mRNA pathways were constructed including 5 circRNAs, 9 miRNAs, and 33 mRNAs (Fig. 6E). Interestingly, the sequencing results did not nd that circRNAs with highly expression meet the requirements of ceRNA network construction. However, it was found that circRNAs were in nitely downregulated in 13 pathways, in which mRNAs lowly expressed with in nite-fold change in 4 pathways ( Fig. 6F), the details of these ceRNA networks were shown in Table S7. In addition hsacirc_027876-hsa-miR-582-3p-FREM1 was believed to be related to the development of CRC chemo-resistance, with hsa-miR-582-3p was upregulated with 5.6-fold change, and the heatmap of this ceRNA network was shown in and biological process (BP), were endomembrane system, 3',5'−cyclic−AMP phosphodiesterase activity and organic hydroxy compound metabolic process, respectively. KEGG pathway analysis showed that the most enriched pathways in these DE-mRNAs involved in the circRNA-miRNA-mRNA networks were terpenoid backbone biosynthesis, TNF signaling pathway, arrhythmogenic right ventricular cardiomyopathy (ARVC) and AMPK signaling pathway (Fig. 7B). KEGG analysis of the two regulatory networks showed that the TNF signaling pathway may be pivotal in the process of CRC drug resistance.

Discussion
CRC is a high-risk digestive system malignant tumor that causes a large number of deaths every year, seriously in uencing health and lives of patients. Chemotherapy is the main treatment approach for patients with advanced or inoperable colorectal cancer [31]. However, chemo-resistance is a huge obstacle for CRC treatment because patients usually possess intrinsic or acquired drug resistance thus leading to poor clinical treatment [32].
At present, a large number of studies have shown that, in addition to mRNAs, non-coding RNAs also related to drug resistance of cancers, including gastric cancer [33], bladder cancer [34], lung adenocarcinoma [35]. Moreover, in CRC drug resistance studies, it was reported that LINC00957 could increase 5uorouracil resistance and result in poor survival of patients [36], hsa_circ_0079662 was related to oxaliplatin resistance via TNF-α pathway [25], microRNA-375-3p was involved in 5-uorouracil resistance by targeting thymidylate synthase [37].
In addition, lncRNAs and circRNAs have been reported to function as miRNAs sponges to regulate target mRNAs expression that related to biological processes in many cancers [38][39][40]. Microarray analysis and high-throughput sequencing have been used to analyze the mRNAs or ncRNAs pro ling during the development of drug resistance in many types of cancers. Zhu et al [41] analyzed the comprehensive expression pro le and ceRNA regulatory networks between mRNAs and ncRNAs in the osteosarcoma chemo-resistance through RNA-sequencing technology. A report by Li et al [42] identi ed deregulated mRNAs and lncRNAs in drug-resistant pancreatic cancer by next-generation RNA sequencing and constructed the co-expression network of differently expressed lncRNAs and mRNAs. Besides, Wang et al [43] explored the expression of chemoresistance-associated miRNAs in breast cancer by microarray analysis. However, the study about pro les of mRNAs and ncRNAs and involved ceRNA networks is seldom in chemo-resistance of CRC and whether different drug resistance of CRC have the same mechanisms also needs to be explored. In order to explore the underline mechanisms of chemo-resistance in CRC, in the present study, we identi ed the pro les of mRNA and ncRNAs between chemo-resistant (HCT8/5-Fu and HCT8/DDP) and chemo-sensitive cells (HCT8) through whole transcriptome sequencing, further compared the expression and function of dysregulated genes, and constructed ceRNA regulatory networks that may be related to different drug resistance in CRC. Our results indicated that compared with the parental cells, there were 1779 mRNAs, 64 miRNAs, 11 circRNAs and 295 lncRNAs were common differentially expressed in two chemo-resistant cell lines. GO enrichment of mRNAs and ncRNAs indicated that there are many same terms in the 10 most enriched CC, MF, BP terms among two chemo-resistant cells. The GO terms at the BP level, such as regulation of cell migration and motility, DNA repair, regulation of signaling, were involved in chemo-resistance [44][45][46]. The results of KEGG pathway showed that dysregulated mRNAs and ncRNAs mainly enriched in MAPK, TGF-β, and Hippo signaling pathway, it was proved that these pathways were taken part in the chemoresistance of CRC [47][48][49].
Non-coding RNAs have been regarded as "dark matter" or "by-product" in the past, subsequently Salmena et al. rst proposed the "ceRNA hypothesis" which suggested that lncRNA can affect mRNAs expression by competitively combining to miRNAs via miRNA response elements [50,51]. The lncRNAs (circRNAs)associated ceRNA mechanisms have been researched in the biological processes of CRC [52,53]. In the present study, the ceRNA regulatory networks in the CRC chemo-resistance were comprehensively constructed based on RNA-sequencing data. A total of 5,767 lncRNA-miRNA-mRNA pathways and 47 circRNA-miRNA-mRNA pathways were predicted. It was further ltered and found that fold change of lncRNAs and mRNAs upregulation were in nity simultaneously in 6 pathways in lncRNA-related ceRNA network, in which AC109322.2-hsa-miR-371a-5p-BTNL3 was considered as the most potential ceRNA networks involved in CRC drug resistance due to its huge fold change in drug resistant cells. In addition, circRNAs and mRNAs were in nitely downregulated with in nite-fold change in 4 pathways in circRNA-associated ceRNA network, these ceRNA networks may play a more signi cant role in the development of chemo-resistance, especially for hsacirc_027876-hsa-miR-582-3p-FREM1 pathway. Subsequently, KEGG pathway analysis of ceRNA networks were conducted, the result indicated TNF signaling pathway was signi cantly enriched in two networks, it might indicated some ceRNA networks in uence CRC drug resistance by regulating TNF signaling pathway.
There are limitations in our study, the validation of the top 10 upregulated mRNAs and ncRNAs in drug resistant CRC cell lines were not performed, and the mechanism of the ceRNA networks were not con rmed, further research still needs to be performed in the future. To our knowledge, this is the rst report to analyze the pro les of mRNAs, miRNAs, lncRNAs and circRNAs between two chemo-resistant CRC cell lines and the parental cell lines through high-throughput sequencing technology. CeRNA regulatory networks were predicted, including lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA, GO and KEGG pathway analyses were used to explore potential functions of differentially expressed genes, of them, TNF signaling pathway was both enriched in two ceRNA networks, demonstrating this pathway may be related to CRC 5-Fu and DDP resistance. Our study might provide novel insights into the progress of chemoresistance and facilitate to improve the strategy for the treatment of advanced colorectal cancer patients.

Conclusion
In conclusion, our ndings showed a comprehensive expression pro le of mRNAs and ncRNAs. We also constructed a lncRNA-miRNA-mRNA and circRNA-miRNA-mRNA ceRNA regulatory networks for drug resistance in CRC. In particular, we found that AC109322.2-hsa-miR-371a-5p-BTNL3 and hsacirc_027876hsa-miR-582-3p-FREM1 were the most potential ceRNA networks involved in drug resistance. KEGG pathway analysis of two ceRNA regulatory networks showed that the TNF signaling pathway may be crucial in the process of CRC drug resistance. Although more mechanisms need to be further investigated, but this study also provide a promising therapeutic strategy for CRC chemo-resistance.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.