Coding and Non-coding RNAs: Molecular Basis of Forest-Insect Outbreaks

Insect population dynamics are closely related to ‘human’ ecological and economic environments, and a central focus of research is outbreaks. However, the lack of molecular-based investigations restricts our understanding of the intrinsic mechanisms responsible for insect outbreaks. In this context, the moth Dendrolimus punctatus Walker can serve as an ideal model species for insect population dynamics research because it undergoes periodic outbreaks. Here, high-throughput whole-transcriptome sequencing was performed using D. punctatus, sampled during latent and outbreak periods, to systemically explore the molecular basis of insect outbreaks and to identify the involved non-coding RNA (ncRNA) regulators, namely microRNAs, long non-coding RNAs, and circular RNAs. Differentially expressed mRNAs of D. punctatus from different outbreak periods were involved in developmental, reproductive, immune, and chemosensory processes; results that were consistent with the physiological differences in D. punctatus during differing outbreak periods. Targets analysis of the non-coding RNAs indicated that long non-coding RNAs could be the primary ncRNA regulators of D. punctatus outbreaks, while circular RNAs mainly regulated synapses and cell junctions. The target genes of differentially expressed microRNAs mainly regulated the metabolic and reproductive pathways during the D. punctatus outbreaks. Developmental, multi-organismal, and reproductive processes, as well as biological adhesion, characterized the competing endogenous RNA network. Chemosensory and immune genes closely related to the outbreak of D. punctatus were further analyzed in detail: from their ncRNA regulators’ analysis, we deduce that both lncRNA and miRNA may play significant roles. This is the first report to examine the molecular basis of coding and non-coding RNAs’ roles in insect outbreaks. The results provide potential biomarkers for control targets in forest insect management, as well as fresh insights into underlying outbreak-related mechanisms, which could be used for improving insect control strategies in the future.


INTRODUCTION
As the most abundant and diverse animals on earth, insects play important roles in ecosystems, and their population dynamics are closely related to humans' ecological and economic environments. Outbreaks are the central focus of insect population research, long attracting the attention of both entomologists and ecologists (Hunter, 1995). Investigations on why insect outbreaks occur have been performed, considering such reasons as nutrition (Quezada García et al., 2015), natural enemies (Berryman et al., 1987), environmental factors (Veran et al., 2015;Fuentealba et al., 2017), host resistance (Elderd et al., 2013), and population interactions (Lu et al., 2010). Yet, in the same habitat, only small numbers of species may undergo an outbreak, while most of the other co-occurring species maintain low and stable population sizes within the community (Alison, 1991). This indicates that outbreaking species may respond to changes in biological and/or abiotic factors differently from that of non-outbreaking species, and the former may have an intrinsic motivation differing from the latter (Myers, 1993). To date, the physiological responses of outbreak species to some factors have been heavily researched (Berryman, 1988;Veran et al., 2015); however, the molecular basis of these responses is still unclear, which limits our understanding of the intrinsic mechanisms responsible for insect outbreaks.
Dendrolimus spp. (Lepidoptera; Lasiocampidae) are major destructive pests of conifer forests (Xiao, 1992) and Dendrolimus punctatus Walker is the most widely distributed species in China (Chen, 1990;Billings, 1991;Luo et al., 2018), making it the most important pine tree defoliator. Notably, D. punctatus is a periodic outbreaking species; once it outbreaks, the caterpillars feeding on pine needles can quickly cause the large-scale destruction of pine forests, in a phenomenon called "Fire without smoke." D. punctatus present different physiological characteristics during different outbreak stages. During the lowdensity latent period, the egg production and female ratio of D. punctatus are maintained at a certain level; during the ascending period, the egg production, female ratios, and population density increase considerably, with the degree of population aggregation also increasing; in the declining period, the body size, emergence rate, and female ratio of pine caterpillars all decrease to very low levels, while their rates of parasitism rise to a high level, additionally, their density drops rapidly to a very low level (Billings, 1991;Zhang et al., 2002). Collectively, these characteristics indicate that population density and several physiological characteristics, such as reproduction, immunity, or chemosensory system, are closely linked to the outbreak stages of D. punctatus. Outbreak of other insects also infer the importance of these physiological characteristics (Guo et al., 2011;Wang et al., 2013). Thus, D. punctatus offers a good model insect for investigating forest outbreak mechanisms, though the genetic mechanisms underlying this phenomenon remain generally unexplored.
With the advent of modern sequencing technology and acquisition of genomic information on outbreaking species, studies of insect outbreak mechanisms can now move from their ecological and physiological levels to understanding their basis at the molecular level (Kang et al., 2004;Wang et al., 2014). Accordingly, it is essential to investigate the molecular mechanisms that drive the occurrence of outbreak insects at multiple stages of an outbreak. Furthermore, biological processes are regulated not only by protein coding genes but also by non-coding RNA (ncRNA), such as microRNA (miRNA), long non-coding RNA (lncRNA), and circular RNA (circRNA) (Costa, 2005;Hansen et al., 2013;Beermann et al., 2016;Hombach and Kretz, 2016). The miRNAs-small ncRNAs of 20-30 nt-are well-known regulators of insect development and reproduction (Wei et al., 2009;He et al., 2016;Belles, 2017;Roy et al., 2018). LncRNAs are ncRNA longer than 200 bp (Amaral et al., 2008;Chen et al., 2016;Wu et al., 2016), which play important epigenetic regulating functions because they can form complexes with chromatin regulators at appropriate genomic regions in the cis and trans forms (Kretz et al., 2012;Mercer and Mattick, 2013). The circRNAs are newly discovered ncRNAs that form rings without the 5 -cap and 3 -polyadenylated tail (Wilusz and Sharp, 2013;Ling-Ling and Li, 2015), and they can originate from exons, introns, or intergenic regions Jeck and Sharpless, 2014;. Both lncRNAs and circRNAs can regulate gene expression through lncRNA/circRNA-miRNA-gene networks and they are also called "miRNA sponges" (Hansen et al., 2013). In particular, circRNA can regulate the splicing, transcription, and expression of genes, and subsequently regulate physiological characteristics Wang et al., 2017). With the development of next-generation sequencing technologies and bioinformatics, it is now possible to explore the molecular basis of dynamic population changes of outbreaking insects and to analyze the various regulators and modifiers of gene expression levels, such as miRNA, lncRNA, and circRNA, as well as the networks formed among different regulators and protein-coding genes (Yvonne et al., 2014;Li et al., 2017a;Wang et al., 2019).
In this comprehensive study, high-throughput wholetranscriptome sequencing was performed using populations of D. punctatus at low-density (latent period) and high-density (outbreak period) to systemically explore the molecular basis of different physiological characteristics during these two periods and to identify the involved regulating lncRNAs, circRNAs, and miRNAs. Furthermore, the chemosensory and immune genes and their regulating ncRNAs, which maybe closely linked to outbreak stages of D. punctatus, were analyzed in detail. This study is the first to not only identify the ncRNAs in D. punctatus but also explore the coding and non-coding RNA-based molecular mechanisms of forest insect outbreaks. These results will enhance our understanding of the unique ecological and evolutionary characteristics of the population dynamics-related mechanisms of outbreak insects.

