# Characterization and expression profiling of microRNAs in response to plant feeding in two host-plant strains of the lepidopteran pest Spodoptera frugiperda

Yves Moné, Sandra Nhim, Sylvie Gimenez, Fabrice Legeai, Imène Séninet, Hugues Parrinello, Nicolas Negre, Emmanuelle d'Alençon

► To cite this version:

Yves Moné, Sandra Nhim, Sylvie Gimenez, Fabrice Legeai, Imène Séninet, et al.. Characterization and expression profiling of microRNAs in response to plant feeding in two host-plant strains of the lepidopteran pest Spodoptera frugiperda. BMC Genomics, 2018, 19 (1), pp.1-15. 10.1186/s12864-018-5119-6. doi:10.1186/s12864-018-5119-6

HAL Id: hal-01926342

https://inria.hal.science/hal-01926342v1

Submitted on 18 Jul 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
RESEARCH ARTICLE

Open Access

# Characterization and expression profiling of microRNAs in response to plant feeding in two host-plant strains of the lepidopteran pest *Spodoptera frugiperda*

<figure><img src="image_2.png" /></figure>

CrossMark

Yves Moné<fnref n="1" />†, Sandra Nhim<fnref n="1" />†, Sylvie Gimenez<fnref n="1" />†, Fabrice Legeai²,³, Imène Seninet<fnref n="1" />, Hugues Parrinello⁴, Nicolas Nègre<fnref n="1" /> and Emmanuelle d'Alençon<fnref n="1" />*

**Abstract**

**Background:** A change in the environment may impair development or survival of living organisms leading them to adapt to the change. The resulting adaptation trait may reverse, or become fixed in the population leading to evolution of species. Deciphering the molecular basis of adaptive traits can thus give evolutionary clues. In phytophagous insects, a change in host-plant range can lead to emergence of new species. Among them, *Spodoptera frugiperda* is a major agricultural lepidopteran pest consisting of two host-plant strains having diverged 3 MA, based on mitochondrial markers. In this paper, we address the role of microRNAs, important gene expression regulators, in response to host-plant change and in adaptive evolution.

**Results:** Using small RNA sequencing, we characterized miRNA repertoires of the corn (C) and rice (R) strains of *S. frugiperda*, expressed during larval development on two different host-plants, corn and rice, in the frame of reciprocal transplant experiments. We provide evidence for 76 and 68 known miRNAs in C and R strains and 139 and 171 novel miRNAs. Based on read counts analysis, 34 of the microRNAs were differentially expressed in the C strain larvae fed on rice as compared to the C strain larvae fed on corn. Twenty one were differentially expressed on rice compared to corn in R strain. Nine were differentially expressed in the R strain compared to C strain when reared on corn. A similar ratio of microRNAs was differentially expressed between strains on rice. We could validate experimentally by QPCR, variation in expression of the most differentially expressed candidates. We used bioinformatics methods to determine the target mRNAs of known microRNAs. Comparison with the mRNA expression profile during similar reciprocal transplant experiment revealed potential mRNA targets of these host-plant regulated miRNAs.

**Conclusions:** In the current study, we performed the first systematic analysis of miRNAs in Lepidopteran pests feeding on host-plants. We identified a set of the differentially expressed miRNAs that respond to the plant diet, or differ constitutively between the two host plant strains. Among the latter, the ones that are also deregulated in response to host-plant are molecular candidates underlying a complex adaptive trait.

**Keywords:** Adaptation, Phenotypic plasticity, Regulation of gene expression, Insect plant interaction, microRNAs

<figure><img src="image_3.png" /></figure>

BMC

© The Author(s). 2018 **Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

<footnote>* Correspondence: emmanuelle.d-alencon@inra.fr</footnote>

<footnote>†Yves Moné, Sandra Nhim and Sylvie Gimenez contributed equally to this work.</footnote>

<footnote n="1">¹DGIM, Univ Montpellier, INRA, Montpellier, France</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
## Background

Successful adaptation to the host-plant is of fundamental importance to herbivorous insects. It requires that the adult female accepts the host-plant for oviposition and that the host allows feeding and proper development of larvae. Understanding of the underlying genetic mechanisms used by insects in response to their host-plants ([1] for review) recently progressed thanks to availability of new reference genomes and transcriptomes of major phyphagous herbivorous insect pests. Genomic sequences analyses highlighted expansion of chemosensory and/or detoxification genes in generalist herbivores compared to specialist ones reflecting their larger diets [2–6]. RNA-seq analyses revealed that generalist herbivores use transcriptional plasticity of various categories of genes in response to their diet [3, 6–9]: detoxification, digestion, cuticular and ribosomal genes. Transcriptional plasticity of specific sets of genes has also been shown in an oliphagous lepidopteran species, *Manduca sexta*, when it is reared on host -or acceptable non-host plants [10]. These data show that a change of host-plant in herbivorous insect requires large scale transcriptional changes involving combinations of various gene family members.

While the nature and role of protein coding genes involved in adaptation to the host-plant begin to emerge, the putative role of non-coding genes has yet to be explored. Among them, microRNAs (miRNAs) are a class of small non coding RNAs (sncRNAs) of approximately 22 nt in length, which act as post-transcriptional regulators of gene expression and are known to help fine-tune complex genetic networks ([11], for review). The mode of action of miRNAs results in relatively weak modulation of less than twofold both at the RNA and protein levels [12]. Two other classes of small non coding RNAs combat the invasion and the expansion of transposable elements (TE), the short interfering RNA (siRNA) pathway suppresses TEs in all tissues of plant and animals, whereas the Piwi-interacting RNA pathway (piRNA) seems more specific to the gonads of metazoans (see [13] for review).

Since miRNAs have been shown to be involved in many physiological or cellular processes such as differentiation, proliferation, apoptosis and development [14], we hypothesized that they may play a role in adaptation of a phytophagous insect to its host-plants and that they may show different expression patterns in different host-plant races of the insect. To test this hypothesis, we used the noctuid moth *Spodoptera frugiperda*, which consists of two host-plant strains, one mostly associated to corn (C strain or SfC), the other to rice (R strain or SfR), and whose genomes are recently available [4]. We performed reciprocal transplant (RT) experiments of the two strains on the two host-plants and isolated and sequenced sncRNAs from feeding larvae. We present the differential expression patterns of miRNAs and their

putative coding genes targets involved i) in phenotypic plasticity (the ability of a single genotype to produce multiple phenotypes in response to variation in the environment) of each strain in response to corn or rice ii) in adaptive evolution or genetic drift, by additional comparison of the two strains on the same host-plant.

## Results

### Deep sequencing of *S. frugiperda* small RNA

To characterize *S. frugiperda* miRNA, small RNA libraries were constructed from whole body of corn-fed and rice-fed larvae of *S. frugiperda* C and R strains. Two independent libraries (biological replicates) were prepared and independently sequenced using Illumina technology. Between 31.5 to 57 millions high quality sequence reads were obtained after adapter trimming in each library (See Additional file 1: Table S1). Their size distribution shows an over-abundance of sequences at 22 nt, typical of miRNAs (Fig. 1a-b). On average, we detected 11.12 +/- 3.67% of sequences of size 22 nt on corn and 5.49 +/- 1.06% on rice. We detected a second peak between 25 and 33 nt that could correspond to expression of piRNAs. After collapsing the reads, we identified between 2 to 5 million unique sequences, 4% of which (4.2+/- 0.2% on corn and 3.71+/- 0.52% on rice) on average corresponding to miRNAs of 22 nt (Fig. 1c, d). However, in term of diversity, piRNAs seem more abundant than miRNAs (Fig. 1c, d) which may reflect abundance and diversity of transposable elements (TE) in the genome of *S. frugiperda* (29.14% and 29.10% of genome coverage by TE in C and R strain, respectively [4]).

