miR-31-5p regulates cold acclimation of the wood-boring beetle Monochamus alternatus via ascaroside signaling

Background Survival to cold stress in insects living in temperate environments requires the deployment of strategies that lead to physiological changes involved in freeze tolerance or freeze avoidance. These strategies may consist of, for instance, the induction of metabolic depression, accumulation of cryoprotectants, or the production of antifreeze proteins, however, little is known about the way such mechanisms are regulated and the signals involved in their activation. Ascarosides are signaling molecules usually known to regulate nematode behavior and development, whose expression was recently found to relate to thermal plasticity in the Japanese pine sawyer beetle Monochamus alternatus. Accumulating evidence also points to miRNAs as another class of regulators differentially expressed in response to cold stress, which are predicted to target genes involved in cold adaptation of insects. Here, we demonstrate a novel pathway involved in insect cold acclimation, through miRNA-mediated regulation of ascaroside function. Results We initially discovered that experimental cold acclimation can enhance the beetle’s cold hardiness. Through screening and functional verification, we found miR-31-5p, upregulated under cold stress, significantly contributes to this enhancement. Mechanistically, miR-31-5p promotes production of an ascaroside (asc-C9) in the beetle by negatively targeting the rate-limiting enzyme, acyl-CoA oxidase in peroxisomal β-oxidation cycles. Feeding experiments with synthetic asc-C9 suggests it may serve as a signal to promote cold acclimation through metabolic depression and accumulation of cryoprotectants with specific gene expression patterns. Conclusions Our results point to important roles of miRNA-mediated regulation of ascaroside function in insect cold adaptation. This enhanced cold tolerance may allow higher survival of M. alternatus in winter and be pivotal in shaping its wide distribution range, greatly expanding the threat of pine wilt disease, and thus can also inspire the development of ascaroside-based pest management strategies. Supplementary information The online version contains supplementary material available at 10.1186/s12915-020-00926-w.


Background
Low temperatures pose a major challenge for insects in a highly seasonal temperate zone environment. To survive cold winters, insects have evolved strategies such as freeze tolerance by having the ability to survive internal ice formation and freeze avoidance by supercooling [1][2][3][4]. These strategies require diverse physiological changes including metabolic depression, accumulation of cryoprotectants, ion and water balance, and production of antifreeze proteins [4][5][6][7]. While much is now known about physiological plasticity involved in insect cold hardiness, signaling processes involving the acquisition of cold tolerance are still poorly understood.
As autumn progresses, insects respond to environmental signals such as shortening day length and decreasing temperatures with dramatic physiological changes that trigger seasonal cold acclimation and initiation of diapause [3,6]. Signaling pathways involving hormonal control [8], energetic regulation [9], and insulin signaling [10,11], which are associated with entry into diapause, have received much attention. As diapause and cold tolerance are often inextricably linked, these signaling pathways also provide some insights for insect cold tolerance. Mechanistically, insects first sense the cold signals, then activate the signaling pathways and reprogram the expression of cold-responsive genes that ultimately leads to metabolite changes [12,13]. Sometimes, these thermal-dependent metabolites are not only a passive target of cold signaling, but could also act as a secondary chemical signal that functions like hormones to initiate cold signaling and ultimately regulates cold-responsive gene expression [9,14]. Although the cryoprotective roles of these metabolites in insect cold adaptation have been widely recognized [6,15,16], our understanding of their potential roles in signaling is still limited.
Cold stress-relevant changes in gene expression have been documented for a variety of insect species such as bumble bees, mosquitoes, and moths [17][18][19]. However, their epigenetic regulatory networks remain largely unknown. MicroRNA (miRNA), small non-coding RNAs of ∼ 21 to 24 nucleotides in length, are a ubiquitous epigenetic regulator at the post-transcriptional level in organisms. Accumulating evidence from sequencing data have also uncovered several subsets of miRNAs that are differentially expressed in response to cold stress and are predicted to target some critical genes for physiological changes, suggesting miRNAs are potentially involved in cold adaptation in insects [20][21][22], while the mechanism by which miRNAs regulate cold adaptation remains unexplored.
The Japanese pine sawyer beetle, Monochamus alternatus Hope (Coleoptera: Cerambycidae), is not only a serious wood-boring pest in itself, but is also the main vector of the invasive pinewood nematode (PWN), Bursaphelenchus xylophilus, which is native to North America and causes destructive pine wilt disease over a wide temperate range including Japan, Korea, China, and Europe [23]. This beetle is a freeze-avoiding insect and has developed several mechanisms such as facultative diapause, pre-emptive acclimation, and seasonal and developmental variation in cold hardiness to mitigate cold winters [24,25]. During both natural and experimental cold acclimation, supercooling point (SCP), which is the lowest temperature at which an abrupt rebound on the thermal curve occurs due to the release of the latent heat of ice crystallization, significantly drops and survival rate increases in this beetle, indicating enhancement of cold hardiness [24], but mechanisms underlying this process have not been well explored.
Intriguingly, ascarosides, widely conserved small signaling molecules in nematodes that contain dideoxy sugar ascarylose and a fatty acid side chain [26], were also found to be produced by this beetle [23]. Moreover, the larval beetle, which is unable to vector PWN, produced asc-C9 that harbors 9 carbons on the side chain and is known as ascr#10 in nematodes [27], and it was shown this ascaroside retards the larva's development to synchronize its life cycle with the nematodes its adult vectors, maintaining their mutualistic relationship [23]. It is noteworthy that abundance of this small molecule correlated positively with low temperature [23]: (i) In the field, asc-C9 is mostly produced by the beetle during the last larval instar, which experiences low-temperature stress during cold winters; (ii) in the laboratory, the amount of asc-C9 increases during an experimental cold acclimation process (4°C). This correlation between ascaroside abundance and low temperature thus suggests a potential role of ascarosides in cold acclimation of this beetle.
Here, we hypothesized that miRNAs and ascarosides might be involved in insect cold signaling. To test this hypothesis, we first demonstrated cold acclimation can enhance the beetle's cold hardiness. Subsequently, we screened miR-31-5p, a miRNA whose expression was upregulated under cold stress, and validated that it significantly facilitates the observed cold tolerance. Mechanistically, miR-31-5p negatively targets the rate-limiting enzyme, acyl-CoA oxidase (ACOX1), in the peroxisomal β-oxidation cycle that shortens fatty acid chains of ascarosides. Feeding experiments with synthetic asc-C9 and gene expression analysis suggested the ascaroside might act as a signal to regulate the beetle's cold tolerance by metabolic suppression and cryoprotectant accumulation. These results demonstrate the important roles of miRNAs and ascarosides in insect cold acclimation.