Sampling, RNA Extraction, and RNA Quantification
Dendrolimus punctatus pupae of different densities during different outbreak stages (low-density latent period and highdensity outbreak period) were collected from Quanzhou, Guilin City, in Guangxi Province, China. The low-density insects was collected from an area where less than 10% of the trees were damaged by D. punctatus, and the insects were rare on the damaged trees; the high-density insects was collected from an area where more than 80% of the trees were damaged by D. punctatus; many insects were visible on these pine tree and so they were easy to collect. These two sampling areas were located in two neighboring villages of the same town (low-density D. punctatus: ca. N25 • 58 4.72 , E111 • 14 47.19 ; high-density D. punctatus: ca. N25 • 49 36.90 , E111 • 21 21.95 ). All the sample specimens were collected from host Masson pine (Pinus massoniana) trees (ca. 10-years old) Approximately 50 pupae were collected from each low-or high-density population, and maintained in our laboratory under a 16-h light:8-h dark photoperiod at 26 ± 2 • C and 50% ± 10% relative humidity. Newly emerged male and female D. punctatus were collected and immediately frozen in liquid nitrogen. Four groups were used: males of low-density, females of low-density, males of high-density, and females of high-density. Three biological replications were prepared for each sample.
Total RNA per sample was extracted using TRIzol (Invitrogen, Carlsbad, CA, United States). RNA degradation and contamination, especially with DNA, were checked using 1.5% agarose gels. The concentration and purity of RNA were both measured by a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States), and its integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies, CA, United States).

Library Preparation for lncRNA-Seq, Clustering, and Sequencing
From each sample, a total of 1.5 µg RNA was used as input material after first removing the ribosomal RNA (rRNA) with the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI, United States). Sequencing libraries were generated using the NEBNext R Ultra TM Directional RNA Library Prep Kit for Illumina R (NEB, United States), by following the manufacturer's recommendations. Unique index codes were added to attribute the sequences to each sample. Briefly, fragmentation was carried out using divalent cations under an elevated temperature in the NEBNext First Strand Synthesis Reaction Buffer (5×). Firststrand cDNA was synthesized using a random hexamer primer and reverse transcriptase; next, second-strand cDNA synthesis was performed using DNA polymerase I and RNase H. The remaining overhangs were converted into blunt ends through exonuclease/polymerase activities. After the adenylation of the DNA fragments' 3 ends, NEB Next adaptors with hairpin-loop structures were ligated to them to prepare for the hybridization. To select insert fragments that were preferentially 150-200 bp in length, library fragments were first purified with AMPure XP beads (Beckman Coulter, Beverly, MA, United States). Then, 3 µL of USER enzyme (NEB) was used with the sizeselected, adaptor-ligated cDNA, at 37 • C for 15 min, before running the PCR. Then, the PCR was performed using Phusion High-Fidelity DNA polymerase, universal PCR primers, and an index (X) primer. Ensuing PCR products were then purified (AMPure XP system) and library quality assessed on an Agilent Bioanalyzer 2100.
The clustering of index-coded samples was done in an acBot Cluster Generation System, with a TruSeq PE Cluster Kitv3-cBot-HS (Illumina), according to the manufacturer's instructions. After this cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and paired-end reads were generated.