In eukaryotic cells, the vast majority of cellular RNA consists of ribosomal RNA (rRNA) ~ 80 to 90% of total RNA for most cells, followed by transfer RNA (tRNA) - 10 to 15%, messenger RNA (mRNA) - 3 to 7%, miRNA - 0.003 to 0.02% [15]. The small ncRNA sequences of SfC or SfR were aligned against these different references, in addition to the sequences of TE copies, which should highlight putative piRNAs or endo-siRNAs (Fig. 1e, f), The most abundant hits (29–29.1%) were found matching to Sf miRNAs precursors (precursors of known or novel miRNAs sequences) or Sf TE copies (15.5–18.4%) corresponding to putative piRNAs or endo-siRNAs. 16–16.7% of them matched to rRNA, and 0.2–2.4% to tRNA, the most abundant RNA classes expected in the cells. These data show that the sRNA sequences were enriched in functional miRNAs or TE-interacting RNAs compared to those resulting from degradation of ribosomal or transfer RNAs. We found also sncRNA sequences matching to miRNA precursors of plants (*Zea* or *Oryza*). It is expected since sncRNAs have been extracted from whole larvae feeding on plants. The lack of sequence homology between miRNA families in plants and animals (with exception of one family [16]) as well
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

**Fig. 1** Size profiling of small non coding RNAs and their homology to different RNA classes or to Transposable Elements (TE). The percentage of sncRNA reads is plotted as a function of their size (between 15 nt to 40 nt corresponding to the size range that has been selected from the gel for library construction), a and c SfC, c and d SfR, in green on corn, in red on rice. CC: SfC on corn, CR: SfC on rice, RC: SfR on rice, RR: SfR on rice. a and b total reads, c and d unique reads. e and f Pie charts representing the average % of reads (total counts from 2 replicates on corn for SfC (e) or SfR (f)) mapping either to SfC or plant miR precursors, or TE (SfC TE copies) as expected for putative endo-siRNA or piRNAs, or mRNA (SfC OGS2.2), or SfC tRNA, SfC rRNA (18S and 28S RNAs)

as differences in their biogenesis and mode of action suggested that miRNAs have evolved independently in both kingdoms from an ancient siRNA mechanism already present in the last common ancestor of all eukaryotes [12]. The presence of plant miRNA reads in larval samples is thus not expected to interfere with the study of animal miRNAs performed in this paper.

## Annotation of *S. frugiperda* miRNA genes
miRDeep2 analysis

To detect miRNA genes, raw sequencing data were analyzed with miRDeep2 software. miRDeep2 [17, 18] maps the sRNA reads to the genome and excises potential miRNA precursors sequences from the genome. The secondary structures of the miRNA precursors are

predicted and their stability is estimated by RNAfold. mirDeep2 uses a probabilistic model of miRNA biogenesis by the Dicer protein to score frequency and compatibility of mapping of the small RNA sequence reads (the signature) on the secondary structure of the miRNA genomic precursors (the structure) as compared to a non-miRNA precursor hairpin [17, 18]. Read stacks correspond to mature miRNA sequences. The score reflects the likelihood of each precursor to be a genuine miRNA. Furthermore, since the algorithm may generate hairpins with read stacks that have no connection with miRNA biology, corresponding to false positives, miRDeep2 estimates the rate of false positive by shuffling the observed combinations of structures and signatures and submit them to the core algorithm. The difference in score
distribution between the genuine combinations and the control ones is used to estimate the number of false positive novel miRNAs for varying score cut-offs. The sequence of mature predicted miRNAs are compared to mature miRNA sequences contained in miRBase (release 21) which allows to sort them in two classes, known or novel depending if they existed or not in miRBase. The genomes of the C strain (v3.1) and the R strain (v1.0) were used as references using a pool of sequence reads from either the C or the R strain, respectively. As shown on Table 1, we obtained 76 and 68 known miRNAs predicted genes (70 and 64 unique ones) and 139 and 171 novel ones (126 and 158 unique ones) in the C and R strain.

We calculated for the 139 novel miRs predicted in SfC, the processing precision frequency, defined as the ratio of reads corresponding exactly to the mature miRNA and miRNA\* sequences, divided by the total amount of reads mapping to the hairpin [19] (See Additional file 2: Excel file S1, tabs “SfC Novel”). A value close to one indicates high precision and a value close to zero indicates the production of very few miRNA/ miRNA\* duplexes, with a cut-off of less than 0.1 corresponding to low processing precision [20]. One hundred thirty six novel miRs (97.8%) have an efficiency of more than 0.1, and 111 (79.8%) have an efficiency of > 0.5 thereby showing a medium to high processing precision in *S. frugiperda* and\*or that the mirDEEP2 prediction generated few false positives.

After filtering for the prediction showing a score > 4 (the lowest score cut-off corresponding to a prediction signal-to-noise ratio *r* > 10, *r* = total miRNA hairpins reported/ mean estimated false positive miRNA hairpin over 100 rounds of permutated control) and a significant randfold *p*-value, we obtained 66 and 59 genes for known miRNAs and 78 and 102 for novel ones in C and in R strain. All miRNA genes prediction can be found in Supplementary excel file 1 (See Additional file 2). Compared to the miRNAs identified from ovary cell lines of *S. frugiperda* (Sf21) [21] whose sequences are not in the miRBase release 21, we found 40 to 36 additional predictions of known miRNAs genes and 103 to 155 novel ones in C and R strain.

Other miRs involved in response to Baculovirus infection have been described from another ovary cell lines of *S. frugiperda*, Sf9, with precursor sequences mapped on *Bombyx mori* genome [22].

The specific nucleotide occurrence was analyzed in the obtained miRNA sequences, *S. frugiperda* showed a dominant bias for uracil (U) at the first nucleotide (Fig. 2). The dominance of U at first position towards 5'end is a conserved feature of miRNAs [23].

## Orthology

To identify orthologous miRNAs between C and R strain a reciprocal blastn was performed between mature sequences of known or novel miRNAs at an e-value < 0.001. We identified 57 orthologs among known miRNAs and 75 among novel ones. The list of orthologs can be found in Additional file 2: Excel file S1. This step allowed comparison of expression levels of miRNAs according to the genetic background (C or R strain).

## miRNA expression variation

### DESeq2 analysis according to host-plant

Expression values generated by miRDeep2 in the RT experiments were used to analyze differential expression of *S. frugiperda* miRNAs according to the host-plants in each strain. In Additional file 3: Figure S1, for SfC in the upper part and SfR in the lower part, the MA-plots show the log2 fold change on rice versus corn over the mean of normalized counts, i.e. the average of counts normalized by size factors. The points with FDR less than 0.05 are colored in red. We also visualized samples (treatments, replicates) by Principal Component Analysis (PCA) as shown in Additional file 4: Figure S2. Most variation was linked to treatment, i.e., change of host-plant, the smallest variation was found between biological replicates.

