*2.1. Identification of Differential Expressed (DE) and Novel Non-Coding RNAs (ncRNAs)*

In our results, 9848 lncRNAs were found: 84 were known lncRNAs, and 9764 were novel lncRNAs found in the control and chilling samples (Table S1). Among them, most of the lncRNAs were lincRNAs (8022, 81.5%), followed by antisense-lncRNAs (919, 9.3%), sense lncRNAs (682, 6.9%), and intronic-lncRNAs (225, 2.3%) (Figure 1A). In addition, 213 novel circRNAs were found, with many emanating from chromosome 8 (Table S1). The majority of circRNAs were over 3000 nt and from intergenic regions, while additional circRNAs were between 400 to 800 nt and derived from exons (Figure 1B, Table S1). In total, 281 miRNAs were found in our libraries with 120 known and 161 novel miRNAs. Most of the novel miRNAs were between 21 and 24 nt. The miRNAs nucleotide bias was

also analyzed in our results, and, intriguingly, we found that the first nucleic acid bases were U and A, while the last was G (Table S1). bias was also analyzed in our results, and, intriguingly, we found that the first nucleic acid bases were U and A, while the last was G (Table S1). bias was also analyzed in our results, and, intriguingly, we found that the first nucleic acid bases were U and A, while the last was G (Table S1).

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 3 of 15

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 3 of 15

**Figure 1.** The four kinds of long non-coding RNAs (lncRNAs) were: long intergenic noncoding RNAs (lincRNAs) (8022, 81.5%), antisense-lncRNAs (919, 9.3%), sense lncRNAs (682, 6.9%), and introniclncRNAs (225, 2.3%) (**A**). **Figure 1.** The four kinds of long non-coding RNAs (lncRNAs) were: long intergenic noncoding RNAs (lincRNAs) (8022, 81.5%), antisense-lncRNAs (919, 9.3%), sense lncRNAs (682, 6.9%), and intronic-lncRNAs (225, 2.3%) (**A**). **Figure 1.** The four kinds of long non-coding RNAs (lncRNAs) were: long intergenic noncoding RNAs (lincRNAs) (8022, 81.5%), antisense-lncRNAs (919, 9.3%), sense lncRNAs (682, 6.9%), and introniclncRNAs (225, 2.3%) (**A**).

We compared the expression profiles of lncRNAs, circRNAs, miRNAs, and mRNAs between the control and chilling groups, and found that 380 lncRNAs, 36 circRNAs, 18 miRNAs, and 4128 mRNAs were differentially expressed (Figure 2, Table S2). Among them, 198 lncRNAs, 28 circRNAs, 3 miRNAs, and 2833 mRNAs were upregulated, whereas 182 lncRNAs, 8 circRNAs, 15 miRNAs, and 1295 mRNAs were downregulated in the chilling sample compared with the control. The differentially expressed non-coding RNAs are listed in Supplementary Table S2. The differentially expressed lncRNAs and mRNAs were widely distributed on the autosomal chromosomes, while the differentially expressed circRNAs were not found in chromosomes 5, 9, and 11, and their number was the largest in chromosome 1. We compared the expression profiles of lncRNAs, circRNAs, miRNAs, and mRNAs between the control and chilling groups, and found that 380 lncRNAs, 36 circRNAs, 18 miRNAs, and 4128 mRNAs were differentially expressed (Figure 2, Table S2). Among them, 198 lncRNAs, 28 circRNAs, 3 miRNAs, and 2833 mRNAs were upregulated, whereas 182 lncRNAs, 8 circRNAs, 15 miRNAs, and 1295 mRNAs were downregulated in the chilling sample compared with the control. The differentially expressed non-coding RNAs are listed in Supplementary Table S2. The differentially expressed lncRNAs and mRNAs were widely distributed on the autosomal chromosomes, while the differentially expressed circRNAs were not found in chromosomes 5, 9, and 11, and their number was the largest in chromosome 1. We compared the expression profiles of lncRNAs, circRNAs, miRNAs, and mRNAs between the control and chilling groups, and found that 380 lncRNAs, 36 circRNAs, 18 miRNAs, and 4128 mRNAs were differentially expressed (Figure 2, Table S2). Among them, 198 lncRNAs, 28 circRNAs, 3 miRNAs, and 2833 mRNAs were upregulated, whereas 182 lncRNAs, 8 circRNAs, 15 miRNAs, and 1295 mRNAs were downregulated in the chilling sample compared with the control. The differentially expressed non-coding RNAs are listed in Supplementary Table S2. The differentially expressed lncRNAs and mRNAs were widely distributed on the autosomal chromosomes, while the differentially expressed circRNAs were not found in chromosomes 5, 9, and 11, and their number was the largest in chromosome 1.

**Figure 2.** Differentially expressed (DE) ncRNAs: 380 lncRNAs, 36 circRNAs, 18 miRNAs, and 4128 mRNAs were found differentially expressed between the control and chilling groups. **Figure 2.** Differentially expressed (DE) ncRNAs: 380 lncRNAs, 36 circRNAs, 18 miRNAs, and 4128 mRNAs were found differentially expressed between the control and chilling groups. **Figure 2.** Differentially expressed (DE) ncRNAs: 380 lncRNAs, 36 circRNAs, 18 miRNAs, and 4128 mRNAs were found differentially expressed between the control and chilling groups.