Library Preparation for Small RNA-Seq, Clustering, and Sequencing
For the small RNA library preparation, 1.5 µg of total RNA per sample was used with an NEBNext R Multiplex Small RNA Library Prep Set for Illumina R (NEB), by following the manufacturer's recommendations, to which index codes were added to attribute the sequences to each sample.
Briefly, the 3 SR and 5 SR adaptors were ligated, and M-MuLV reverse transcriptase was used to synthesize the first stand chain. Then, PCR amplifications were performed using Long AmpTaq 2 × Master Mix, an SR primer, and an index (X) primer. PCR products were purified on 8% polyacrylamide gels (run at 100 V, for 80 min). DNA fragments of 140-160 bp (i.e., length of small ncRNA plus the 3 and 5 adaptors) were recovered and dissolved in an 8-µL elution buffer. Finally, all the PCR products were purified by the AMPure XP system, and library quality assessed.
The clustering of these index-coded samples was performed in a cBot Cluster Generation System, with a TruSeq PE Cluster Kit v4-cBot-HS (Illumina), according to the manufacturer's instructions. After this cluster generation, each library preparation was sequenced on an Illumina HiSeq 2500 platform and 50-bp single-end reads generated.

Quality Control
Raw data (i.e., raw reads) in the 'fastq' format were first processed using in-house perl scripts. In this step, clean data were obtained by removing from the raw data those reads containing adapters and poly-Ns, as well as any low quality reads. Additionally, the Q20, Q30, GC-content, and sequence duplication level of the clean data were calculated. For small RNA-Seq data, reads were trimmed and cleaned by removing sequences shorter than 15 nt or longer than 35 nt. Finally, at least 16.28 Gb of clean data for each lncRNA sequencing sample, with a Q30 > 92.99%, and more than 14.30 Mb of clean data for each small RNA-Seq sample, with a Q30 > 99.43%, were obtained. All downstream analyses were thus performed using only high-quality clean data.

lncRNA Identification
The transcriptome was assembled using StringTie (v1.3.1) software 1 , based on the reads mapped to the reference genome of D. punctatus (Daehwan et al., 2015;Pertea et al., 2016). The assembled transcripts were then annotated using the gff compare program 2 . The unknown transcripts were used to screen for putative lncRNAs. To sort the ncRNA candidates among the unknown transcripts, four computational approaches, namely CPC/CNCI/Pfam/CPAT, were combined and used (Lei et al., 2007;Kai et al., 2013;Liang et al., 2013;Finn et al., 2014). Those transcripts exceeding 200 nt in length, and having more than two exons, were selected as lncRNA candidates (Kelley and Rinn, 2012;Lv et al., 2013); these were further screened using CPC/CNCI/Pfam/CPAT, which has the ability to distinguish not only proteincoding from non-coding genes but also different types of lncRNAs-including lincRNA, intronic lncRNA, anti-sense lncRNA, sense lncRNA.

miRNA Identification
Use the Bowtie tools software (Langmead and Salzberg, 2012), clean reads were used as queries against the Silva, GtRNAdb, Rfam, and Repbase databases for sequence alignments, and to filter out rRNA, transfer RNA, small nuclear RNA, small nucleolar RNA, other ncRNAs as well as any repeats. The remaining reads were used to detect miRNAs through comparisons with D. punctatus genome (Zhang et al., 2020) and known miRNAs from the miRBase (v21) 3 , and novel miRNA were predicted by miRDeep3 (Friedländer et al., 2012).

CircRNA Identification and Target Prediction
CircRNA Identifier (CIRI) was used to predict the circRNAs . CIRI scans the SAM files, twice. In the first scan, CIRI detects junction reads of circRNA candidates, by using paired-end mapping and GT-AG splicing signals. Then, it rescans the SAM alignment to eliminate false-positive candidates derived from incorrectly mapped reads of homologous genes or repetitive sequences. We used the miRanda (Doron et al., 2008) and RNA hybrid (Rehmsmeier et al., 2004) tools for predicting circRNA target miRNAs and the corresponding target genes of miRNAs.

Expression Analysis
StringTie (v1.3.1) was used to calculate the fragments per kilobase of exon per million fragments mapped (FPKM) values, for both lncRNAs and coding genes, in each sample (Pertea et al., 2016). Gene FPKMs were computed by summing the FPKMs of transcripts in each gene group; FPKMs were calculated based on the length of the fragments and the read counts mapped to this fragment. The expression levels of miRNA and circRNA were each quantified using TPM (transcript per million) (Fahlgren et al., 2007).

Differential Expression Analysis
The differential expression analysis was performed using the 'DESeq' (v1.10.1) (Anders and Huber, 2012). DESeq provides statistical routines for determining differential expression in digital gene expression data, by using a model based on the negative binomial distribution. The resulting P-values are adjusted using the Benjamini and Hochberg's approach (Haynes, 2013) for controlling the false discovery rate (FDR). Genes with an adjusted P-value < 0.01 and absolute value of log 2 (foldchange) > 1 were designated as differentially expressed.

Gene Functional Annotation
The mRNAs and targets of ncRNAs were queried in a BLAST algorithm-based search against the NR, SWISSPROT, KEGG, and KOG databases (using cut-off threshold of 1e-5), from which the most similar sequence targets were selected for functional annotations. Then, the molecular functions, biological processes, and cellular components of the genes were assigned a gene ontology (GO) annotation by Blast2GO (Stefan et al., 2008). The GO enrichment analysis of differentially expressed genes (DEGs) was carried out using the 'topGO' (Alexa and Rahnenfuhrer, 2010), and KOBAS software (Xie et al., 2011) was used to test for the statistical enrichment of these DEGs in KEGG pathways.

Competing Endogenous RNA (ceRNA) Network Analysis
Differentially expressed ncRNAs and mRNAs between lowand high-density D. punctatus groups were also analyzed. The miRNA response elements were first identified, by screening the circRNA, lncRNA, and mRNA sequences. Then, the STRING database (Cytoscape 3.4.0) was used to identify protein-protein interactions among the products of these DEGs.