Based on read counts, 34 known or novel miRNAs (out of 144, 23.6%) were differentially expressed (FDR < 0.05) in the C strain larvae fed on rice as compared to the C strain larvae fed on corn (Details of DESeq2 analyses can be found in Additional file 5: Excel file S2.Twenty one known or novel miRNAs (out of 161, 13%) were differentially expressed on rice compared to corn in R strain. For example, in the C strain, known miR-34-5p and miR-190-5p were overexpressed on rice (in red on Fig. 3a), while novel miR-375-5p was overexpressed on corn (in green). In the R strain (Fig. 3b), novel tca-miR-375-5p and known miR-190-5p were overexpressed on rice as compared to corn (in red). Novel mmu-miR-155-3p was overexpressed on corn compared to rice (in green). We used TaqMan RT-qPCR miRNA assays on total RNA extracted from the

<table><caption>**Table 1** Number of miRNAs genes predicted in the *Spodoptera frugiperda* genome</caption><thead><tr><th>Species</th><th>Known</th><th>New compared to Kakumani<sup>a</sup></th><th>Known and unique</th><th>Known Score > 4 (+randfold yes)</th><th>Novel</th><th>New compared to Kakumani<sup>a</sup></th><th>Novel And unique</th><th>Novel Score > 4 (+randfold yes)</th></tr></thead><tbody><tr><td>SfC</td><td>76</td><td>40</td><td>68</td><td>66 (66)</td><td>139</td><td>103</td><td>126</td><td>92 (78)</td></tr><tr><td>SfR</td><td>68</td><td>36</td><td>64</td><td>61 (59)</td><td>171</td><td>155</td><td>158</td><td>115 (102)</td></tr></tbody></table>

<footnote><sup>a</sup>Since the *S. frugiperda* microRNAs predicted by Kakumani et al., [21] are not registered in miRBase, we checked whether each of our predicted miRs had been predicted in Sf21 genome</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 2** Base composition of known or novel mature miRNAs</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 3** Differential expression of miRNAs genes on rice compared to corn in *Spodoptera frugiperda* larvae (L4 instar), after rearing for 3 generations on whole plants. The miRNAs showing a significant differential expression after DESeq2 analysis (log2foldchange > 1 or < 1 and FDR < 0.05) are shown. a In SfC b In SfR. In red, miRNAs up-regulated on rice, in green miRNAs up-regulated on corn</figcaption></figure>

RT experiments to validate expression variation based on read counts, for some of the most differentially expressed miRNAs (Table 2). We could confirm that known miR-34-5p and miR-190-5p were significantly overexpressed on rice compared to corn in the C strain and that known miR-190-5p and novel tca-miR-375-5p were significantly overexpressed on rice as compared to corn in the R strain (Table 2).

### miRNA expression variation between strains

We compared expression values between the two strains when they were reared on the same host plant to detect differences linked to the genetic background. In Additional file 6: Figure S3, on corn (top panel) and on rice (bottom panel), the MA-plots show the log2 fold change in SfR versus SfC over the mean of normalized counts. Again by PCA analysis of samples, more variation was found between strains than between biological replicates (See Additional file 7: Figure S4).

Based on reads counts, nine known or novel miRNAs (out of 129, 6.97%) were differentially expressed (FDR < 0.05) in the R strain compared to C strain when reared on corn (See Additional file 8: Excel file S3). A similar number of miRNAs (10 out of 129, 7.75%) were differentially expressed between strains on rice.

For example, known miR-34-5p and novel dme-miR-275-3p were overexpressed in SfR compared to SfC on corn (in red on top of Fig. 4), the latter being also overexpressed in SfR compared to SfC on rice (in red on bottom of Fig. 4). Using TaqMan RT-qPCR miRNA assays, we could confirm upregulation of miR-34-5p in SfR
<table><caption>**Table 2** Experimental validation of variation in miRNA expression</caption><thead><tr><th>Experimental condition</th><th>MicroRNA</th><th>Relative expression</th><th>Std Error</th><th>95% C.I.</th><th>P(H1)</th><th>Result</th></tr></thead><tbody><tr><td rowspan="3">SfC (Rice/Corn)</td><td>miR-31 (Ref)</td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>miR-34</td><td>1.763</td><td>(1.080-3.062)</td><td>1.003-3.513</td><td>0.004</td><td>UP</td></tr><tr><td>miR-190</td><td>1.549</td><td>(1.320-1825)</td><td>1.235-2.032</td><td>0.001</td><td>UP</td></tr><tr><td rowspan="3">SfR (Rice/Corn)</td><td>miR-31 (Ref)</td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>miR-190</td><td>1.218</td><td>0.671-1.520</td><td>0.900-1.669</td><td>0.05</td><td>UP</td></tr><tr><td>tca-miR375*</td><td>5.802</td><td>(3.115-11,105)</td><td>1.511-14.148</td><td>0.002</td><td>UP</td></tr><tr><td rowspan="3">On Corn (SfR/SfC)</td><td>miR-31(Ref)</td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>miR-34</td><td>1.868</td><td>1.395-2.466</td><td>1.282-2.811</td><td>0.003</td><td>UP</td></tr><tr><td>dme-miR-275*</td><td>0.965</td><td>0.821-1.121</td><td>0.726-1.251</td><td>0.587</td><td></td></tr><tr><td rowspan="2">On Rice (SfR/SfC)</td><td>miR-31(Ref)</td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>dme-miR-275*</td><td>1.434</td><td>1.236-1.545</td><td>1.194-1.907</td><td>0.001</td><td>UP</td></tr></tbody></table>

The relative expression of miR genes depending on the host-plant (on rice compared to corn) in either the C or the R strain, or depending on the genetic background (SfR compared to SfC) was calculated according to [57]

P(H1) Probability of alternate hypothesis that difference between sample and control groups is due only to chance. The miRs with a star correspond to * Novel miRs having conserved seed region only

compared to SfC on corn, and that of the novel dme-- miR-275-3p in SfR compared to SfC on rice (Table 2), however not on corn, maybe due to the presence of a RNA molecule competing with the microRNA for the TaqMan probe in this condition.

<figure><img src="image_14.png" /></figure>

## Expression differences both between strains and between plants

The expression differences between host-strains may result from genetic drift and also possibly from divergent selection by the environment, in this case the different host-plants. To identify the ones putatively involved in adaptation to the host-plant, we looked for the miRNA genes of the known class that were differentially expressed both between strains and between plants (FDR < 0.05). As shown on Fig. 5a, four miRs (miR-10-5p, miR-34-5p, miR-263a-5p, miR278-5p) were differentially expressed both constitutively between strains on corn, and within SfC when reared on different plants. The two miRNAs genes that were differentially expressed between strains on rice (miR-279b-3p and miR-308-3p) were also differentially expressed in SfC on different plants. On Fig. 5b, we found that among the 4 miRNAs that were differentially expressed between SfR and SfC on corn, only one, miR-278-5p, was also differentially expressed when SfR was reared on rice compared to corn. None of the two miRNAs that were differentially expressed between strains on rice was differentially expressed in SfR according to the host-plant. Most constitutive differences between strains are involved in interaction with the plant in SfC, less of them in SfR, suggesting that adaptation to the plant played a more pronounced role in SfC evolution than in SfR.

## Potential target genes regulated by miRNA

mRNA targets of known miRNAs were searched using TargetScan, miRanda, Rna22 and miRmap (See Methods) against the 3' UTR of the C strain gene set, OGS2.2. The target gene list of the seven DE miRNAs identified in this study (miR-34, miR-190, miR-1a-5p, miR-998, miR-278, miR-263a-5p, miR-10-5p) can be found in Additional file 9: Excel file S4. Among the gene targets predicted by
<figure><img src="image_15.png" /></figure>

TargetScan, we filtered out those showing variation in their expression in the RT experiments (according to [7]). We also checked whether they were predicted as miR targets by at least one of the three other softwares. The list of target genes that are differentially expressed in SfC on different host plants can be found in Additional file 10: Excel file S5. The target genes that are differentially expressed in SfR compared to SfC on corn are listed in Additional file 11: Excel file S6.

## Potential regulated target genes

To reduce the number of putative false positive coding gene targets among those predicted by TargetScan, we assumed that true targets should be expressed in the same experimental conditions as miRNA genes, and downregulated when miRNA genes were overexpressed. To identify these candidate genes, we used RNA-Seq results obtained in the same reciprocal transplant experiments after two generations on plants [7]. We limited the analysis of target genes to the Official Gene Set of the C strain, which is the reference and has been manually curated [4], and to RT experimental conditions in which the RNA-Seq data of putative target genes were available in duplicates in [7]. We thus provide a detailed analysis of targets genes in the following conditions: 1) SfC reared on corn and on rice 2) SfC and SfR when reared on corn. We focused on the most differentially expressed miRNAs of “known” class (FDR < 0.05 and absolute value of log2 fold change > 1.3) in these conditions. The complete list of differentially expressed (FDR < 0.05) coding gene targets of these known miRNA genes, with their log2 fold change and annotation can be found in Additional file 10: Excel file S5 and Additional file 11: Excel file S6. miR-34 and miR-190 are up-regulated in SfC when reared on rice, a non-host plant. Among down-regulated targets of miR-34, we found members of different gene families (Table 3), some of which are also targets of miR-190: First, a representative of the takeout gene family, then three genes encoding cuticle proteins. A gene encoding a subunit of the 26S proteasome was specifically targeted by miR34. As miR-190 specific targets, were found a cuticle protein gene, an acyl-CoA desaturase gene, and a member of the Osiris gene family, osi9a [7].