Cold acclimation enhances cold hardiness
To confirm whether low-temperature acclimation can enhance the cold hardiness of Monochamus alternatus beetles, we tested both SCPs and temperature causing 50% mortality (LT 50 ). We made daily SCP measurements during 10 days' acclimation to 4°C, and measured mortality of larvae after 10 days of acclimation by exposing them to a series of low temperatures (− 20°C, − 15°C, − 10°C, and − 5°C) for 6 h, then determining LT 50 for each treatment by probit regression. The SCPs of lowtemperature-acclimated beetles dramatically decreased compared to those of control beetles after 4 days (Fig. 1a). Further, LT 50 of acclimated beetles was 4.9°C lower than that of controls ( Fig. 1b; LT 50 : control: − 11.7 ± 0.5°C, acclimated: − 16.6 ± 0.4°C; extra sum-of squares F test: F 1, 20 = 57.260, P < 0.001). Collectively, these results suggest that acclimation to 4°C enhances cold hardiness in this beetle.

miRNAs involved in the cold acclimation
We performed high-throughput sequencing of small RNAs of the acclimated and control last-instar beetle larvae (Fig. 2a). About 202-231 known miRNAs and 75-90 novel miRNAs were identified from 3 replicate control libraries and 3 replicate cold acclimated libraries, respectively (Additional file 1: Fig. S1A). Through comparative analysis, we found a total of 20 miRNAs to be differentially expressed (DE) between the acclimated and control beetles (Fig. 2b, Additional file 1: Fig. S1B, Additional file 2: Table S1; P value < 0.05 and fold change > 1.5). Their expression patterns were subsequently validated by quantitative PCR (qPCR) (Fig. 2c, Additional file 1: Fig. S2). Among the 8 markedly upregulated miRNAs expressed during low-temperature acclimation, miR-31-5p exhibited the greatest extent of expression upregulation (RPM in RNA-seq: 2.91-fold changes; qPCR: 2.58-fold changes) as well as the highest abundance (RPM in RNA-seq: 1369.30 ± 216.53 in control and 3984.73 ± 143.16 in the acclimated beetles) (Fig. 2c, Additional file 2: Table S2). This expression pattern indicated a potential role of miR-31-5p in cold acclimation of M. alternatus beetles.