Validation by Real Time Quantitative PCR (qRT-PCR)
To verify our RNA-Seq results, a total of thirteen DEGs, three lncRNAs, three miRNAs, and three circRNAs, were randomly selected to detect their respective expression levels by qRT-PCR. The RNAs used for this validation were the same as those of Illumina-sequenced RNAs. PimeScipt RT reagent Kit (TaKaRa, Dalian, China) was used to synthase cDNA (for mRNA, lncRNA and circRNA) with random hexamers. Mir-X TM miRNA First Strand Synthesis and the SYBR R qRT-PCR Kit (Takara, Dalian, China) were used to synthesize the DNA template for miRNA, for which miRNA primers were designed according to the Kit's manual. The qRT-PCR was carried out with the SYBR Green PCR kit (TaKaRa, Dalian, China). Beta-actin (for mRNA, lncRNA, and circRNA) and U6 (for miRNA) served as controls. All primers used for qRT-PCR can be found in Supplementary Table S1. All qRT-PCR reactions were performed in a Roche LightCycler 480 (Stratagene, La Jolla, CA, USA), using this program: 2 min at 95 • C, 40 cycles of 20 s at 95 • C, 20 s at 58 • C, and 20 s at 72 • C; finally, a melting curve analysis (58 • C to 95 • C) was conducted to evaluate the specificity of obtained PCR products. Ct values were calculated in Roche qPCR software (v1.5.1) by applying the second derivative method. Each PCR reaction was done in triplicate. Data are presented as the log 2 (fold-change) in expression and were compared with RNA-Seq results.

Overview of Sequencing and ncRNA Predictions
In all, 209.24 Gb of clean data with a Q30 > 92.99% were obtained (Supplementary Table S2), and this data mapped onto the D. punctatus reference genome assembly, by using HISAT2 (Daehwan et al., 2015;Supplementary Table S3). A total of 55,684 lncRNAs were predicted, most of them being lincRNAs (41,324; 74.2%), followed by intronic-lncRNAs (7,384; 13.3%), antisense-lncRNAs (4,120; 7.4%), and sense lncRNAs (2,856; 5.1%) ( Figure 1A). For the lengths of lncRNA and mRNA, they had similar distribution patterns (Supplementary Figures S1A,B) but more of the mRNAs contained four or more exons (Supplementary Figures S1C,D). Most of the ORFs were shorter than 100 bp from the lncRNA, yet longer than 200 bp from the mRNA (Supplementary Figures S1E,F). Additionally, the mRNAs were characterized by more isoforms than were lncRNAs (Supplementary Figure S1G), indicating the former are more complex than the latter. Finally, the expression levels of lncRNAs and mRNAs were generally similar (Supplementary Figure S1H).
Overall, 15,698 circRNAs were predicted in D. punctatus according to CIRI, with most of these derived from exons < 3,000 bp whereas most of those derived from intergenic regions were longer than 3,000 bp ( Figure 1B). Compared with mRNAs, the expression levels of circRNAs were higher (Supplementary Figure S2).
For the miRNAs, 193.67 Mb of clean reads were acquired, with at least 14.30 Mb per sample, and their Q30 was > 99.43%. From them, any reads less than 15 nt and longer than 35 nt were first removed, the remaining reads were mapped onto the D. punctatus reference data (Supplementary Table S4). This yielded 3,738 miRNAs, of which most were between 20 nt and 22 nt ( Figure 1C); however, the length distributions of known and novel miRNAs differed (Supplementary Figures S3A,B). The first nucleotide bias analysis indicated that as the total length of a miRNA increased, so did the bias for U serving as the first nucleotide (Supplementary Figure S3C). Additionally, there was evidence of miRNA nucleotide bias at the last position for U and C (Supplementary Figure S3D).

Expression Differences and Functional Analyses of Coding RNAs Between Low and High Density D. punctatus
To analyze the molecular responses of D. punctatus during different stages of outbreaks, we compared and selected the DEGs using DESeq (Anders and Huber, 2012). The DEG numbers indicated that the differences between sexes in D. punctatus (occurring at both low and high densities) were much greater than those found between densities (both male and female) ( Table 1).
Then we focused on analyzing gene expression differences in D. punctatus from different outbreak stages. An MA plot of the differences between low-and high-density insects showed the more DEGs in males than females (Supplementary Figures S4A,B). The GO annotation of DEGs in females at low-versus high-density mainly focused on multicellular organism process, developmental process, multi-organism process, reproduction, and immune system process (Figure 2A). The top-10 enriched GO terms for DEGs between low-and high-density female insects revealed that, apart from basic metabolic processes, also included were olfactoryrelated sensory perception of smell, regulation of postsynaptic membrane potential, chromatin silencing at rDNA ( Table 2). These differences indicated that D. punctatus females from high-density and low-density populations mainly differ on their development, reproduction, immunity, and olfaction. To detect, in greater detail, further possible differences between low-and high-density insects, we analyzed the up-and down-regulated genes separately. According to their GO enrichment results, this  activity was very different between low-and high-density females of D. punctatus: the latter's up-regulated genes were involved in biological regulation, localization, response to stimulus, signaling, developmental process, multi-organism process, and reproduction ( Figure 2B); however, the down-regulated genes of these high-density females were involved in multicellular organism process, development process, multi-organism process, reproduction, reproduction process, and immune system process ( Figure 2C). Thus, the low-density D. punctatus females have better developmental, reproductive, and immune capabilities. For males of D. punctatus, a GO term analysis of DEGs between those at low-and high-density revealed they differed significantly in developmental process, multi-organism process, reproduction, and immune system process, as well as the presynaptic process involved in chemical synaptic transmission ( Figure 2D). Up-and down-regulated DEGs between lowand high-density males of D. punctatus showed different characteristics from the females. The up-regulated genes were involved in biological regulation, localization, response to stimulus, and signaling (much like females), but no developmental and reproductive processes were evidently enriched. Additionally, the presynaptic process involved in chemical synaptic transmission was represented to a higher degree in up-regulated DEGs than among the total genes ( Figure 2E). The down-regulated genes were similar to those of females, in that they were mainly involved in multicellular organismal process, developmental process, multi-organism process, reproduction, and immune system process ( Figure 2F). These results indicated the low-density males of D. punctatus put more energy into development, production, and immune functions than high-density conspecific individuals, while the males at high-density allocate more energy into basic metabolism and olfactory functions.