Among miRNAs that are overexpressed in SfC on corn compared to rice, three, miR-998, miR-263a-5p and miR-10-5p share two targets that encode transporters (facilitated trehalose transporter Tret1 and a sodium- and of chloride-dependent glycine transporter 1) that are downregulated on corn. In addition, a gene encoding a spermine oxidase enzyme is the target of both miR-998 and miR-263a-5p.

When we compared the two strains on the same host-plant corn (See Additional file 11: Excel file S6), we found that miR-34 was up-regulated in SfR and that miR-263a-5p, miR10-5p and miR-278-5p were up-regulated in SfC. Among the down-regulated targets of miR-34, we found a representative of the transcription factor family AP-2 (Table 3). Among the down-regulated targets of miR-263a-5p, were found a protein of unknown function and a cuticular protein. Among the targets of miR10-5p, two genes coding for transporters (proton-coupled folate transporter, facilitated trehalose transporter Tret1) were found.
<table><caption>**Table 3** List of target genes of known miRs and their differential expression</caption><thead><tr><th rowspan="2">MiR<br>Name</th><th rowspan="2">Sf strain</th><th rowspan="2">Regulated on</th><th colspan="4">Putative target genes</th><th rowspan="2">Log2FC</th></tr><tr><th>Log2FC</th><th><sup>a</sup>Predicted Targets</th><th>Down or up regulated targets</th><th>Down-regulated</th></tr></thead><tbody><tr><td>miR-34</td><td>Sf corn</td><td></td><td>2.32</td><td>321</td><td>38</td><td>28</td><td>take-out<br>GSSPFG00021718001-RA</td><td>-4.22</td></tr><tr><td>miR-190</td><td></td><td></td><td>2.07</td><td>479</td><td>46</td><td>25</td><td>cuticular proteins<br>GSSPFG00000626001-RA<br>GSSPFG00013845001-RA<br>GSSPFG00010502001-RA</td><td>-3.24<br>-1.6<br>-1.21</td></tr><tr><td>miR-1a-5p</td><td></td><td></td><td>1.33</td><td>524</td><td>46</td><td>18</td><td>subunit DSS1 of 26S proteasome<br>GSSPFG00027010001-RA</td><td>-2.14</td></tr><tr><td>miR-998</td><td></td><td></td><td>-3.33</td><td>229</td><td>31</td><td>10</td><td>cuticular protein<br>GSSPFG00005155001-RA<br>acyl-CoA desaturase<br>GSSPFG00006314001-RA</td><td>-3.54<br>-6.54</td></tr><tr><td>miR-278</td><td></td><td></td><td>-2.50</td><td>91</td><td>8</td><td>4</td><td>osiris 9A<br>GSSPFG00012033001-RA</td><td>-4.71</td></tr><tr><td>miR-263a-5p</td><td></td><td></td><td>-2.18</td><td>436</td><td>56</td><td>27</td><td>Serine proteinase<br>GSSPFG00014426001-RA</td><td>-3.38</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>sodium- and chloride-dependent glycine transporter 1-like<br>GSSPFG00018036001-RA</td><td>2.61</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>facilitated trehalose transporter Tret1<br>GSSPFG00023681001-RA</td><td>1.31</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>Spermine oxidase<br>GSSPFG00025432001-RA</td><td>3.02</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>integrin alpha-PS2<br>GSSPFG00025987001</td><td>1.15</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>sodium- and chloride-dependent glycine transporter 1-like<br>GSSPFG00018036001-RA</td><td>2.61</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>facilitated trehalose transporter Tret1<br>GSSPFG00023681001-RA</td><td>1.31</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>Spermine oxidase<br>GSSPFG00025432001-RA</td><td>3.02</td></tr></tbody></table>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>
<figure><img src="image_19.png" /></figure>

of different gene families, some of which are also targets of miR-190: First, a representative of the takeout gene family, putatively involved in feeding behavior and response to starvation as in *D. melanogaster* [26]. This representative is targeted both by miR-34 and miR-190. Second, three genes encoding cuticle proteins whose downregulation may reflect slower development of SfC

on rice compared to corn [7]. Another cuticle component is also targeted by miR-190. Among miR-34 specific targets, we found a gene encoding a subunit of the 26S proteasome, which can be involved in protein degradation in response to oxidative stress. In the case of phytophagous insects, oxidative stress can be generated by prooxidant allelochemicals produced by host-plants.
Among miR-190 specific targets, an acyl-CoA desaturase gene, necessary for fatty acid biosynthesis, may be repressed because rice is a poor food for SfC, and a member (*Osiris* 9) of the Osiris gene family, putatively involved in response to plant toxins as in *Drosophila sechellia* [27]. A recent analysis of the conserved patterns of Osiris gene expression in different insect species, suggests that Osiris genes may play a central role in insect adaptive evolution [28].

The three miRNAs miR-998, miR-263a-5p and miR-10-5p that are overexpressed in SfC on corn share two targets that encode transporters (facilitated trehalose transporter Tret1 and a sodium- and of chloride-dependent glycine transporter 1) that are downregulated on corn (up-regulated on rice). In most insects, trehalose (a-D-glucopyranosyl-(1,1)-a-D-glucopyranoside) is the main haemolymph sugar. In Drosophila, Tret 1 is necessary for the transport of trehalose produced in the fat body and its uptake into other tissues that require a carbon source, and thereby regulates trehalose levels in the hemolymph [29]. The glycine transporter GLYT1, by controlling the re-uptake of glycine at synapses [30], regulates neurotransmission, where glycine plays the role of inhibitory neurotransmitter. In addition a gene encoding a spermine oxidase enzyme is the target of both miR-998 and miR-263a-5p. Polyamines (PA), comprising spermine (Spm), spermidine (Spd) and putrescine (Put), are ubiquitous polybasic molecules, with many important biological functions, like cell growth, differentiation, and apoptosis. They interact reversibly with nucleic acids, regulating chromatin status and gene expression, and modulating ion-channels’ function and stability (reviewed in [31]). Polyamine oxidases (PAO) include spermine oxidase. These enzymes, containing a flavin adenine dinucleotide (FAD), catalyze the oxidation of polyamines, and lead to formation of hydrogen peroxide. They regulate cellular polyamine concentration.