Screening for gene targets of miR-31-5p
We conducted comparative transcriptome analysis between the acclimated and control last-instar beetle larvae by RNA-seq. In total, we found 15,486 differentially expressed transcripts (DETs) (fold change > 1.5, P value < 0.05), accounting for 25.8% of the total number of transcripts (59,971) in the reference transcriptome, of Fig. 1 Cold acclimation enhances beetle cold hardiness. a Supercooling points (SCPs) of last-instar larval beetles were significantly lower after 4 days of acclimation at 4°C than those of beetles held at 25°C (n = 20 for each treatment per day). b Last-instar larval beetles acclimated to 4°C for 10 days had significantly reduced mortality relative to those held at 25°C (n = 10 for each treatment). Lethal temperature causing 50% mortality (LT 50 ) was determined by probit regression. The data in a and b are shown as mean ± s.e.m. Two-tailed Student's t test was used to test significant differences. *P < 0.05; **P < 0.01; ***P < 0.001. AC, low-temperature acclimation (4°C); Control (25°C). The last larval beetles used in a and b were collected at 20-22 Oct. 2016 and 20-25 Oct. 2019 which 8439 were upregulated and 7047 were downregulated (Fig. 3a, Additional file 3: Dataset S1). KEGG analysis of these DETs showed 21 clusters were enriched (P < 0.05) (Fig. 3b, Additional file 3: Dataset S2). Among these clusters, five pathways including biosynthesis of secondary metabolites, glycine, serine and threonine metabolism, pentose and glucuronate interconversions, cysteine and methionine metabolism, and neuroactive ligand-receptor interaction were significantly enriched at the level of P < 0.01 (Fig. 3b). Among these pathways, 197 significantly enriched genes were designated as candidate genes (Additional file 3: Dataset S2). We employed the algorithms miRanda [28] and RNAhybrid [29] to predict the target genes of miR-31-5p and ultimately identified six putative target genes, namely aminomethyltransferase (amt), hairy, succinate-CoA ligase (suc), alpha actinin (α-actinin), acyl-CoA oxidase 1 (acox1), and proline dehydrogenase (pdh).
Furthermore, we predicted the target sites of mal-miR-31-5p in acox1 orthologs from M. alternatus and other model animals including Homo sapiens, Mus musculus, Drosophila melanogaster, and Caenorhabditis elegans (Additional file 1: Fig. S3A). The results showed that most (except for Mal-acox1b, Cel-acox1.2, and Cel-acox1.3) of the selected genes harbored the target site of mal-miR-31-5p in either coding sequence (CDS) regions (M. alternatus, D. melanogaster, and C. elegans) or 3′ UTR regions (H. sapiens and M. musculus) (Additional file 1: Fig. S3B). Though the relative position of these target sites in CDS was variant (Additional file 1: Fig. S3C), it is at least suggested the miR-31-5p-acox1 axis may be conserved across the Animal Kingdom.
The presence of a putative binding site in the coding sequences suggests a possible direct regulation between miR-31-5p and acox1 mRNA in M. alternatus (Fig. 4d, Additional file 1: Fig. S3B). We tested this hypothesis by performing luciferase assays in Drosophila S2 cells. Luciferase activity from the constructs containing the predicted target site (wild type, WT) from coding sequences (CDS) of acox1 decreased significantly after injection of agomir-31-5p compared to those injected with agomircontrol. Furthermore, luciferase activity of the constructs harboring a seed sequence-mutated binding site had no difference after injection of agomir-31-5p compared to those injected with agomir-control ( Fig. 4d; one-way ANOVA with Tukey's HSD: F 3, 23 = 6.555, P = 0.013). Collectively, these results suggested miR-31-5p directly targeted acox1 in a negative manner.
To further characterize specific gene expression patterns after asc-C9 treatment, we conducted comparative transcriptome analysis between the asc-C9-fed and control last-instar beetle larvae by RNA-seq. In total, we found 1601 differentially expressed transcripts (DETs) (fold change > 1.5, P value < 0.05), of which 446 were upregulated and 1155 were downregulated after asc-C9 treatment (Additional file 1: Fig. S6A, Additional file 3: Dataset S3). Gene Ontology (GO) analysis of these DETs showed "metabolic process" with 100 genes accounting for 54.05% of the annotated DETs, was the most frequent GO term in the category of biological process (Fig. 6c, Additional file 3: Dataset S4), indicating the dramatic effect of asc-C9 on beetle's metabolism. Among these enriched 100 genes in metabolic process, 78 were downregulated and 22 were upregulated (Additional file 3: Dataset S4). In addition, most genes, which were classified into the 12 significantly enriched sub-GO terms (P < 0.05), were downregulated (percentage of downregulated Fig. 6 Asc-C9 acts as a cold signal to promote beetle's cold acclimation. a The effect of asc-C9 on the supercooling points (SCPs) of the larval beetles (n = 20). b Probit regression between mortality rate and temperature of asc-C9-fed and control beetles (n = 10 for each treatment). c Distribution of Gene Ontology (GO) terms with > 1% of annotated genes in the category of biological process enriched from the DETs between asc-C9-fed and control beetles. d Simplified schematic diagram of the synthesis and metabolism of these three cryoprotectants modified from Barth et al. [33]. The genes in red are significantly upregulated and those in blue are significantly downregulated, and those in black have no changes in expression after treatment of asc-C9. Genes legends: Tpp, trehalose-phosphate phosphatase; Gp, glycogen phosphorylase; A1e, aldose-1-epimerase; Akr, aldo-keto reductase; Gpi, glucose-6-phosphate isomerase; Pfk, 6-phospho-fructokinase; Suc, succinyl-CoA ligase; Eno, enolase; Idh, isocitrate dehydrogenase; FBPase, fructose-1,6-bisphosphatase; Pgm, phosphoglycerate mutase. G1P, glucose-1-phosphate; GAP, glyceraldehyde-3-phosphate. The data were shown as mean ± s.e.m. Two-tailed Student's t test was used to test significant differences. *P < 0.05; **P < 0.01; ***P < 0.001. C9, asc-C9-fed beetles (4°C); Control beetles (25°C). The last larval beetles in a and b were collected at 6-8 Oct. 2018 and 20-25 Oct. 2019, respectively genes in each sub-term ranges from 50 to 100%) (Additional file 1: Fig. S6B; Additional file 3: Dataset S5), suggesting a strong metabolic depression.

Discussion
In this study, we found a novel function of thermaldependent ascarosides that can serve as a secondary cold signal to facilitate cold acclimation in M. alternatus. Mechanistically, upon sensing low temperature, asc-C9 production was promoted by miR-31-5p through negatively targeting a rate-limiting enzyme, ACOX1, in the peroxisomal β-oxidation cycle for shortening of the side chains of the ascarosides. Asc-C9 appears to further regulate cold physiology by altering gene expression linked to metabolic depression and cryoprotectant accumulation, ultimately enhance cold hardiness in this beetle (Fig. 7). These results reveal a novel molecular mechanism by which a miRNA-mediated small signaling molecule promotes cold acclimation in an insect.
Since the first report of the involvement of ascarosides in dauer formation in C. elegans [34], ascarosides have been shown to be multi-functional signaling molecules regulating male attraction, hermaphrodite repulsion, olfactory plasticity, aggregation, and lipid metabolism in nematodes [35,36]. In addition, ascarosides were also shown to coordinate the life cycles between the pinewood nematode and its vector beetle, thus maintaining the invasive symbiosis and promoting the dispersal of pinewood nematodes [23]. Here we further demonstrated the first reported function of these small molecules in cold acclimation in M. alternatus.
Insect cold hardiness is often achieved by physiological changes [5,6], which are mostly linked to specific gene expression patterns [33,37]. Here, we characterized a gene expression pattern of metabolic depression and accumulation of cryoprotectants after asc-C9 treatment. Metabolic depression is essential for conserving energy stores and thus could enhance resistance to cold stress [9,14,16]. Low molecular weight polyols and sugars including glycerol, trehalose, and sorbitol have been reported as cryoprotective agents in many insects [4,8]. The functions of these compounds include both colligative depression of supercooling point (SCP) and non-colligative stabilization of membranes and proteins [6]. The gene expression patterns of metabolic depression and cryoprotectant accumulation documented here thus are typical for cold hardy insects, and further supported the important role of ascaroside in enhancement of the beetle's cold hardiness. Since low winter temperature is a limiting factor to the distribution and potential dispersal areas in this beetle [38], the enhanced cold tolerance allowing higher survival of M. alternatus in winter and might have shaped its wide distribution range that occupies a large latitudinal span across Indomalayan and Palearctic realms, which has dramatic implications for the future spread of PWN and pine wilt disease [39]. It is well known that chemical signals, especially the glandular-produced hormones like juvenile hormone (JH) and ecdysone, are largely involved in cold stress responses in insects [14]. Additionally, several metabolites such as soluble sugars, the tetrapyrrole intermediate Mg-protoporphrin (Mg-ProtoIX), and ROS have also been suggested to be chemical signals involved in cold adaptation in plants [40][41][42][43]. However, metabolitederived stress signals are almost completely unknown in insects. Here, we provide experimental lines of evidences that an ascaroside is likely to be a metabolic signal for cold acclimation, and this is supported by the following two reasons: First, asc-C9, a small molecule derived from metabolites, increased in abundance in response to decreased temperatures, illustrating the feature that metabolic signals have thermal-dependent plasticity; second, feeding with trace synthetic asc-C9 can significantly altered gene expression patterns of metabolic depression and cryoprotectant accumulation in our study, and retards the development of the beetle by reducing ecdysone, slowing it down for pupation [23]. Lastly, many analogs of ascarosides with different side chain lengths have been identified and quantified in insects from several orders (unpublished data from the authors), implying these small signaling molecule might widely participate in cold adaptation in other insects. miRNA-mediated stress adaptation has been documented in several studies [20]. Here, we obtained the first miRNA library of M. alternatus beetles and described a special miRNA-mediated cold adaptation response where miR-31-5p regulates insect cold acclimation via regulation of the production of a cold signal (ascaroside) rather than directly targeting the stress responsive genes or their transcription factors [44]. Specifically, upregulated miR-31-5p under cold stress promoted the accumulation of asc-C9 by negatively targeting acox1, the first and rate-limiting enzyme gene in peroxisomal β-oxidation, resulting in enhanced cold hardiness by metabolic modification. Targeting a signal makes the miRNA more effectively promote cold acclimation compared to that of targeting a single responsive gene, as the signal can synchronously coordinate a series of pathways underlying this process [45].

Conclusions
Overall, we addressed an important and under-explored aspect of insect cold acclimation that a miRNAmediated ascaroside serves as a metabolic cold signal and promotes insect cold hardiness via shaping the expression of genes involved in metabolic depression and accumulation of cryoprotectants. Effective cold adaptation enables both the vector insect and pinewood nematode to have a wide distribution range and leads to further range expansion of pine wilt disease, which may be exacerbated by climate change. These results thus can also inspire the development of miRNA or ascarosidebased pest management for this important quarantine pest. Furthermore, the possible involvement of ascarosides in cold adaptation of other insects is worthy of further exploration.

Beetle collection and rearing
Final (5th) instar larval Monochamus alternatus were collected from host trees of Pinus massoniana, in Fuyang, Zhejiang province, in late autumn (i.e., 1 October-15 November) from 2016 to 2018, when the larval beetles begin to face lower temperatures and initiate cold acclimation, but have not yet completed cold acclimation nor entered into diapause. Additionally, the amount of asc-C9 of these larval beetles peaks during this period [23]. Thus, this instar is the ideal stage for the experimental cold acclimation and ascaroside analysis. Beetles were brought to the laboratory and reared on artificial diet (200 g of sawdust, 15 g of agar, and 550 ml of distilled water) in a 10-ml tube at a 12:12 h light:dark (L:D) cycle at 25°C placed in a climate chamber. Fresh diet was provided every week. These beetles were used for all of the subsequent experiments unless otherwise noted.
A previous study has shown that M. alternatus has seasonal variation in super cooling points [24]. Meanwhile, our study proved the amount of ascaroside (asc-C9) produced by the last-instar beetle larvae significantly changes with cold acclimation in the laboratory (artificial cold acclimation) and at different time points depending on the collection date during late autumn (i.e., natural cold acclimation) in field (Additional file 1: Fig. S4A, B). To eliminate the effect of these variations, we standardized the beetle larvae used in each experiment by only using beetles that were field collected at the same time.

Experimental acclimation, measurement of SCPs, and mortality rate
To examine the effect of low temperature on the beetle's cold hardiness, 200 final (5th instar) larval beetles with weights ranging from 300 to 500 mg were individually placed in a climate chamber at 4°C for experimental cold acclimation [23,24], with 200 individuals at 25°C for non-acclimation as control. Except for the temperature, the other parameters were the same as their rearing conditions. SCP is the lowest temperature at which an abrupt rebound on the thermal curve occurs due to the release of the latent heat of ice crystallization in an insect [24]. Twenty individuals were sampled for their super cooling point (SCP) measurement from both cold acclimation experiments and asc-C9-feeding experiments every day for 10 days. SCPs were measured using an insect SCP tester (SUN-V, Pengcheng Electronics, Beijing, China) with a refrigerator (BC/BD 326C, Haier, China). Specifically, the larvae was placed in a 1.5-ml tube and attached to a thermocouple of the tester, then fixed and enclosed by a piece of cotton [50]. The thermocouple with the larva was placed into a freezing chamber that was cooled gradually at a rate of ca. 1°C/min during measurements. The larva's body temperature decreased in a nonlinear way, the SCP being represented on the recorder by a sudden spike in the temperature of the thermocouple. Data were automatically recorded using customized software (Pengcheng Electronics, Beijing). Twenty replicates were used for each treatment.
Mortality experiments were conducted as previously described [24,50] with modification. Specifically, larval beetles from each treatment were exposed to one of four different cold constant temperatures (− 20°C, − 15°C, − 10°C, and − 5°C) in refrigerators for 6 h, respectively.
After that, all larvae were placed in a climate chamber (25°C, light:dark = 12:12 h) to recover for 1 day. Death was assessed by the absence of mandibular or body movement when larvae were stimulated with a needle. Finally, the lethal temperature causing 50% mortality (LT 50 ) for each treatment was determined by probit regression. Nine to twelve individuals were tested in each treatment. Each treatment had three replicates.

High-throughput sequencing of small RNAs and RNA-seq
To examine the miRNAs and mRNAs involved in the cold acclimation and ascaroside feeding experiments, the acclimated (4°C)/control (25°C) and asc-C9-fed/control larval beetles (weight: 400 ± 20 mg) were respectively collected after 10 days of treatment for high-throughput sequencing of small RNAs and RNA-seq, respectively. Three replicates were used for each treatment.
Total RNA from the whole body was extracted by TRIzol Reagent (Invitrogen, USA). Small RNA and mRNA libraries were constructed using a TruSeq small RNA sample preparation kit (Illumina, USA) and RNA Library Prep Kit (Illumina, USA), respectively. The libraries were sequenced on the platform of Illumina HiSeq2500 at Bionova Biotech Company (Beijing). Raw data of the libraries were deposited in the NCBI Sequence Read Archive (accession number: PRJNA522884).
After quality control, bioinformatics analyses were performed. For miRNAs, the remaining reads were clustered based on sequence similarity and dominant reads were analyzed by bowtie (version 2.3.4) [51], against Genbank and Rfam to discard rRNA, scRNA, snoRNA, snRNA, and tRNA. Subsequently, known miRNAs were identified by mapping to the miRNA sequences from insects in miRBase 22.1, and novel miRNAs were predicted using miRDeep2 (version 2.0.1.2) [52]. Differential expression analysis of miRNAs was performed using the DEGseq method [53], and differentially expressed miRNAs were finally selected with fold change > 1.5 and P < 0.05. With respect to RNA-seq, the filtered data were assembled by Trinity (version 2.6.5) [54], to obtain reference transcripts due to the reference transcriptome (PRJNA374773). The differentially expressed transcripts (DETs) were analyzed using Cuffdiff software (version 7) [55], and differentially expressed mRNAs were finally selected with fold change > 1.5 and P < 0.05. GO and KEGG enrichments were performed with cut-off P values of 0.01 and 0.05, respectively (Additional file 3: Dataset S6).
Predicted targets of miR-31-5p miRanda and RNAhybrid were employed to predict the targets of miR-31-5p in the significantly enriched genes (P < 0.01; Additional file 3: Dataset S2) by KEGG analysis of cold related DETs (Additional file 3: Dataset S1). An alignment score greater than 154 was used in miRanda [28]. The parameters of RNAhybrid were set as 1 hit per target, free energy threshold was − 25 kcal mol − 1 , helix constraint was 2 to 7, max bulge loop length was 2, max internal loop length was 6, and no G:U in seed is included [56]. The overlapped genes predicted through both methods were selected.
qPCR for mRNA and miRNA Total RNA was isolated from the whole body by TRIzol (Invitrogen, USA). Two micrograms of total RNA from each sample was reverse-transcribed using a miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, China) and the FastQuant RT Kit with gDNase (Tiangen) according to the manufacturers' instructions for miRNA and mRNA, respectively. The relative expression of miRNAs and mRNAs was respectively quantified by the miRcute miRNA qPCR Detection Kit (Tiangen) and a Real Master Mix Kit (SYBR Green) (Tiangen) with an MX3000P Thermal cycler (Stratagene, USA) instrument. The programs of thermal cycling were set following the manufacturer's instructions for miRNA: 95°C for 15 min, 40 cycles at 94°C for 20 s, 60°C for 34 s; mRNA: 95°C for 1.5 min, 40 cycles at 94°C for 15 s, 55°C for 20 s, 68°C for 60 s. To check for the specificity, melting curves were analyzed for each data point. The primers used here are included in Additional file 3: Dataset S7. U6 snRNA and β-actin were respectively selected as reference genes from 5 miRNAs and 1 snRNA, and 3 housekeeping genes examined following expression stability analysis using the Excel 2007 add-in NormFinder software (version 0.953) [57]. The expression values were calculated using the 2 −ΔΔCT method [58]. Six biological replicates and three technical replicates were performed for statistical analysis.