Functions of Non-coding RNAs in the Differential Characteristics of Low-and High-Density D. punctatus
To analyze the functions of ncRNAs (lncRNA, miRNA, and circRNA) in regulating the outbreak of D. punctatus, the differentially expressed ncRNAs were first identified. Using a cutoff of fold change ≥ 2 coupled to an FDR < 0.05, we distinguished 473 and 954 differentially expressed lncRNAs between sexes of D. punctatus from low-and high-density populations, respectively; and likewise, 138 and 169 differentially expressed lncRNAs between low-vs.  high-density D. punctatus in females and males, respectively ( Table 1). According to a Volcano plot of DEGs targeted by mRNAs and lncRNAs, the ratio values of differently expressed lncRNAs to mRNAs between low-and high-density areas (Supplementary Figures S5A,B) was higher than those between females and males (Supplementary Figures S5C,D); this indicated a stronger regulatory role of lncRNAs in generating differences found between low-and high-density insects than those between sexes. Applying the same cutoff criteria as above, 77 and 80 differentially expressed circRNAs between sexes of D. punctatus from low-and high-density areas, and 39 and 81 differentially expressed circRNAs between low-vs. high-density D. punctatus in females and males were identified, respectively (Table 1). Likewise, for differentially expressed miRNAs, 581 and 288 of them were found between female and male D. punctatus from low-and high-density populations, respectively, with 347 and 239 identified between low-vs. high-density D. punctatus in the females and males, respectively (Table 1).
Next, the targets of differently expressed ncRNAs between low-and high-density D. punctatus insects were further analyzed to explore the regulating roles of ncRNAs during this insect's outbreaks. The lncRNAs may regulate the associated mRNA genes of neighboring (cis target) and overlapping (trans targets) (Li et al., 2017a). Hence, both the cis and trans targets of lncRNAs were predicted and analyzed. For females, the cis targets of differentially expressed lncRNA between low-and high-density insects were mainly involved in developmental and multi-organism processes, reproduction, and behavior, while the trans targets were involved in biological regulation, response to stimulus, signaling, and cellular component organization or biogenesis (Figures 3A,B). The differently expressed mRNAs (Figure 2A) and the targets of differently expressed lncRNAs between low-and high-density females of D. punctatus were both involved in metabolism, development, and reproduction; however, their enriched GO terms revealed some divergence, in that mRNAs were more specific to the immunity term, whereas lncRNAs' targets were more specific to the behavior term. For males, the cis targets of differentially expressed lncRNAs between low-and highdensity insects were mainly involved in reproduction and the presynaptic process involved in chemical synaptic transmission, while the trans targets participated in multicellular organismal, developmental, and immune system processes, as well as reproduction and cell killing (Figures 3C,D). In sum, the GO terms of differently expressed mRNAs and targets of differently expressed lncRNAs from male insects at low-and high-density were coincident, indicating a primary regulator role for lncRNA in D. punctatus outbreaks.
The respective functions of miRNAs related the differing characteristics of low-and high-density populations of D. punctatus were then analyzed. In females, there were fewer target genes of differentially expressed miRNAs between lowand high-density D. punctatus than targets of all the miRNAs in the immune system process (Figure 3E). Because miRNAs play negative regulatory roles on their target genes (Ghildiyal and Zamore, 2009), these results indicated that miRNA play more regulatory roles in the low-density insect. In males, fewer target genes of differentially expressed miRNAs between lowand high-density D. punctatus populations were found than targets of all the miRNAs in both the multi-organism process and reproduction terms ( Figure 3F). This indicated that these miRNAs function in metabolic and reproductive adjustments that occur between low-and high-density males of D. punctatus.
The targeted genes of differentially expressed circRNAs between low-and high-density populations of D. punctatus were also analyzed. In females, these mainly participated in the biological process category of development, the cellular component category of synapse and cell junction, and the molecular function categories of transporter, signal transducer, and molecular transducer activities ( Figure 3G). In males, they were primarily involved in the cellular component category of synapse and cell junction, as well as molecular function categories of protein binding and signal transducer, molecular transducer, nucleic acid-binding transcription factor, and transcription factor activities ( Figure 3H). These results revealed the different targeted genes of the differentially expressed circRNAs; interestingly, the category synapse and cell junction was enriched in both comparisons.
To verify the RNA-Seq data, we selected thirteen mRNA, three lncRNA, three miRNA, and three circRNA randomly. Their respective fold-change values in expression between lowand high-density D. punctatus were tested by qRT-PCR. These results verified the reliability of the initial set of RNA-Seq results (Supplementary Figure S6).