When we compared the two strains on the same host-plant corn (See Additional file 11: Excel file S6), we found that miR-34 was up-regulated in SfR and that miR-263a-5p, miR10-5p and miR-278-5p were up-regulated in SfC. Among the down-regulated targets of miR-34, we found a representative of the transcription factor family AP-2. From [32], AP-2 is expressed during Drosophila embryogenesis in the maxillary segment and neural structures, whereas during larval development, it is expressed in the central nervous system (CNS) and the leg, antennal and labial imaginal disks [33, 34]. In a Drosophila AP-2 mutant, defects in proboscis development and leg-joint formation have been described [35, 36]. We found a homolog of the graves disease carrier protein, a protein of as yet uncharacterized function that belongs to the mitochondrial metabolite carrier family (which includes the ADP/ATP translocator, the phosphate carrier

and the hydrogen ion uncoupling protein). Among the down-regulated targets of miR-263a-5p, were found a protein of unknown function and a cuticular protein. Among the down-regulated targets of miR-10-5p, two genes coding for transporters (proton-coupled folate transporter, facilitated trehalose transporter Tret1) were identified. Folates are a family of B9 vitamins found in nature primarily as 5-methyltetrahydrofolate (5-methylTHF). 5-MethylTHF provides the methyl group for the synthesis of methionine from homocysteine and therefore is necessary for the formation of S-adenosylmethionine, which is required for a variety of methylation reactions [37]. In addition, folates are the sources of a methylene moiety for the de novo synthesis of thymidylate from deoxyuridylate and two formate moieties for the de novo synthesis of the purine ring. PCFT has been extensively studied in humans due to the importance of folates in cancer progression, however it has not been studied in insects. In humans, PCFT is expressed at the acidic microenvironment of the apical brush-border membrane of the proximal small intestine and allows the intestinal absorption of folates. In humans and mice, loss of function mutations in PCFT, results in severe systemic folate deficiency with anemia, sometimes pancytopenia, hypo-immunoglobulinemia and gastrointestinal defects [38].

This analysis of putative targets of DE miRNA reveals that they control important pathways necessary for cellular or organismal homeostasis during insect-plant interaction. Although we provide transcriptomic evidence of downregulation of putative targets, their experimental validation will require additional efforts like the use of biotin tagged miRNAs to capture them and/or translation profiling [39, 40], which will be the focus of future work.

For expression analysis of the predicted target genes, insect samples were collected after two generations on plants whereas for the study of miRNA expression they were collected after three generations. This was to ensure dilution of miRNAs putatively maternally inherited (it is the case for miR-34 in *D. melanogaster* for instance [41]) and related to the artificial diet on which the insect fed before the reciprocal transplant experiment on plants started. We considered that in controlled experimental conditions of growth, gene expression during each successive insect developmental cycle is reproducible and comparable from one generation to another (Since *S. frugiperda* is a quarantine organism in France, we performed the RT experiment in large incubators with controlled hygrometry, temperature and light conditions). If the plant exerted a selection pressure, gene allele frequencies may have changed but not significantly between two generations since the laboratory strain that we used has a limited genetic polymorphism.

For differential expression analysis of miR genes and their targets, two biological replicates of the reciprocal transplant experiment have been performed. Since the first
use of RNA-Seq to analyze transcriptomic data [42], we
know more on the parameters necessary to optimize detec-
tion of differentially expressed genes. Lamarre et al. re-
cently showed that the read depth in part compensates for
the number of replicates to increase the ratio of differen-
tially expressed genes detected [43]: with 20 million reads
(for 20,000 DE genes - corresponding to a coverage of
1000 reads per gene) and 2 replicates, 85% of DE genes are
found, and 8 replicates are needed to reach a ratio of 100%,
while with 2.5 million reads (~ coverage 125 reads/gene)
and 2 replicates, only 15% of DE genes are detected and
the use of 8 replicates increases the ratio of DE genes to
only 60%. In our miR differential expression analysis, we
performed only two biological replicates. However, the read
depth in each SfC library (Additional file 1: Table S1) was
42 million on average, with more than 20% of the se-
quences having the expected size for miRs (21 to 23 nt),
this corresponds to an average coverage of > 39,688 reads
per miR gene, a higher coverage than the maximal
one used in the study of Lamarre et al.. Lamarre et
al. also showed that the rate of false positives ob-
tained with DESeq2 is minimal with 2 replicates and
increases with the number or replicates. According to
them, 70% of true positives can be detected with two
replicates and this number increases with the number
of replicates. We conclude that our experimental de-
sign enabled detection of a reasonable number of reli-
able DE miR candidates although more repetitions
may be necessary to deepen the study.

In the same line, our differential expression analysis
was done from RNA extracted from whole body of
the larvae feeding on plants. We are aware that this
experimental design may underestimate the number
of miRs showing variation in their expression due to
the fact that overexpression in one tissue may be
masked by down regulation in another one. It will be
the subject of future effort to look for expression
variation tissue by tissue.

The miRNAs expression differences that we uncovered
between the host-strains may result from genetic drift
and also possibly from divergent selection by the envir-
onment, in this case the different host-plants. To identify
the latter, we searched for the miRNAs that were differ-
entially expressed both between strains on the same
plant and within strains in response to the host-plant.
Interestingly, we found that all the miRNAs expression
difference between strains on plants (6/6, miR-10-5p,
miR-34-5p, miR-263a-5p, miR-278-5p, miR279b-3p,
miR-308-3p) were also involved in response to the
host-plant in SfC. By contrast, one miRNA only
(miR-278-5p) out of the six that showed constitutive dif-
ference between strains, showed also variation in re-
sponse to the host-plant in SfR. By measuring the fitness
of the insects on plants, we had found better survival of

the C strain on corn compared to rice [7] suggesting
that SfC was adapted to corn. The miRNAs that are both
deregulated between strains and in response to plant
may be involved in this adaptive evolution of SfC.

MiRNAs have been shown to be involved in response to
various environmental changes, like response to starvation
in *C. elegans* [44], to freezing and anoxia stress in the
freeze tolerant fly *Eurosta solidaginis* [45], to thermal plas-
ticity of the Senegalese sole [46], to drought in Tobacco
[47]. In the case of *Spodoptera frugiperda*, by regulating
phenotypic plasticity, miRNAs may have also played a role
in evolutionary adaptation as has been discussed in the
case of human miRNAs [48]. Further work is needed to
show that the miR candidates identified in this study are
directly involved in fitness of the insect on its host-plant,
however the possibility of transgenerational inheritance of
these molecules in male or female gametes [41, 49, 50]
suggests that they could facilitate transmission of complex
adaptive traits from parents to offspring.

Conclusions

In the current study, we performed the first systematic
analysis of miRNAs in Lepidopteran pests reared on
whole host and non-host plants. We identified a set of the
differentially expressed miRNAs that respond to the plant
diet, or differ constitutively between the two host plant
strains. The analysis of the putative targets of these DE
miRNAs revealed that they control important pathways
necessary for cellular or organismal homeostasis during
insect-plant interaction. Since two classes of non coding
RNAs, siRNAs and miRNAs have been used to control in-
sect development [51, 52], this study of regulatory mole-
cules of insect - plant interaction can bring clues on novel
environment friendly biological control of crop pests.

Methods

Biological material