Western blot
Polyclonal antibody against ACOX1 was developed through injection of the KLH protein-coupled polypeptide "CTQKSPLNAEQVNKS" into rabbits, which is specific for ACOX1, but lacks in ACOX1b and ACOX3 [59]. After 45 days, the antibody in serum was collected and purified to 4.8 mg/mL with 90+% purity (Beijing Protein Innovation, China), and the specificity of the antibody was evaluated (Additional file 1: Fig. S7). For protein analysis, the whole body of a larval beetle was homogenized, then lysated in RIPA Lysis Buffer (Beyotime, China) with 1% proteinase inhibitor cocktail (Thermo Scientific, USA) on ice for 2 h. The protein concentration was measured by BCA Protein Assay Kit (Thermo Scientific). After denaturation in boiling water for 5 min, the proteins were separated using Mini-Protean TGX (Bio-Rad, USA) and transferred onto polyvinylidene difluoride membranes (Millipore). Membranes were blocked for non-specific binding using blocking buffer (Thermo Scientific), then incubated with the primary antibody (rabbit anti-ACOX1 serum or rabbit anti-a-TUBULIN, 1:5000) in blocking reagents at room temperature for 2 h. After incubation, the membranes were washed using TBST three times (15 min/each). The washed membranes were incubated with secondary antibody (anti-rabbit IgG, CST, 1:10,000) for 2 h at room temperature and were washed again. Pierce ECL Western blotting Substrate (Thermo Scientific) was added onto the membranes, and the protein stripes were imaged and analyzed on a gel image analyzing system. Three replicates were used for each treatment. The original Western blot figure for each experiment was stored in Additional file 1: Fig. S8.