Construction of ceRNA Networks in D. punctatus
The circRNAs, lncRNAs, and mRNAs may bind to the same miRNA response elements competitively in a regulatory network, which is referred to as a ceRNA network (Salmena et al., 2011). Figure 4A shows the total numbers of differentially expressed coding and ncRNAs in different situations. The ratios of ncRNAs targeting DEGs between low-and high-density populations are higher than those between males and females, indicating the regulatory functions of ncRNAs were stronger in the outbreak process than between sexes. Then, to derive a ceRNA network, we screened the relationships between miRNA vs. mRNA, miRNAs vs. circRNA and miRNAs vs. lncRNA (Supplementary Table S5), using the following parameters: number of miRNAs interacting with the candidate ceRNAs > 5 and P < 0.05 (Li et al., 2014). From this, a total of 55,707,481 ceRNA relationships with 41,508 nodes were obtained. With such a complex network, it was necessary to extract the critical RNAs and perform a follow-up analysis: score for all the nodes in the network were obtained and the top 10% of nodes were kept for further examination. A GO annotation of these critical RNAs indicated their involvement in developmental, multi-organism, and reproduction processes, as well as biological adhesion (Figure 4B), and these GO terms also figured prominently among the DEGs distinguished between low-and high-density D. punctatus populations. Taken together, these results indicated that an operating ceRNA regulatory network is crucial during insect outbreaks.

Differently Expressed Chemosensory and Immune Genes Between Low-and High-Density D. punctatus and Their ncRNA Regulators
The chemosensory and immune system closely related to the physiological characteristics during different outbreak stages of D. punctatus, as indicated in the previous research (Billings, 1991;Zhang et al., 2002). Our transcriptomic results also illustrated their critical roles in the outbreak process of this defoliating insect. Accordingly, we further analyzed these two gene families in low-and high-density D. punctatus. Heatmaps (Figure 5) showed that relative expression levels of chemosensory and immune genes not only differed greatly between sexes, but also between individuals at low-and high-density, especially in male insects.
Then, we investigated in depth the differently expressed chemosensory and immune genes between low-and highdensity D. punctatus and their ncRNA regulators. In females, the expression of five chemosensory genes was statistically different between low-and high-density D. punctatus ( Figure 6A): it was greater for DpunCSP8 and DpunOBP20 in low-density females, but greater for DpunOBP22, DpunOBP46, and DpunOR56 in high-density females. 14 lncRNAs contained DpuncCSP8 as their target (Supplementary Figure S7A), yet only MSTRG.120627.1 was expressed significantly higher in high-density females ( Figure 6B). Eight, six, and six lncRNAs respectively contained DpuncOBP20, DpuncOBP22, and DpuncOBP46 as targets (Supplementary Figures S7B-D), but none of these lncRNAs were expressed in a significantly different way between low-and high-density females. Six miRNAs contained DpuncOR56 as their target, of which three were significantly down-regulated in females at high-density ( Figure 6C, Supplementary Figure S7E), the reverse of DpunOR56's expression. Three immune genes, Serpin-4L, DpunGap11, and DpunTalin1, were found expressed significantly higher in high-density than in low-density females of D. punctatus (Figure 6D). No ncRNA regulator could be found for Serpin-4L. Eight lncRNAs contained DpunGap11 as a target (Supplementary Figure S7F), with only MSTRG.146927.1 expressed at a significantly higher level in the high-density than low-density females ( Figure 6E). DpunTalin1 was target of one lncRNA (Supplementary Figure S7G), one miRNA (Supplementary Figure S7H), and seven circRNAs (Supplementary Figure S7I); however, none of these ncRNAs were expressed in a significantly different way between low-and high-density females. Next, we considered the situations of low-vs. high-density males of D. punctatus, finding more differently expressed chemosensory or immune genes in males than females between low-and high-density populations of D. punctatus (Figure 7). Ten chemosensory genes were expressed at significantly differently levels between low-and high-density males of D. punctatus ( Figure 7A): DpunCSP7, Dpund17, DpunIR21, and DpunIR25a were expressed more in lowdensity males, whereas DpunCSP19, DpunOBP1, DpunPBP1, next-OR60, and DpunSNMP1 were all expressed more in high-density males. Among these chemosensory genes, an lncRNA regulator, MSTRG.120685.1, of both DpunCSP17 ( Figure 7B) and DpunCSP19 (Figure 7C), had significantly more expression in high-density males, as did another lncRNA regulator, MSTRG.120685.1, of next-IR1 ( Figure 7D); DpunPBP1 was regulated by 30 lncRNAs (Figure S8E), of which eight were expressed significantly more in highdensity males (Figure 7E). The ncRNA regulators of other chemosensory genes did not undergo a significantly different level of expression between low-and high-density males of D. punctatus (Supplementary Figure S8). Eleven immune genes were expressed significantly and differently between low-and high-density male D. punctatus (Figure 7F), of which six were down-regulated and five were up-regulated in high-density males. Among these immune genes, the lncRNA regulator, MSTRG.156626.2, of DpunCASP3 ( Figure 7G), was expressed significantly more in lowdensity males, as was the miRNA regulator, asu-miR-34-5p, of DpunKn (Figure 7H). Similarly, another miRNA regulator, unconservative_000344F_pilon_5149037, of DpunCLIP-serpin 3 (Figure 7I), had significantly greater expression in highdensity males, as did the lncRNA regulator MSTRG.226902.1 of DpunSerpin easter-L ( Figure 7J) and both lncRNA regulators, MSTRG.165722.1 and MSTRG.165710.3, of DpunPGRP2 ( Figure 7K). Finally, DpunHemocytin and DpunHemocytin-2 were both regulated by a miRNA that was expressed significantly higher in low-density males (Figures 7L,M). The ncRNA regulators of other immune genes did show any significantly expressed differences between low-and high-density males of D. punctatus (Supplementary Figure S9).