Two laboratory strains of *S. frugiperda* were used in this
study: the corn-strain, originated from French
Guadeloupe and the rice-strain, originated from Florida,
USA. Each strain was reared on its principal and alterna-
tive host plant (*Zea mays* L. cv B73 or *Oryza sativa* L.
*Tajikica* cv Arelate) under controlled conditions
(temperature: 24 °C, photoperiod 16:8 light:dark, relative
humidity: 65%). Eggs from *S. frugiperda* were deposited
on plant leaves and fed ad libitum. After three generations
on plants, larvae of fourth instars (L4) were collected.
Therefore, larvae selected for small RNA extraction and
sequencing comprised four groups: C strain reared on
corn (CC), R strain reared on corn (RC), C strain reared
on rice (CR) and R strain reared on rice (RR). Each experi-
ment comprised two independent replicates.
## Small RNA extraction and sequencing

Small RNA was extracted from 12 L4 larvae using the mirVana miRNA Isolation Kit (Ambion) according to the manufacturer's instructions. Extracted small RNA was quantified using a NanoDrop 2000 spectrophotometer (Thermo Scientific) and analyzed by PAGE and silver-staining. Small RNA samples from each replicate of CC, RC, CR and RR larvae were used for cDNA library preparation using the TruSeq *Small RNA Sample Prep Kit (Illumina) and sequencing on Illumina HiSeq 2500. A size selection targeting snRNAs in the range of 15 to 40 nt was done. Library preparation and sequencing were performed by the platform MGX-Montpellier GenomiX (Montpellier, France).

Computational analysis of small RNA sequencing data and miRNA identification.

High throughput sequencing generates small RNA reads of 50 nucleotides in length (single reads). Raw reads were trimmed to remove adapter sequences using cutadapt software (version 1.4.1) [53]. For downstream analysis, only reads having length between 15 and 40 bases were considered.

## miRNA genes annotation

miRDeep2 algorithm [18] was used to detect miRNA from small RNA deep sequencing data. This algorithm uses a probabilistic model to score the fit of sequenced RNA to the biological model of miRNA biogenesis [17]. Briefly, reads are aligned to the *S. frugiperda* genome (Corn variant assembly version 3.1 or Rice variant assembly version 1.0), and only reads that do not map more than five times to the genome were used for miRNA detection. Then, using the read mappings as guidelines, potential miRNA precursors are excised from the genome and the miRDeep2 core algorithm scores their likelihood to be a real miRNA precursor. The output is a scored list of known and novel miRNA in the deep sequenced sample. Known miRNAs were identified by similarity to miRNA sequences from miRBase database (release 21).

## Orthology analysis

To explore the conservation of miRNAs between C and R strains, a reciprocal blastn was performed to search for orthologs [54, 55]. BLAST for miRNA mature sequences was run with blastn with default parameters except for a lower e-value threshold of 1e-3.

## Analysis of differential miRNA expression

miRDeep2 algorithm also provides read counts for the detected miRNA. To assess changes in miRNA expression between *S. frugiperda* variants and host-plant conditions, the read counts data for known and novel miRNA were used as input for the R package DESeq2

[56]. DESeq2 uses negative binomial generalized linear models to test for differential expression. An adjusted *p*-value for multiple testing was computed with the Benjamini-Hochberg procedure to control false discovery rate (FDR). Results with a FDR < 0.05 were considered statistically significant.

## Experimental validation of differential expression

The miRNA expression levels were quantified using TaqMan small RNA assay system from Life Technologies. Briefly, total RNA from samples was isolated using Trizol. One to 10 ng of total RNA was used for reverse transcription using a specific RT primer with the following conditions, i.e., 16 °C: 30 min, 42 °C: 30 min, 85 °C: 5 min. Subsequently, the cDNA was used for qRT analysis with TaqMan probes according to the manufacturer's instructions. For all qRT based TaqMan assays, the qRT-PCR quantifications were performed on two biological replicates (pools of 12 larvae) of the reciprocal transplant experiments, with 3 technical replicates. During the qRT analysis, the 2-ΔΔCT method was employed and each Ct value of the test miR was normalized to that of an endogenous miRNA (*Sf mir-31 5p*) whose expression remained stable in the different experimental conditions that were tested in C and R strain. To check that the miR expression level differed significantly between two experimental conditions, we used a pair wise fixed reallocation randomization statistical test [57] (2000 iterations, *p*-value< 0.001) which avoids making any assumptions about distributions compared to standard parametric tests such as analysis of variance or t-tests.

## Detection and functional annotation of potential target genes regulated by miRNA

TargetScan [58–60] predicts biological targets of miRNAs by searching for the presence of conserved 8-mer and 7-mer sites that match seed region of each miRNA by calculating thermodynamic free energy using the RNAFold package [61]. Predictions are ranked using the site number, site type, and site context. TargetScan (version 5.0) was run with default parameters using mature known miRNAs and the 3'UTR of predicted coding genes set OGS2.2_UTR3 of the C strain available on the webportal [62]. Predictions obtained by TargetScan were confirmed by at least one of the three following softwares: MiRanda, Rna22 and miRmap. MiRanda (version v3.3a) [63] allows one wobble pairing in the seed region when it is compensated by matches in the 3' end of the miRNA, it calculates the binding energy of the duplex structure and its position within the 3'UTR, it was used with the same parameters as in [64]. Rna22 (version v2) [65] is a tool based on a search for patterns that are statistically significant miRNA motifs created after a sequence analysis of known mature miRNAs, it was used
with default parameters. MiRmap [66] is a web-based application [67] that combines many thermodynamic, evolutionary, probabilistic and sequence-based features.

## Additional files

**Additional file 1: Table S1.** Number of sequence reads in each small non coding RNAs library. (DOCX 14 kb)

**Additional file 2: Excel file S1.** Predictions of miR genes by mirDEEP2 in SfC and SfR and orthology table. (XLSX 92 kb)

**Additional file 3: Figure S1.** MA-plots showing the relative expression of known or novel miR according to the host-plant (Rice compared to corn) in each strain. Top panel, in SfC, bottom panel In SfR. (PPTX 5002 kb)

**Additional file 4: Figure S2.** Variation between samples (treatments, replicates) of larvae exposed to different plants displayed by Principal Component Analysis (PCA). Top panel: SfC, bottom panel: SfR. (PPTX 4001 kb)

**Additional file 5: Excel file S2.** Relative expression resulting from DESEQ2 analysis within strains according to the host-plant. In red novel miR, in black known miRs. (XLSX 53 kb)

**Additional file 6: Figure S3.** MA-plots showing the relative expression of known or novel miR according to the genetic background. Top panel, relative expression analyzed by DESEQ2 in SfR compared to SfC on corn, bottom panel, relative expression in SfR compared to SfC on rice. (PPTX 4997 kb)

**Additional file 7: Figure S4.** Variation between samples (treatments, replicates) of larvae exposed to different plants displayed by Principal Component Analysis (PCA). Top panel: SfR compared to SfC on corn, bottom panel: SfR compared to SfC on rice. (PPTX 3008 kb)

**Additional file 8: Excel file S3.** Relative expression resulting from DESEQ2 analysis according to the genetic background on the same host-plant. In red novel miR, in black known miRs. (XLSX 39 kb)

**Additional file 9: Excel file S4.** List of gene target predictions of known differentially expressed miRs. (XLSX 1578 kb)

**Additional file 10: Excel file S5.** Target genes of known miRs and their relative expression in SfC on rice compared to corn. (XLSX 35 kb)

**Additional file 11: Excel file S6.** Target genes of known miRs and their relative expression in SfR compared to SfC on corn. (XLSX 19 kb)

## Acknowledgments