Injection of dsRNAs, antagomir, and agomir
Double-stranded RNAs (dsRNAs) were synthesized using primers ligated with T7 RNA polymerase promoter sequences (Additional file 3: Dataset S7, S8) following the protocol of the T7 RiboMax Express RNAi System Kit (Promega, Madison, USA). Twenty larval beetles were injected with dsAcox1 or dsGFP (control) every 72 h, respectively. All the injections were performed using a 7000-series modified microliter syringe (Hamilton, Bonaduz, Switzerland) with 30 μg of dsRNA into the intersegmental membrane between the third and fourth abdominal segments, where the dsRNAs could be transferred to the target genes through the hemolymph. After rearing larvae for 7 days at standard conditions, effectiveness of the knock-down was evaluated by qPCR and Western blot (Additional file 1: Fig.  S7), and the amount of ascaroside (asc-C9) and SCPs were also measured.
To determine the effect of miR-31-5p-acox1 axis on the SCPs and LT 50 of beetles, we performed rescue experiments. We first injected the beetles with antagomir-31-5p or antagomir-control, and then rescued the antagomir-31-5p-treated beetles by RNAi by injecting dsAcox1 or exposing them to experimental lowtemperature acclimation at 4°C (AC) for 7 days. Finally, the SCPs were measured and compared between beetles from control, antagomir-31-5p, dsAcox1, and AC groups. Eleven to seventeen replicates were used for each treatment.