DISCUSSION
Insect species that experience outbreaks undergo extreme fluctuations in population density and consequently are capable of exerting substantial ecological and economic effects over large areas (Baltensweiler and Fischlin, 1988;Pener and Simpson, 2009). Understanding the large-scale and long-term dynamic mechanisms of outbreaking insects is therefore necessary to develop effective management strategies to mitigate these occurrences. Hence, in this study we employed whole-transcriptome sequencing and bioinformatics technology to comprehensively analyze the outbreak-related functions of coding and ncRNAs profiles and ceRNA networks in a typical forest insect pest, the defoliating moth D. punctatus.
Currently, molecular mechanism-related studies of D. punctatus outbreak characteristics remain limited to mRNAs Yang C.H. et al., 2017;Zhang et al., 2016Zhang et al., , 2018, while the study of ncRNAs in forest insects is generally limited. Our study is the first to identify and quantify expression levels of lncRNA, circRNA, and miRNA in D. punctatus. We identified 55,684 novel lncRNAs, in addition to 5,698 novel  circRNAs and 3,521 novel miRNAs. Furthermore, differentially expressed coding and non-coding RNAs were identified between the low-and high-density D. punctatus populations, for which functional annotations revealed their specific regulatory roles in this pest's outbreaks.
Firstly, the differentially expressed mRNAs related to D. punctatus population density were investigated. The differences arising between low-and high-population densities in males and females were various. In females, the DEGs were mainly involved in metabolism, development, reproduction and immunity, while in males, in addition to metabolism, development, reproduction, and immunity, the DEGs contributed to synapse junctions. Reproduction and immunity differences between low and high density of both sexes are consistent with other research in which fertility and parasitism rates differed during the non-outbreak and outbreak stages of D. punctatus (Zhang et al., 2002). The synapse is involved in neuron-related activities, such as learning and memory (Priyamvada et al., 2009;Fiumara et al., 2015). We posit that the differences in synaptic capabilities found between lowand high-density male insects reflect their varied neural activities. In this way, the molecular basis of differences in physiological traits of low-and high-density D. punctatus could be outlined.
The ncRNAs appear to play important roles in cellular functions, particularly in organismal responses to biotic and abiotic stresses (Gomes et al., 2013), and our results suggest the regulatory functions of ncRNAs in insect outbreaks are also significant. The lncRNAs participate in the regulation of different biological processes, albeit in disparate ways (Long et al., 2017;Marchese et al., 2017). Our results indicate that lncRNAs have regulatory functions during the outbreaking process of D. punctatus. In females of D. punctatus, the targets of differentially expressed lncRNAs between low-and high-density populations were mostly involved in metabolism, development, reproduction, and behavior. These functions are similar to those of their differently expressed mRNAs. Additionally, in males, the functions of the targets of differentially expressed lncRNAs were consistent with those of differently expressed mRNAs. Consequently, we conclude that lncRNAs are the primary ncRNA regulators of molecular differences that arise in this insect when its population is at low versus high density. The regulatory functions of lncRNAs during insect reproduction have been demonstrated (Wen et al., 2016;Maeda et al., 2018), and our work here indicates that reproductive regulation of lncRNAs can also affect different stages of insect outbreaks. Immune-related regulatory functions of lncRNAs are extensively researched in animals and humans (Atianand et al., 2017;Chen et al., 2017;Mowel et al., 2018), but a focus on insects has been limited. Our results demonstrated that immune-related regulatory functions of lncRNAs also operate in insects, and this may, in part, explain the resulting different immune capabilities of low-and highdensity individuals of D. punctatus.
The circRNA revealed the existence of different regulatory functions between sexes in low-and high-density D. punctatus populations, but they both participated in synapse and cell junction. This result is consistent with previous findings in which synapse regulation was determined to be a vital function of circRNAs (Westholm et al., 2014;You et al., 2015). As newly found ncRNAs, research on the functioning of circRNAs remains limited to humans and a few vertebrate species (Li et al., 2017b;FIGURE 7 | Differently expressed chemosensory and immune genes and their significant differently expressed non-coding RNA regulators between low-vs. high-density males of Dendrolimus punctatus. (A) Differently expressed chemosensory genes; (B-E) significant differently expressed non-coding RNA regulators of chemosensory genes in (A). (F) Differently expressed immune genes; (G-M) significant differently expressed non-coding RNA regulators of immune genes in (F). Xiu et al., 2019;Zhang et al., 2019). The regulatory functions of circRNAs in insects during changes in their population density deserve further research.
Being important epigenetic regulators of mRNA translation, miRNAs' involvement is implicated in many critical biological processes, such as learning, reproduction, development, metabolism, immunity, and aging (Priyamvada et al., 2009;Nan et al., 2012;Fiumara et al., 2015;Mehta and Baltimore, 2016;Ling et al., 2017;Meuti et al., 2018;Roy et al., 2018). However, miRNAs often operate in networks (Guo et al., 2016), and genome-wide analyses of miRNAs can help to reveal the complexity inherent to these networks (Li Y. et al., 2017). Our whole-transcriptome sequencing of both coding and non-coding RNAs could facilitate the identification of miRNA and miRNA-mRNA interacting pairs in D. punctatus. Functional annotations of the target genes of differentially expressed miRNAs indicated they are chiefly involved in regulating the immune and reproductive pathways during outbreaks of D. punctatus. Further functional analyses are essential to confirm the roles of miRNAs in insect outbreaks.
The "ceRNA hypothesis, " proposed by Salmena et al. (2011), describes how lncRNAs or circRNAs could inhibit miRNAs to positively regulate targeted mRNAs, by unifying the functions of several coding and ncRNAs. No ceRNAs were previously reported in D. punctatus, and so here we analyzed the ceRNA network related to D. punctatus outbreaks for the first time. Key genes analysis of the D. punctatus population density-related ceRNA network indicated that the developmental process, multiorganism process, reproduction process, and biological adhesion were components of its ceRNA network. These findings provide a theoretical basis for further investigations of the outbreak-related mechanisms of D. punctatus and other similar insects.
Besides the basic metabolic processes, our results indicated that chemosensory and immune genes play important roles in outbreak of D. punctatus. We found that chemosensory and immune genes that were differently expressed between lowand high-density D. punctatus characterized the males more than females. The lncRNAs were the main regulators for these DEGs, yet miRNA regulators were also found in several of these genes. Most of these lncRNA regulators had the same expression pattern as their mRNA targets, but some of them also showed reverse expression patterns vis-à-vis their mRNA targets, such as lncRNA MSTRG.120685.1 and DpunCSP17. These results suggest a complex regulation effect imparted by lncRNAs (Rinn and Chang, 2012). miRNA can regulate the olfactory activities of insects (Guo et al., 2018), and we also found that olfactory receptor was regulated by miRNAs in D. punctatus; but the regulatory roles of miRNA for immune gene were more comment in D. punctatus. Most of these miRNAs' expression patterns were the opposite of their mRNA targets, a finding consistent with other research (Ghildiyal and Zamore, 2009). Nevertheless, some unexpected cases were also found in our study. For example, the expression levels of miRNA regulators of DpunHemocytin and DpunHemocytin-2 as well as those these two genes were all down-regulated in the high-density population of D. punctatus insects. The reasons for these unexpected results may be inaccurate target predictions (Riffo-Campos et al., 2016) or that the miRNA acts as a subsidiary regulators of these genes, highlighting the need for further mechanistic research. Our work indicates that ncRNAs are crucial regulators during outbreaks of D. punctatus, but the mechanisms underpinning interactions between differently expressed genes and their ncRNA regulators merits further investigation.