#### *2.2. GO and KEGG Pathway Analyses of ncRNAs 2.2. GO and KEGG Pathway Analyses of ncRNAs 2.2. GO and KEGG Pathway Analyses of ncRNAs*

To explore the potential functions of the differential expressed non-coding RNAs, we performed GO and KEGG analyses. To our knowledge, the lncRNAs could regulate the expression of neighboring and overlapping coding genes; hence, lncRNAs likely regulate related mRNA genes [30]. To explore the potential functions of the differential expressed non-coding RNAs, we performed GO and KEGG analyses. To our knowledge, the lncRNAs could regulate the expression of neighboring and overlapping coding genes; hence, lncRNAs likely regulate related mRNA genes [30]. To explore the potential functions of the differential expressed non-coding RNAs, we performed GO and KEGG analyses. To our knowledge, the lncRNAs could regulate the expression of neighboring

The cis and trans targets of the differentially expressed lncRNAs were both analyzed, and the most relevant GO terms associated with biological processes and molecular functions contained many

The cis and trans targets of the differentially expressed lncRNAs were both analyzed, and the most

(CI).

and overlapping coding genes; hence, lncRNAs likely regulate related mRNA genes [30]. The cis and trans targets of the differentially expressed lncRNAs were both analyzed, and the most relevant GO terms associated with biological processes and molecular functions contained many important response regulators and key enzymes involved in chilling injury. (Figure 3A, Table S3). The KEGG analysis results showed that the most frequently predicted pathways were involved in oxidative phosphorylation, carbon metabolism, ubiquitin-mediated proteolysis for the *cis*-acting targets of lncRNAs; in contrast, the trans-acting targets of lncRNAs were mainly involved in glyoxylate and dicarboxylate metabolism, carbon fixation in photosynthetic organisms, and carbon metabolism, indicating their specific regulation functions in chilling injury in bell pepper (Tables S3 and S4). *Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 4 of 15 important response regulators and key enzymes involved in chilling injury. (Figure 3A, Table S3). The KEGG analysis results showed that the most frequently predicted pathways were involved in oxidative phosphorylation, carbon metabolism, ubiquitin-mediated proteolysis for the *cis*-acting targets of lncRNAs; in contrast, the trans-acting targets of lncRNAs were mainly involved in glyoxylate and dicarboxylate metabolism, carbon fixation in photosynthetic organisms, and carbon metabolism, indicating their specific regulation functions in chilling injury in bell pepper (Table S3 and S4).

**Figure 3.** The gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis of the ncRNAs targets indicated that several important enzymes were involved in the chilling injury **Figure 3.** The gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) analysis of the ncRNAs targets indicated that several important enzymes were involved in the chilling injury (CI).

The function of the differentially expressed miRNAs was also parsed. The most enriched GO terms were related to biological processes, such as response to oxygen-containing compound, response to abiotic stimulus, signal transduction, and hormone-mediated signaling pathway. However, the most relevant GO terms associated with molecular functions were protein binding, protein kinase activity, signaling receptor activity, and ATPase activity (Table S3). KEGG pathway analysis suggested that mRNAs were remarkably enriched in the pathways involved in RNA degradation and transport, peroxisome, mRNA surveillance pathway, plant hormone signal The function of the differentially expressed miRNAs was also parsed. The most enriched GO terms were related to biological processes, such as response to oxygen-containing compound, response to abiotic stimulus, signal transduction, and hormone-mediated signaling pathway. However, the most relevant GO terms associated with molecular functions were protein binding, protein kinase activity, signaling receptor activity, and ATPase activity (Table S3). KEGG pathway analysis suggested that mRNAs were remarkably enriched in the pathways involved in RNA degradation and transport, peroxisome, mRNA surveillance pathway, plant hormone signal transduction, phenylpropanoid biosynthesis, and folate biosynthesis (Table S4).

transduction, phenylpropanoid biosynthesis, and folate biosynthesis (Table S4). So far, the majority of circRNAs has not been functionally annotated [31]. To explore the potential function of differentially expressed circRNAs, GO and KEGG pathway analyses of circRNAs were performed. Our data showed that the most relevant GO terms associated with biological processes were response to temperature stimulus, response to stress, signaling, cell communication, regulation of cellular process, regulation of primary metabolic process, and signal So far, the majority of circRNAs has not been functionally annotated [31]. To explore the potential function of differentially expressed circRNAs, GO and KEGG pathway analyses of circRNAs were performed. Our data showed that the most relevant GO terms associated with biological processes were response to temperature stimulus, response to stress, signaling, cell communication, regulation of cellular process, regulation of primary metabolic process, and signal transduction (Figure 3B, Table S3). However, the KEGG pathways only included mRNA surveillance pathway and spliceosome (Table S4).

transduction (Figure 3B, Table S3). However, the KEGG pathways only included mRNA surveillance