Extraction and analysis of ascaroside
Ascarosides from beetles and nematodes were extracted as previously described [23,32,60] with modification. For the larval beetles, it was first confirmed that nematodes were absent from the beetle larvae as determined through microscopic observation. To extract the surface ascarosides, the whole body of each beetle was extracted for 30 min in 1 ml of ethanol [23]. The extracts were purified through filter paper, lyophilized in a vacuum centrifugal concentrator (RVC2-18CD, Sigma), and resuspended in 150 μL of 50% (vol/vol) methanol in water. Subsequently, the liquids were transferred into a liquid chromatography-mass spectrometry (LC-MS) vial and stored at 4°C until chemical analysis. There were eight replicates for each experiment.
To get the exo-metabolome extracts of C. elegans, 100 ml of liquid cultures inoculated with 10,000 mixed stage nematodes was kept at 20°C and 150 rpm. Concentrated Escherichia coli OP50 bacteria pellet from an overnight culture in Lysogeny Broth (LB) medium at 37°C and 170 rpm was provided as food. After 6 days, the liquid cultures were centrifuged at 800g for 2 min. After removing the worms at the bottom, the supernatant was centrifuged again at 2700g for 10 min and the harvested supernatants were filtered, lyophilized, and extracted in methanol for 12 h. The extracts were again filtered and lyophilized, then reconstituted in 150 μL of 50% (vol/vol) methanol for HPLC-MS analysis.
The ascaroside analysis was performed using an Agilent 1290 Infinity LC/6460 triple quadrupole MS, equipped with a ZORBAX SB-Aq column (2.1 × 100 mm, 3.5 μm, Agilent, USA) following the previous methods [23]. Asc-C9 was monitored in the negative ion mode (the characteristic product ion monitored was m/z 73 and m/z 303.2). The concentration was then determined using standard curves. Asc-C9 was synthesized following the methods of Zhao et al. [23].