CONCLUSION
In this study, we elucidated the coding and non-coding profiles of control and outbreak insects from a whole-transcriptome analysis. The results indicated that the physiological differences between low-and high-density D. punctatus matched well to the functions of identified DEGs. Accordingly, we may deduce that reproduction, immunity, and multi-organism processes, and the chemosensory system, are the key signal pathways involved in this insect's outbreaks. Analysis of three kinds of ncRNAs suggested that lncRNAs are the primary regulators, while miRNAs function mainly in the immune and reproduction adjustments made at different outbreak stages in D. punctatus. Further analysis of those chemosensory and immune genes closely related to outbreaking D. punctatus showed that its males contain more DEGs than do females between low-and high-density insects, and by analyzing their non-coding RNA regulators we deduced that lncRNA and miRNA play important roles during the outbreak process. To the best of our knowledge, this is the first report to have examined lncRNAs, circRNAs, miRNAs, and mRNAs expression levels and their functions in an outbreaking insect species. These findings increase our understanding of ceRNA networks and can generally aid and guide further exploration of their regulatory functions in key physiological processes required for insect outbreaks to emerge. However, the outbreak of an insect population is a complicated process: the regulatory relationships among the factors we identified above, along with initiating factors of this process, remain unclear, and so further functional-related research on these key DEGs and their non-coding regulators is necessary. Our work provides new and timely insights into the mechanisms underlying insect outbreaks that could inform control strategies to mitigate them.

DATA AVAILABILITY STATEMENT
Raw reads from sequencing are deposited in the Sequence Read Archive (SRA) database with NCBI accessiona SRR10481893-SRR10481916.

AUTHOR CONTRIBUTIONS
ZZ and SZ designed the research. SZ, SS, ZY, XK, and FL, collated the sample, performed the research and/or analyzed the data. SZ wrote the manuscript. ZZ revised the manuscript. All authors reviewed the manuscript.

ACKNOWLEDGMENTS
We thank Charlesworth Author Services for English editing of this work.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.00369/ full#supplementary-material  FIGURE S4 | MA plot of the differences between low-and high-density Dendrolimus punctatus in females (A) and males (B). FIGURE S5 | Volcano plot of mRNA and lncRNA-targeted DEGs in Dendrolimus punctatus between low-vs. high-density population in females (A), males (B), and between sexes in low-(C) and high-(D) density populations.
FIGURE S6 | Real-time PCR validation of the RNA-Seq data. The x-axis shows the RNA names, and the y-axis is their log 2 (fold change) based on the ratio of high-density and low-density insects.