We are grateful to Clotilde Gibard and Gaëtan Clabots for maintaining the insect collections of the DGIMI laboratory in Montpellier. We thank Clotilde Gibard for her help to rear Spodoptera frugiperda on whole plants during reciprocal transplant experiments.

## Funding

This work was supported by a grant from the French National Research Agency (ANR-12-BSV7-0004-01; http://www.agence-nationale-recherche.fr/) for E. A. including a post-doctoral fellowship for Y.M. H.P. acknowledges financial support from the "France Génomique" National infrastructure, funded as part of "Investissement d'Avenir" (ANR-10-INBS-09). The funding bodies had no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

## Availability of data and materials

The datasets generated during this study are available from NCBI SRA (SRP132165).

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP132165

https://www.ncbi.nlm.nih.gov/sra/?term=SRP132165

The Genbank accession numbers of microRNA precursors identified in this study are MH340548 - MH341001.

## Authors' contributions

Conception and design: YM, NN, EA Acquisition of data, analysis and interpretation: YM, SN, SG, FL, IS, HP Drafting and revision of the manuscript: YM, NN, EA. All authors have read and approved the manuscript.

## Ethics approval and consent to participate

Not applicable.

## Consent for publication

Not applicable.

## Competing interests

The authors declare that they have no competing interests.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Author details

¹DGIMI, Univ Montpellier, INRA, Montpellier, France. ²INRA, UMR Institut de Génétique, Environnement et Protection des Plantes (IGEPP), Bioinformatics Platform for Agroecosystems Arthropods (BIPAA), Campus Beaulieu, Rennes, France. ³INRIA, IRISA, Genscale, Campus Beaulieu, Rennes, France. ⁴MGX, CNRS, INSERM, Univ Montpellier, Montpellier, France.

Received: 12 September 2018 Accepted: 26 September 2018

Published online: 06 November 2018

## References

1. Simon JC, d'Alencon E, Guy E, Jacquin-Joly E, Jaquieri J, Nouhaud P, Peccoud J, Sugio A, Streiff R. Genomics of adaptation to host-plants in herbivorous insects. Brief Funct Genomics. 2015;14(6):413-23.

2. Xu W, Papanicolaou A, Zhang HJ, Anderson A. Expansion of a bitter taste receptor family in a polyphagous insect herbivore. Sci Rep. 2016;6:23666.

3. Pearce SL, Clarke DF, East PD, Elfekih S, Gordon KHJ, Germiin LS, McGaughran A, Oakeshott JG, Papanikolaou A, Perera OP, et al. Genomic innovations, transcriptional plasticity and gene loss underlying the evolution and divergence of two highly polyphagous and invasive Helicoverpa pest species. BMC Biol. 2017;15(1):63.

4. Gouin A, Bretaudeau A, Nam K, Gimenez S, Aury JM, Duvic B, Hilliou F, Durand N, Montagné N, Darboux I, et al. Two genomes of highly polyphagous lepidopteran pests (Spodoptera frugiperda, Noctuidae) with different host-plant ranges. Sci Rep. 2017;7(1):11816.

5. Cheng T, Wu J, Wu Y, Chilukuri RV, Huang L, Yamamoto K, Feng L, Li W, Chen Z, Guo H, et al. Genomic adaptation to polyphagy and insecticides in a major east Asian noctuid pest. Nat Ecol Evol. 2017;1(11):1747-56.

6. Mathers TC, Chen Y, Kaithakottil G, Legeai F, Mugford ST, Baa-Puyoulet P, Bretaudeau A, Clavijo B, Colella S, Collin O, et al. Rapid transcriptional plasticity of duplicated gene clusters enables a clonally reproducing aphid to colonise diverse plant species. Genome Biol. 2017;18(1):27.

7. Orsucci M, Moné Y, Audiot P, Gimenez S, Nhim S, Nait-Saidi R, Frayssinet M, Dumont G, Pommier A, Boudon JP, et al. Transcriptional plasticity evolution in two strains of Spodoptera frugiperda (Lepidoptera: Noctuidae) feeding on alternative host-plants: BioRxiv; 2018. https://doi.org/10.1101/263186.

8. Celorio-Mancera MDLP, Heckel DG, Vogel H. Transcriptional analysis of physiological pathways in a generalist herbivore: responses to different host plants and plant structures by the cotton bollworm, Helicoverpa armigera. Entomol Exp Appl. 2012;144:123-33.

9. de la Paz Celorio-Mancera M, Wheat CW, Vogel H, Söderlind L, Janz N, Nylin S. Mechanisms of macroevolution: polyphagous plasticity in butterfly larvae revealed by RNA-Seq. Mol Ecol. 2013;22:4884-95.

10. Koenig C, Bretschneider A, Heckel DG, Grosse-Wilde E, Hansson BS, Vogel H. The plastic response of Manduca sexta to host and non-host plants. Insect Biochem Mol Biol. 2015;63:72-85.

11. Cora D, Re A, Caselle M, Bussolino F. MicroRNA-mediated regulatory circuits: outlook and perspectives. Phys Biol. 2017;14(4):045001.

12. Moran Y, Agron M, Praher D, Technau U. The evolutionary origin of plant and animal microRNAs. Nat Ecol Evol. 2017;1(3):27.

13. Ghildiyal M, Zamore PD. Small silencing RNAs: an expanding universe. Nat Rev Genet. 2009;10(2):94-108.
14. Catalanotto C, Cogoni C, Zardo G. MicroRNA in Control of Gene Expression: An Overview of Nuclear Functions. Int J Mol Sci. 2016;17(10):1712.

15. Palazzo AF, Lee ES. Non-coding RNA: what is functional and what is junk? Front Genet. 2015;6:2.

16. Arteaga-Vazquez M, Caballero-Perez J, Vielle-Calzada JP. A family of microRNAs present in plants and animals. Plant Cell. 2006;18(12):3355-69.

17. Friedlander MR, Chen W, Adamidi C, Maaskola J, Einspanier R, Knespel S, Rajewsky N. Discovering microRNAs from deep sequencing data using miRDeep. Nat Biotechnol. 2008;26(4):407-15.

18. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37-52.

19. Ma Z, Coruh C, Axtell MJ. Arabidopsis lyrata small RNAs: transient MIRNA and small interfering RNA loci within the Arabidopsis genus. Plant Cell. 2010;22(4):1090-103.

20. Pradhan B, Naqvi AR, Saraf S, Mukherjee SK, Dey N. Prediction and characterization of tomato leaf curl New Delhi virus (ToLCNDV) responsive novel microRNAs in Solanum lycopersicum. Virus Res. 2015;195:183-95.

21. Kakumani PK, Chinnappan M, Singh AK, Malhotra P, Mukherjee SK, Bhatnagar RK. Identification and characteristics of microRNAs from army worm, Spodoptera frugiperda cell line Sf21. PLoS One. 2015;10(2):e0116988.

22. Mehrabadi M, Hussain M, Asgari S. MicroRNAome of Spodoptera frugiperda cells (Sf9) and its alteration following baculovirus infection. J Gen Virol. 2013;94(Pt 6):1385-97.

23. Lau NC, Lim LP, Weinstein EG, Bartel DP. An abundant class of tiny RNAs with probable regulatory roles in Caenorhabditis elegans. Science. 2001; 294(5543):858-62.

24. Kertesz M, Iovino N, Unnerstall U, Gaul U, Segal E. The role of site accessibility in microRNA target recognition. Nat Genet. 2007;39(10):1278-84.

25. Zuker M. Mfold web server for nucleic acid folding and hybridization prediction. Nucleic Acids Res. 2003;31(13):3406-15.

26. Sarov-Blat L, So WV, Liu L, Rosbash M. The Drosophila takeout gene is a novel molecular link between circadian rhythms and feeding behavior. Cell. 2000;101(6):647-56.

27. Andrade Lopez JM, Lanno SM, Auerbach JM, Moskowitz EC, Sligar LA, Wittkopp PJ, Coolon JD. Genetic basis of octanoic acid resistance in Drosophila sechellia: functional analysis of a fine-mapped region. Mol Ecol. 2017;26(4):1148-60.

28. Smith CR, Morandin C, Noureddine M, Pant S. Conserved roles of Osiris genes in insect development, polymorphism, and protection. J Evol Biol. 2018;31(4):S16-S29.

29. Kanamori Y, Saito A, Hagiwara-Komoda Y, Tanaka D, Mitsumasu K, Kikuta S, Watanabe M, Cornette R, Kikawada T, Okuda T. The trehalose transporter 1 gene sequence is conserved in insects and encodes proteins with different kinetic properties involved in trehalose import into peripheral tissues. Insect Biochem Mol Biol. 2010;40(1):30-7.

30. Fernandez-Sanchez E, Diez-Guerra FJ, Cubelos B, Gimenez C, Zafra F. Mechanisms of endoplasmic-reticulum export of glycine transporter-1 (GLYT1). Biochem J. 2008;409(3):669-81.

31. Miller-Fleming L, Olin-Sandoval V, Campbell K, Ralser M. Remaining mysteries of molecular biology: the role of polyamines in the cell. J Mol Biol. 2015;427(21):3389-406.

32. Eckert D, Buhl S, Weber S, Jager R, Schorle H. The AP-2 family of transcription factors. Genome Biol. 2005;6(13):246.

33. Bauer R, McGuffin ME, Mattox W, Tainsky MA. Cloning and characterization of the Drosophila homologue of the AP-2 transcription factor. Oncogene. 1998;17(15):1911-22.

34. Monge I, Mitchell PJ. DAP-2, the Drosophila homolog of transcription factor AP-2. Mech Dev. 1998;76(1-2):191-5.

35. Kerber B, Monge I, Mueller M, Mitchell PJ, Cohen SM. The AP-2 transcription factor is required for joint formation and cell survival in Drosophila leg development. Development. 2001;128(8):1231-8.

36. Monge I, Krishnamurthy R, Sims D, Hirth F, Spengler M, Kammermeier L, Reichert H, Mitchell PJ. Drosophila transcription factor AP-2 in proboscis, leg and brain central complex development. Development. 2001;128(8):1239-52.

37. Visentin M, Diop-Bove N, Zhao R, Goldman ID. The intestinal absorption of folates. Annu Rev Physiol. 2014;76:251-74.

38. Zhao R, Goldman ID. The proton-coupled folate transporter: physiological and pharmacological roles. Curr Opin Pharmacol. 2013;13(6):875-80.

39. Tan SM, Lieberman J. Capture and identification of miRNA targets by biotin Pulldown and RNA-seq. Methods Mol Biol. 2016;1358:211-28.

40. Thomson DW, Bracken CP, Goodall GJ. Experimental strategies for microRNA target identification. Nucleic Acids Res. 2011;39(16):6845-53.

41. Soni K, Choudhary A, Patowary A, Singh AR, Bhatia S, Sivasubbu S, Chandrasekaran S, Pillai B. miR-34 is maternally inherited in Drosophila melanogaster and Danio rerio. Nucleic Acids Res. 2013;41(8):4470-80.

42. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10:57-63.

43. Lamarre S, Frasse P, Zouine M, Labourdette D, Sainderichin E, Hu G, Le Berre-Anton V, Bouzayen M, Maza E. Optimization of an RNA-Seq differential gene expression analysis depending on biological replicate number and library size. Front Plant Sci. 2018;9:108.

44. Garcia-Segura L, Abreu-Goodger C, Hernandez-Mendoza A, Dimitrova Dinkova TD, Padilla-Noriega L, Perez-Andrade ME, Miranda-Rios J. High-throughput profiling of Caenorhabditis elegans starvation-responsive microRNAs. PLoS One. 2015;10(11):e0142262.

45. Lyons PJ, Storey KB, Morin P Jr. Expression of miRNAs in response to freezing and anoxia stresses in the freeze tolerant fly Eurosta solidaginis. Cryobiology. 2015;71(1):97-102.

46. Campos C, Sundaram AY, Valente LM, Conceicao LE, Engrola S, Fernandes JM. Thermal plasticity of the miRNA transcriptome during Senegalese sole development. BMC Genomics. 2014;15:525.

47. Yin F, Qin C, Gao J, Liu M, Luo X, Zhang W, Liu H, Liao X, Shen Y, Mao L, et al. Genome-wide identification and analysis of drought-responsive genes and microRNAs in tobacco. Int J Mol Sci. 2015;16(3):5714-40.

48. Voskarides K. Plasticity vs mutation. The role of microRNAs in human adaptation. Mech Ageing Dev. 2017;163:36-9.

49. Rodgers AB, Morgan CP, Leu NA, Bale TL. Transgenerational epigenetic programming via sperm microRNA recapitulates effects of paternal stress. Proc Natl Acad Sci U S A. 2015;112(44):13699-704.

50. Rassoulzadegan M, Grandjean V, Gounon P, Vincent S, Gillot I, Cuzin F. RNA-mediated non-mendelian inheritance of an epigenetic change in the mouse. Nature. 2006;441:469-74.

51. Zhang YL, Huang QX, Yin GH, Lee S, Jia RZ, Liu ZX, Yu NT, Pennerman KK, Chen X, Guo AP. Identification of microRNAs by small RNA deep sequencing for synthetic microRNA mimics to control Spodoptera exigua. Gene. 2015;557(2):215-21.

52. Nandety RS, Kuo YW, Nouri S, Falk BW. Emerging strategies for RNA interference (RNAi) applications in insects. Bioengineered. 2015;6(1):8-19.

53. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBioNetjournal. 2011;17(1):10-2.

54. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25(17):3389-402.

55. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.

56. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.

57. Pfaffl MW, Horgan GW, Dempfle L. Relative expression software tool (REST) for group-wise comparison and statistical analysis of relative expression results in real-time PCR. Nucleic Acids Res. 2002;30(9):e36.

58. Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of mammalian microRNA targets. Cell. 2003;115(7):787-98.

59. Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120(1):15-20.

60. TargetScanHuman Prediction of microRNA targets. https://genesmitedu/targetscan.

61. Hofacker IL. Vienna RNA secondary structure server. Nucleic Acids Res. 2003; 31(13):3429-31.

62. BIPAA Bioinformatics Platform for Agroecosystem Arthropods. In: https://bipaagenouestorg/is/lepidodb/spodoptera_frugiperda/

63. Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS. MicroRNA targets in Drosophila. Genome Biol. 2003;5(1):R1.

64. John B, Enright AJ, Aravin A, Tuschl T, Sander C, Marks DS. Human MicroRNA targets. PLoS Biol. 2004;2(11):e363.

65. Miranda KC, Huynh T, Tay Y, Ang YS, Tam WL, Thomson AM, Lim B, Rigoutsos I. A pattern-based method for the identification of MicroRNA binding sites and their corresponding heteroduplexes. Cell. 2006;126(6):1203-17.

66. Vejnar CE, Blum M, Zdobnov EM. miRmap web: comprehensive microRNA target prediction online. Nucleic Acids Res. 2013;41(Web Server issue):W165-8.

67. MiRmap. In. https://mirmap.ezlab.org.