In vitro luciferase reporter gene assays
The 429-bp sequence (Additional file 3: Dataset S8) of CDS surrounding the predicted miR-31-5p target sites (CACTAGGTTGATTTTTCTTGCCG) was cloned into the psiCHECK-2 vector (Promega, Madison, USA) using an in-fusion cloning system with overlapping primers (Additional file 3: Dataset S7). The fragments of plasmid vector and gene (cDNA) were amplified using Phusion High-Fidelity DNA polymerase (NEB) in a 50 μl reaction, using the following procedure: 98°C for 30 s, 30 cycles at 98°C for 5 s, 55°C for 20 s, and 72°C for 20 s (the fragments) or 2 min (psiCHECK-2 vector plasmid), then extension at 72°C for 10 min. The products of vector were treated by DpnI enzyme to digest the remaining circular plasmids. Afterward, the PCR products of psiCHECK-2 vector and CDS were ligated using a Gibson Assembly Cloning Kit (NEB) at 50°C for 2 h and transfected into DH5α competent cells for positive clone screening. Site mutagenesis of the nucleotides complementary to the miR-31-5p seed sequence in the recombinant plasmid was performed using the Fast Mutagenesis System (TransGen, China) with mutation from "TCTTGCC" to "CACCAAA." All of the recombinant plasmids were extracted using a Plasmid Mini Kit (Qiagen, German). Five hundred microliters of 8 × 10 4 cells/mL − 1 Drosophila S2 cells were seeded into each well of a 24-well plate with complete Schneider's Drosophila medium. After 24 h, Drosophila S2 cells were co-transfected with 0.4 μg of the luciferase reporter plasmids (WT or MT) and 100 nmol of the agomir-31-5p or agomir-control (Riobio) using 1.5 μL Attractene Transfection Reagent (Qiagen). After 24 h transfection, the activities of the firefly and Renilla luciferases were measured with a Dual-Luciferase Reporter Assay (Promega) using a luminometer (Promega). Eight replicates were done for each treatment.