#### pathway and spliceosome (Table S4). *2.3. Comparative Parsing of LncRNAs and mRNAs and Function Analysis of DE mRNAs*

*2.3. Comparative Parsing of LncRNAs and mRNAs and Function Analysis of DE mRNAs*  There are many differences between lncRNAs and mRNAs, including their lengths, exon numbers, open reading frames, and expression levels. The number of lncRNAs decreased with the increase of the length (<3000 nt), whereas, the length of mRNAs was distributed from 400 to ≥3000 nt and had two peaks at 400 nt and ≥3000 nt. The number of the corresponding exons of lncRNAs was There are many differences between lncRNAs and mRNAs, including their lengths, exon numbers, open reading frames, and expression levels. The number of lncRNAs decreased with the increase of the length (<3000 nt), whereas, the length of mRNAs was distributed from 400 to ≥3000 nt and had two peaks at 400 nt and ≥3000 nt. The number of the corresponding exons of lncRNAs was much less than that of mRNAs and mainly below 10; the mRNAs contained many exons, from 1 to > 30. The length of the corresponding open reading frames of lncRNAs was mainly between 50 and 300 nt, while the

much less than that of mRNAs and mainly below 10; the mRNAs contained many exons, from 1 to > 30. The length of the corresponding open reading frames of lncRNAs was mainly between 50 and 300 nt, while the length of the corresponding open reading frames of mRNAs was mainly between 100 length of the corresponding open reading frames of mRNAs was mainly between 100 and 1100 nt (Table S5). An interactive analysis of the expression of lncRNAs and mRNAs was also conducted, and their distribution on the different chromosomes was described (Figure 4). *Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 5 of 15 and 1100 nt (Table S5). An interactive analysis of the expression of lncRNAs and mRNAs was also conducted, and their distribution on the different chromosomes was described (Figure 4).

**Figure** 4 Distribution of the expression of lncRNAs and mRNAs on the different chromosomes. **Figure 4.** Distribution of the expression of lncRNAs and mRNAs on the different chromosomes.

For differentially expressed mRNAs, the most relevant GO terms associated with biological processes were defense response signaling pathway, response to auxin, response to abiotic stimulus, response to cold, and so on. However, the most relevant GO terms associated with molecular functions were protein kinase activity, transmembrane receptor protein kinase activity, protein serine/threonine kinase activity, signal transducer activity, and so on. KEGG pathway analysis indicated that the most frequently predicted pathways were involved in plant hormone signal transduction, phenylalanine metabolism, carbon metabolism, galactose metabolism, and so on. We found that many differentially expressed mRNAs were involved in chilling-related processes and could be divided into four different groups. The first group was transcription factors, including WRKY, MYB, bHLH, ERF, and NAC transcription factors, the second group was enzymes involved in bio-oxidation and oxidative phosphorylation, such as serine/threonine-protein kinase, polyphenol oxidase, catalase, peroxidase, lipoxygenase, and ATPase, the third group was cell wall metabolism, such as β-galactosidase, cellulose synthase, chitinase, pectate lyase, pectinesterase, and polygalacturonase, the fourth group was plant hormone-related processes, such as ethylene synthesis-related 1-aminocyclopropane-1-carboxlic acid synthase (ACS) and 1-aminocyclopropane-1-carboxlic acid oxygenase (ACO), abscisic acid receptor, gibberellin 2-β-dioxygenase, IAA-amino acid hydrolase, and salicylic acid-binding protein (Table S6). For differentially expressed mRNAs, the most relevant GO terms associated with biological processes were defense response signaling pathway, response to auxin, response to abiotic stimulus, response to cold, and so on. However, the most relevant GO terms associated with molecular functions were protein kinase activity, transmembrane receptor protein kinase activity, protein serine/threonine kinase activity, signal transducer activity, and so on. KEGG pathway analysis indicated that the most frequently predicted pathways were involved in plant hormone signal transduction, phenylalanine metabolism, carbon metabolism, galactose metabolism, and so on. We found that many differentially expressed mRNAs were involved in chilling-related processes and could be divided into four different groups. The first group was transcription factors, including WRKY, MYB, bHLH, ERF, and NAC transcription factors, the second group was enzymes involved in bio-oxidation and oxidative phosphorylation, such as serine/threonine-protein kinase, polyphenol oxidase, catalase, peroxidase, lipoxygenase, and ATPase, the third group was cell wall metabolism, such as β-galactosidase, cellulose synthase, chitinase, pectate lyase, pectinesterase, and polygalacturonase, the fourth group was plant hormone-related processes, such as ethylene synthesis-related 1-aminocyclopropane-1-carboxlic acid synthase (ACS) and 1-aminocyclopropane-1-carboxlic acid oxygenase (ACO), abscisic acid receptor, gibberellin 2-β-dioxygenase, IAA-amino acid hydrolase, and salicylic acid-binding protein (Table S6).

It is reported that both lncRNAs and circRNAs can interact with miRNAs through microRNA

*2.4. Construction of the Competing Endogenous RNAs (ceRNAs) Network* 