RNA immunoprecipitation assay
Polyclonal antibody against AGO1 was developed according to the methods of anti-ACOX1 antibody. The polypeptide used is "DITHPPAGDSRKPSC," and the specificity of the antibody was also evaluated by Western blot. RIP assay was performed using the Magna RIP Quad kit (Millipore, USA) with the polyclonal antibody against AGO-1 protein (GenScript, Nanjing, China) and a negative control IgG (Cell Signaling Technology). Specifically, six last-instar larval beetles were injected twice with agomir-31-5p or agomir-control, respectively. Seven days following injections, the beetle was homogenized in ice-cold RIP lysis buffer, then stored at − 80°C overnight. Next, 50 μL magnetic beads were incubated with 5 μg of antibody (anti-AGO1 or IgG) to form bead-antibody complex. The frozen lysates were thawed and centrifuged, and the supernatant was divided into three equal parts. Two parts of which were incubated with the bead-antibody complex (anti-AGO1 or IgG) at 4°C overnight, and the last one considered as "input" was stored at − 80°C until RNA extraction. Subsequently, the RNA was extracted by Trizol reagent (Life Technology) and reverse-transcribed into cDNA using the highcapacity RNA-to-cDNA Kit (Applied Biosystems). The precipitated mRNA of acox1 from both agomir-31-5pand agomir-control-treated lysates were quantified by qPCR. The IgG control and "input" were assayed to normalize the relative expression and ensure the specificity of the RNA-protein interactions. Six replicates were used for each treatment.
Heterologous expression of Mal-ACOX1 in C. elegans To validate the role that Mal-acox1 plays in the regulation of asc-C9, the amount of ascarosides were compared among C. elegans wild type strain N2 (Bristol), peroxisomal β-oxidation mutants acox1 (ok2257) I strain VC1785, and the rescue type strains (acox1 (ok2257) :: Cel-p:: Mal-acox1) with heterologous expression of Mal-ACOX1 in C. elegans. Eight replicates were used for each treatment.
To get the rescue strains, the promoter of Cel-acox1 (pCel) and CDS of Ma-acox1 were cloned into the pPD95_ 77 plasmid vector using the in-fusion cloning system (Additional file 3: Dataset S7, S8). The recombinant plasmids were injected into the gonads of mutant acox1 (ok2257) I strain VC1785 with the green fluorescent protein (GFP) plasmids. F2 and subsequent offspring with GFP were screened by a fluorescent microscope. Finally, two lines (line 1 and line 2) were selected for ascaroside analysis.

Ascaroside feeding experiment
To test the effect of ascaroside on beetle's SCPs and mortality, ascaroside feeding experiment was conducted as previously described [23]. Four hundred microliters of asc-C9, asc-△C6 (negative control), or sterilized water (control) was added into the artificial diets (200 g of sawdust and 15 g of agar were added to 550 ml of distilled water and mixed, and 400 μl of asc-C9 at 0.1 nM was added). As previously described, asc-C9 and asc-△C6 were synthesized and tested at 0.1 nM and 3.97 nM, respectively [23]. Then, larval beetles were reared on ascaroside added or control diets in a climate chamber at a 12:12 h light:dark (L:D) cycle at 25°C. Twenty individuals were sampled for SCP measurement from each treatment every day for 10 days. Samples from the 10th day were used for qPCR and mortality tests. Six and three replicates were used for each treatment in qPCR and mortality tests, respectively.

Statistical analyses
For mortality experiments, a probit regression was used to calculate LT 50 values and differences were compared by extra sum-of-squares F test. Data analyses were performed using GraphPad Prism (version 6.01) (GraphPad software, Inc., San Diego, CA, USA) [61]. For other experiments, Levene's test was used to test the assumption of homoscedasticity. Student's t test was used for two-group comparisons. One-way ANOVA with Tukey's HSD multiple comparison was performed for more than two-group comparisons. All statistical analyses were performed by IBM SPSS 18.0 software. Differences were considered significant at P < 0.05.