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

Article

# Development of a Species Diagnostic Molecular Tool for an Invasive Pest, *Mythimna loreyi*, Using LAMP

Hwa Yeun Nam<fnref n="1" />, Min Kwon<fnref n="1" />, Hyun Ju Kim<fnref n="2" /> and Juil Kim<fnref n="1" />,<fnref n="3" />,*

Received: 20 October 2020; Accepted: 17 November 2020; Published: 19 November 2020

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

check for
updates

**Simple Summary:** *Mythimna loreyi* is a serious pest of grain crops and reduces yields in maize plantations. *M. loreyi* is a native species in East Asia for a long time. However, this species has recently emerged as a migration (or invasive from other Asian countries) pest of some cereal crops in Korea. Little is known about its basic biology, ecology, and it is difficult to identify the morphologically similar species, *Mythimna separate*, which occur at the cornfield in the larvae stage. Species diagnosis methods for invasive pests have been developed and utilized for this reason. Currently, the molecular biology method for diagnosing *M. loreyi* species is only using the mtCO1 universal primer (LCO1490, HCO2198) and process PCR and sequencing to compare the degree of homology. However, this method requires a lot of time and effort. In this study, we developed the LAMP (loop-mediated isothermal amplification) assay for rapid, simple, and effective species diagnosis. By analyzing the mitochondrial (mt) genome, the species-specific sequence was found at the coding region of the NADH dehydrogenase subunit 5 gene. A broad range of DNA concentration was workable in LAMP assay, in which the minimum detectable DNA concentration was 100 pg. DNA releasing method was applied, which took five minutes of incubation at 95 °C without the DNA extraction process, and only some pieces of tissue from larvae and adult samples were needed. The incidence of invasive pests is gradually diversifying. Therefore, this simple and accurate LAMP assay is possibly used in the intensive field monitoring for invasive pests and integrated management of *Mythimna loreyi*.

**Abstract:** The *Mythimna loreyi* (Duponchel) is one of the well-known invasive noctuid pests in Africa, Australia, and many Asian countries. However, it is difficult to identify the invasive and morphologically similar species, *Mythimna separate*, which occur at the cornfield in the larvae stage. Currently, the molecular biology method for diagnosing *M. loreyi* species is only using the mtCO1 universal primer (LCO1490, HCO2198), which requires a lot of time and effort, such as DNA extraction, PCR, electrophoresis, and sequencing. In this study, the LAMP assay was developed for rapid, simple, effective species identification. By analyzing the mitochondrial genome, the species-specific sequence was found at the coding region of the NADH dehydrogenase subunit 5 gene. Based on this unique sequence, four LAMP primers and two loop primers were designed. The F3 and B3 primers were able to diagnose species-specific, in general, and multiplex PCR and specifically reacted within the inner primers in LAMP assay. The optimal incubation condition of the LAMP assay was 61 °C for 60 min with four LAMP primers, though additional loop primers, BF and LF, did not significantly shorten the amplification time. The broad range of DNA concentration was workable in LAMP assay, in which the minimum detectable DNA concentration was 100 pg. DNA releasing method was applied, which took five minutes of incubation at 95 °C without the DNA extraction process. Only some pieces of tissue of larvae and adult samples were needed to extract DNA. The incidence of invasive pests is

<footnote n="1">1 Highland Agriculture Research Institute, National Institute of Crop Science, Rural Development Administration, Pyeongchang 25342, Korea; jessienam89@gmail.com (H.Y.N.); mkwon@korea.kr (M.K.)</footnote>

<footnote n="2">2 Crop foundation Division, National Institute of Crop Science, Rural Development Administration, Wanju 55365, Korea; yaehyunj@korea.kr</footnote>

<footnote n="3">3 Program of Applied Biology, Division of Bio-resource Sciences, College of Agriculture and Life Science, Kangwon National University, Chuncheon 24341, Korea</footnote>

<footnote>* Correspondence: forweek@kangwon.ac.kr; Tel.: +82-33-250-6438</footnote>
gradually diversifying. Therefore, this simple and accurate LAMP assay is possibly applied in the
intensive field monitoring for invasive pests and integrated management of *Mythimna loreyi*.

**Keywords: Mythimna loreyi; rice armyworm; invasive pest; LAMP; diagnostic PCR**

**1. Introduction**

The *Mythimna loreyi* (Duponchel) (often called the cosmopolitan) is a noctuid pest of grain crops found in Africa, Australia, the Near East, and the Middle East and undergoes multiple generations per year [1–4]. *M. loreyi* feeds on various host plants, including rice, wheat, maize, sugarcane, barley, sorghum, and others, which have a large effect on female fecundity. The fecundity of female moths is greatest when the larvae feed on maize in Egypt [5]. Since *M. loreyi* is facilitated to breed, some researchers focus on the identification of products secreted by the adult corpora allata [6]. Not only the physiological understanding of this species [7–10] but also ecology-based developmental characteristics [11] and flight performance [12] have been studied. For the biological control, *M. loreyi* densovirus (MIDNVs) was isolated in Egypt and characterized [13,14]. The outbreak of this pest and the damage to crops have been proliferated, particularly in some Asian countries. In Japan, *M. loreyi* typically occurs together with *M. separate* (Walker) that has significant negative effects on crop production [15]. Besides, *M. loreyi* has begun to occur and damage host plants together with *M. separate* [16].

During a couple of years in 2019 and 2020, there have been reports that cornfield has been damaged by the larvae of *M. loreyi* in Korea. This indicates that there is a possibility that *M. loreyi* can change into a sporadic pest, which can cause serious damage to crops. As per these cases, *M. loreyi* is occasionally damaged with its sister species, *M. separate*. However, it is difficult to distinguish two species at the cornfield in the larvae stage. Only one molecular diagnostic tool has been studied to distinguish *M. loreyi*, which is based on the sequencing of part of the mitochondrial COI gene [17]. However, only one mutation exists within the 658 bp amplicon. Moreover, high sequence similarity has been shown between *M. loreyi* and *M. separata*, which indicates the limitation of diagnosis of this pest.

In this study, we developed a simpler technique, termed loop-mediated isothermal amplification assay (LAMP). It is also widely used for the rapid and accurate identification of pest species [18–20]. The LAMP is a rapid, simple, effective, and specific amplification of DNA compared to real-time PCR based on the mitochondrial gene. It is performed under isothermal conditions that require a set of four primers, a strand-displacing DNA polymerase, and a water bath or heat block to maintain the temperature at about 65 °C following a one-time denaturation at 95 °C [21] or one-step incubation at about 65 °C [22].

Following the first infestation of *M. loreyi* in Korea, there is great demand from agricultural research, extension services, and farmers for diagnostic methods for these species. Therefore, we present a LAMP-based method for specimens collected in Korea and other sequences from GenBank. This method should be useful in assisting the effective pest management of *M. loreyi*.

**2. Materials and Methods**

**2.1. *Sample Collection and Mitochondrial Genome Sequencing***

The larval stage of *Mythimna loreyi* Korean populations was collected from Hadong (35°02′17″ N, 127°47′12″ E) in a cornfield, 2019. Some larvae reared in the lab for morphological conformation in the adult stage, and the genomic DNA of several of the individual larva was directly extracted with DNAzol (Molecular Research Center, Cincinnati, OH, USA) and quantified by Nanodrop (NanoDrop Technologies, Wilmington, DE, USA). Besides, the genomic DNA of over 20 larvae or adults was extracted (population pooled genomic DNA). Universal primers (LCO1490 and HCO2198) were used
with *M. loreyi*'s individual and pooled genomic DNA as templates in 20 µL PCR reaction containing 1 U TOYOBO KOD—FX TaqTM (Toyobo Life Science, Osaka, Japan), 2X buffer (with 15 mM MgCl2), 0.2 mM each dNTP, 0.5 µM each primer, and 100 ng genomic DNA [23]. The PCR products were directly sequenced (chromatogram) to verify the nucleotide polymorphism, and no mutation was found in intraspecies. For mitochondrial genome sequencing, the Miseq platform was used, and more than 1 Gb was sequenced. To assemble these data, the CLC Assembly Cell package (version 4.2.1) was used. After trimming raw data using CLC quality trim (ver. 4.21), the assembly was accomplished using the CLC de novo assembler with dnaLCW. Assembled sequences were confirmed by BLASTZ [24]. The GeSeq program was used for annotation [25], and the result was manually checked based on the alignment of other Noctuidae species mitochondrial genomes using MEGA 7 [26].

## 2.2. Phylogenetic Analysis and Primer Design

Molecular phylogenetic analysis of mitochondrion genomes was inferred by using the maximum likelihood method implemented by MEGA 7 with bootstrapping [26,27]. Mitochondrial genome sequences of other Noctuidae species were downloaded from GenBank, NCBI. For comparative analysis, mitochondrial genomes were aligned using mVISTA [28,29]. Based on the global alignment result, partial sequences were re-aligned for LAMP primer design using PrimerExplorer V5.

## 2.3. LAMP and PCR

WarmStart® LAMP Kit (New England Biolabs, Ipswich, UK) was used for the LAMP assay. The general protocol of LAMP was performed following the manufacture's guidelines in a 25 µL reaction mixture. For the general PCR, TOYOBO KOD—FX TaqTM (Toyobo Life Science, Osaka, Japan) was used in this study. Appropriate primers were used with the following PCR amplification protocol: a 2 min denaturing step at 94 °C and a PCR amplification cycle consisting of denaturing at 94 °C for 20 s, annealing at 60 °C for 20 s, and extension at 68 °C for 30 s, which was repeated 35 times. The amplified DNA fragments were separated using 1.5% agarose gel electrophoresis and visualized with SYBR green (Life Technologies, Grand Island, NY, USA). The larval stage of *M. loreyi* samples, which were collected from Hadong (Korea) in a cornfield (2019), was used. Pheromone traps were used for adults' sample collection, such as *M. separate*, *Agrotis segetum*, *Spodoptera frugiperda*, *Spodoptera exigua*, *Spodoptera litura*, and *Helicoverpa armigera*. Traps were set in Pyeongchang (37°40'53" N, 128°43'49" E), Hongchen (37°43'35" N, 128°24'33" E), and Gangneung (37°36'56" N, 128°45'59" E) [30]. DNA samples were prepared using DNAzol (Molecular Research Center, Cincinnati, OH, USA) from trapped adults. We used universal primers (LCO1490 and HCO2198) with each DNA sample and pooled genomic DNA as templates in a 20 µL PCR reaction containing 1 U TOYOBO KOD—FX TaqTM (Toyobo Life Science, Osaka, Japan), 2X buffer (with 15 mM MgCl2), 0.2 mM each dNTP, 0.5 µM each primer, and 100 ng genomic DNA [23]. The PCR products were directly sequenced to verify the mutation of each population using RAW_F3S and RAW_B3 primer set (MACROGEN). Three biological DNA samples that were collected from the fields were used in each LAMP and PCR to validate the reliability of the LAMP condition.
3. Results

3.1. *Mitochondrial Genome Sequencing and Primer Design*

A 15,320 bp of mitochondrial genome verified after trimming from about 2.2 Gb (7,312,504 reads) nucleotide sequences information was obtained through *Miseq*. The mitochondrial genome of *Mythimna loreyi* was assembled (GenBank MT506351). The mitochondrial genome included 13 protein-coding genes: NADH dehydrogenase components (complex I, ND), cytochrome oxidase subunits (complex VI, COX), cytochrome oxidase b (CYPB) and two ATP synthases, and two ribosomal RNA genes and 22 transfer RNAs (Supplementary Figure S1).

As a result of MegaBLAST, the most homologous species was an allied species, *Mythimna separate*, which showed 93.7% similarity. The genus of Spodoptera [31], such as *Spodoptera exigua*, *S. litura*, *S. frugiperda*, which are possibly found together with *M. loreyi* in the cornfields, showed about 89 to 90% homology based on the mitochondrial genome sequence. *Agrotis segetum* [32], which, in particular, occurs and damage at the corn seedling stage, showed 90.6% similarity to *M. loreyi* based on mitochondrial genome sequence (data not shown).

The phylogenetic relationship between 15 mitochondrial genomes of 14 species was examined (Figure 1A) to verify a specific nucleotide sequence that only *M. loreyi* possessed among the related species with high gene similarity or a morphologically similar pest. The result of the phylogenetic relationship was mostly similar to the megaBLAST result. Based on the mVISTA alignment results, conserved regions among Noctuidae species and variable regions were observed (Figure 1B). By combining the two results, the partial sequence in eight species of nine ND5 mitochondrial genome was re-aligned to design the specific primer of the *M. loreyi*. Finally, four essential primers and two loop primers were designed (Figure 2 and Table 1).

<figure><img src="image_36.png" /><figcaption>**Figure 1.** Comparison of entire mitochondrial genomes of some Noctuidae pests, including newly sequenced *Mythimna loreyi*. (A) Phylogenetic relationship inferred using maximum likelihood under MEGA7. (B) Schematic diagram of the genes and their flanking regions, showing the sequence diversity in mVISTA. UTR, D-loop denotes untranslated region and displacement-loop, respectively. Eight green colored mitochondrial genome sequences were re-aligned for primer design with that of the target species, *M. loreyi* (Figure 2).</figcaption></figure>
<figure><img src="image_37.png" /></figure>

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

<figure><img src="image_39.png" /><figcaption>**Figure 2.** Location of primers and primer binding regions on partial sequence of some Noctuidae pests' mtDNA for species identification of *M. loreyi*. Inner primer, FIP, consists of F1c (complementary sequences of F1) and F2. Another inner primer, BIP, is also composed of B1 and B2c (complementary sequences of B2). Essential four LAMP primers (F3, FIP, BIP, and B3) generate the dumbbell structure, and two loop primers, LF and LB, possibly accelerate the LAMP (loop-mediated isothermal amplification) reaction. Used primer information is documented in Table 1.</figcaption></figure>

<table><caption>**Table 1.** Primer list for LAMP (loop-mediated isothermal amplification) and PCR in this study.</caption><thead><tr><th>Purpose</th><th>Primers</th><th>Sequence (5' – > 3')</th></tr></thead><tbody><tr><td colspan="3">for LAMP</td></tr><tr><td>RAW_F3</td><td>GTAAATTTATTACAAGAATAAATCCCC 1</td><td></td></tr><tr><td>RAW_B3</td><td>CTTCTACTTTAGAACTGCGGGA</td><td></td></tr><tr><td>RAW_FIP</td><td>TGGGGTAATTTCATATATAAATGATAATATATAATCTAATTCCCCCTATAAAACG</td><td></td></tr><tr><td>RAW_BIP</td><td>TAGCATGAGTTAATAAATGAAAAAAAGGTTTAATAATAAGAATTTTAAGAATGGGT</td><td></td></tr><tr><td>RAW_LF</td><td>GTGTATATCATTATATGATAATACCTCGCAGATATAATAAAGCCTTTAATATAGCATGAGTATAAATGAAAAAGAAAGCTAATCAGGCATCTTCTTATTAATTCTTATTAATACCAATTAATGTTCTTAAGAGATAGAAAGCATAATCTTT</td><td></td></tr><tr><td>RAW_LB</td><td>TTGTATATCATTATATGATAATACCTCGCAGATATAATAAAGCCTTTAATATAGCATGAGTATAAATGAAAAAGAAAGCTAATCAGGCATCTTCTTATTAATTCTTATTAATACCAATTAATGTTCTTAAGAGATAGAAAGCATAATCTTT</td><td></td></tr><tr><td colspan="3">for PCR</td></tr><tr><td>Spo_ace1UF2</td><td>AGGATGAAGAGAAATTTATAGAGGAT</td><td></td></tr><tr><td>Spo_ace1UR1</td><td>TCACCAAACACTGTATCTATAATTTG</td><td></td></tr><tr><td>LCO1490</td><td>GGTCAACAAATCATAAAGATATTGG</td><td></td></tr><tr><td>HCO2198</td><td>TAAACTTCAGGCTGACCAAAAAATCA</td><td></td></tr><tr><td>RAW_F3S</td><td>CCCAAACCCTCTATATAATTCTCT</td><td></td></tr></tbody></table>

Among the six primers, F3 is the specific primer that enables the diagnosis of *M. loreyi*, and only *M. loreyi* had the CCCC sequences in four priming regions, which are marked in the red box. A total of four populations that were collected from different regions had the same nucleotide sequence. Therefore, it could be sufficiently used for species diagnosis. The basic species diagnosis primer production strategy is the same as the previously reported species diagnosis development method of *S. frugiperda* [33]. We targeted the open reading frame (ORF) region because the ORF region is more intraspecies conserved than the non-coding region.

<footnote n="1">¹ 'G's at the 5' end depicted in red were added to adjust the primer melting temperature.</footnote>
3.2. *Diagnostic LAMP and PCR*

As previously reported, the sensitivity of LAMP may vary depending on the temperature and
reaction time [33]. Therefore, the reaction was performed at 65, 63, and 61 °C to find the optimal
reaction conditions (Figure 3).

<figure><img src="image_40.png" /><figcaption>**Figure 3.** The sensitivity of the LAMP assay results in three temperature conditions, such as 65, 63, and 61 °C, for *M. loreyi* species detected under (A1, B1, and C1) visible light, (A2, B2, and C2) ultraviolet light with SYBR Green, and (A3, B3, and C3) gel electrophoresis. The original pink color of the reaction mixture turned yellow in a positive reaction when the product was formed but remained pink in negative reactions. (D) Conventional and multiplex PCR to distinguish *M. loreyi*. The 794 bp amplicon was amplified only in *M. loreyi*, and the conserved partial sequence of ace1 type acetylcholinesterase gene was targeted as an internal reference. Abbreviations are NTC (non-template control), RAW (rice armyworm *Mythimna loreyi*), OAW (oriental armyworm *Mythimna separata*), TM (Turnip moth *Agrotis segetum*), CBW (cotton bollworm *Helicoverpa armigera*), BAW (beet armyworm *Spodoptera exigua*), FAW (fall armyworm *Spodoptera frugiperda*), and TCW (tobacco cutworm *Spodoptera litura*).</figcaption></figure>

As the amount of template DNA was quantified as 50 ng and reacted at each temperature with 25 µL reaction volume, the diagnosis result was confirmed in 120 min at 65 °C and 90 min at 63 °C. The diagnostic level of reaction did not occur when the reaction performed less than the corresponding time in each temperature. Despite the relatively low temperature, at 61 °C, the diagnostic level of reaction was confirmed in only 60 min, and the false-positive reaction did not appear only once in the results of more than three repetition tests. Despite the relatively low temperature, at 61 °C, the diagnostic level of reaction was confirmed in only 60 min, and the false-positive reaction did not appear only once in three repetition tests (Figure 3C). Under 100 bp sized bands were not the LAMP products and were aggregated primers and generated even in the negative control.

In the LAMP assay, to confirm the diagnosis primer F3, it was reacted with reverse primer B3 through PCR, and the 794 bp of PCR product was identified (Figure 3(D1)). Besides, a species-specific reaction was confirmed, as we performed the multiplex PCR with a universal positive control primer set that targets ace1-type acetylcholinesterase, which can produce a positive reaction in all samples (Figure 3(D2)). The two-loop primers were tested for possible enhancement of the LAMP reaction, as suggested by Nagamine et al. [22].

As a result of reacting two-loop primers with each or two together at 61 °C for 60 min, there was no difference between the reaction result of using only four primers and shortening the reaction time (Figure 4). Even when 10 µL of the reaction solution was verified by electrophoresis after the LAMP reaction, there was no difference in band intensity. It was possible to reliably diagnose up to 100 pg when reacting using four LAMP primers without adding a loop primer (Figure 5A). Therefore, *M. loreyi* species diagnosis LAMP method that we developed could diagnose under various DNA concentration conditions from 100 ng to 100 pg. To increase the usability in the field, we cut a part of the tissue of the adult antenna or leg and put it in 30 µL distilled water, which reacts at 95 °C for 5 min, securing
the template DNA without a separate DNA extraction process (Figure 5(B4)). As we measured each sample of DNA concentration obtained through this DNA releasing method, each sample showed various measurements. The species-specific diagnosis was possible in the same method as using the template DNA, which was obtained through separate DNA extraction as a positive control since it was in the range of LAMP reaction (Figure 5B).

<figure><img src="image_41.png" /><figcaption>**Figure 4.** The LAMP assay results with (A) 4 primers and additional loop primers, (B) loop forward, LF, or (C) loop backward, LB, or (D) two-loop primers, LF and LB under (A1, B1, C1, and D1) visible light, (A2, B2, C2, and D2) ultraviolet light with SYBR Green, and (A3, B3, C3, and D3) gel electrophoresis. Abbreviations are NTC (non-template control), RAW (rice armyworm *Mythimna loreyi*), OAW (oriental armyworm *Mythimna separata*), TM (Turnip moth *Agrotis segetum*), CBW (cotton bollworm *Helicoverpa armigera*), BAW (beet armyworm *Spodoptera exigua*), FAW (fall armyworm *Spodoptera frugiperda*), and TCW (tobacco cutworm *Spodoptera litura*).</figcaption></figure>

**4. Discussion**

Invasive pests, such as *S. frugiperda*, are increasing worldwide due to global warming and climate change [34,35]. Currently, *M. loreyi* is a serious grain pest in Africa, Australia, the Near East, and the Middle East and undergoes multiple generations per year [1–4]. Besides, *M. loreyi* is a species that originates from China and results in damage upon spreading gradually. In Korea, *M. loreyi* has been reported for a long time and maybe a native pest to Korea. Recently, there are several reports that the larvae of *M. loreyi* have damaged the cornfield in Korea. There is a possibility that *M. loreyi* can change into a sporadic pest, which can cause serious damage to crops in Korea.

As with other invasive pests found in the field, they are often similar in morphology to allied species, making the investigation of initial density and pest management difficult. Species diagnosis methods for invasive pests, such as *S. frugiperda*, have been developed and utilized for this reason [33]. Currently, the molecular biology method for diagnosing *M. loreyi* species is only using the mtCO1 universal primer (LCO1490, HCO2198) and process PCR and sequencing to compare the degree of homology. However, this method requires a lot of time and effort, such as DNA extraction, PCR, electrophoresis, and sequencing. Therefore, we developed a molecular diagnosis method that could diagnose species within a short time without a separate DNA extraction process, and only a heat block is needed that can control temperature (Figure 5). The basis diagnosis strategy is very similar to the *S. frugiperda* species diagnosis method [33]. This method can complete all experimental procedures and verify the results within 1 h and 30 min right after obtaining a sample. In this study, the results only specified adult samples, but it is possible to use larvae.
<figure><img src="image_42.png" /><figcaption>**Figure 5.** (A) Identification of the detection limit of genomic DNA in the LAMP assay from 100 ng to 100 fg under (A1 and B1) visible light, (A2 and B2) ultraviolet light with SYBR Green, and (A3 and B3) gel electrophoresis. (B) The sensitivity of the LAMP assay results with the DNA releasing technique from insect tissue. Around 10 mg of the adult leg (or antenna) was incubated at 95 °C for 5 min (B4). Abbreviations are NTC (non-template control), RAW (rice armyworm *Mythimna loreyi*), OAW (oriental armyworm *Mythimna separata*), TM (Turnip moth *Agrotis segetum*), CBW (cotton bollworm *Helicoverpa armigera*), BAW (beet armyworm *Spodoptera exigua*), FAW (fall armyworm *Spodoptera frugiperda*), and PC (positive control, isolated DNA from *M. loreyi*).</figcaption></figure>

The simplicity, accuracy, and adaptability for high throughput of the LAMP assay are distinct advantages [21,36,37]. Moreover, recently LAMP is utilized in various fields, such as many ecology studies, medical aspects, outside of the lab, and can be applied to diagnose plant viruses in insect body and insecticide-resistant gene mutation [38,39]. Besides, the diagnostic primer used for LAMP can be used for various diagnostic methods because it was possible to apply in general PCR and multiplex PCR (Figure 3). Therefore, it is feasible to diagnose a larger sample with positive control in the form of multiplex PCR, which is sufficiently modified and used in the laboratory. There are advantages and disadvantages to find a species diagnostic marker in the mitochondrial genome as well as in part of the genomic DNA. First of all, the disadvantage is that LAMP primer production is limited because there are many parts of AT-rich. It is complicated to design a primer that requires a minimum of four primers, and primers FIP and BIP are high-performance liquid chromatography (HPLC) purified primers. As an advantage, LAMP has higher amplification and efficiency and sensitivity compared to real-time PCR, and the results can be visually monitored by the naked eye either through turbidity or color change by fluorescent intercalating dye (Syber Green I) [40]. The number of copies is large, and it is possible to diagnose with a small amount of DNA or DNA releasing method, whose amplification can be accomplished with heating block. On the contrary, a species diagnosis marker in the genomic DNA should be found in exon rather than intron because it has a large variation. But it is difficult to develop marker in exon because it is often quite conserved within an allied species. Besides, the DNA releasing method can be used, but the efficiency is low when the copies of the gene are small [41].
**5. Conclusions**

In this study, a species diagnosis marker was designed within the mitochondrial genome to combine it with a DNA releasing method that is highly applicable to diagnose species in the field. Moreover, a significantly efficient method was developed, which targeted the ace1 gene as a positive control of the species to compare. This simple and accurate diagnosis using LAMP assay could be possibly applied in the intensive field to monitor and for the pest management of *M. loreyi*.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-4450/11/11/817/s1, Figure S1: Organization of the mitochondrial genome of *Mythimna loreyi* from Korea (GenBank MT506351). ND: NADH dehydrogenase components (Complex I) in yellow. COX: cytochrome oxidase subunits (Complex VI) in pink. ATP synthase in green. CYPB: cytochrome oxidase b in purple. Ribosomal RNA genes in red, tRNA genes in blue. Noncoding regions are not colored.

**Author Contributions:** Conceptualization, H.Y.N. and J.K.; methodology, H.Y.N. and J.K.; software, J.K.; validation, J.K.; formal analysis, J.K.; investigation, H.Y.N. and J.K.; resources, H.Y.N., M.K., H.J.K., and J.K.; data curation, H.Y.N. and J.K.; writing—Original draft preparation, H.Y.N. and J.K.; writing—Review and editing, H.Y.N. and J.K.; visualization, H.Y.N. and J.K.; supervision, J.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** The Cooperative Research Program supported this study for Agriculture Science and Technology Development (Project No. PJ01509301), the Rural Development Administration, Republic of Korea.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Calora, F.B. Revision of the species *Leucania*-complex occurring in the Philippines (Lepidoptera, Noctuidae, Hadeninae). *Philipp. Agric.* **1966**, *50*, 633–723.

2. Chandler, K.J.; Benson, A.J. Evaluation of armyworm infestation in North Queensland sugarcane on crops. *Proc. Aust. Soc. Sug. Cane Technol.* **1991**, *13*, 79–82.

3. Edwards, E.D. A second sugarcane armyworm (*Leucania loreyi* (Duponchel) from Australia and the identity of *L. loreyimima* Rungs (Lepidoptera: Noctuidae)). *J. Aust. Entomol. Soc.* **1992**, *31*, 105–108. [CrossRef]

4. Ganesha, S.; Rajabale, S. The *Mythimna* spp. (Lepidoptera: Noctuidae) complex on sugarcane in Mauritius. *Proc. S. Afr. Sug. Technol. Ass.* **1996**, *70*, 15–17.

5. El-Sherif, S.I. On the biology of *Leucania loreyi*, dup. (Lepidoptera, Noctuidae). *J. Appl. Entomol.* **1972**, *71*, 104–111. [CrossRef]

6. Ho, H.Y.; Tu, M.P.; Chang, C.Y.; Yin, C.M.; Kou, R. Identification of in vitro release products of corpora allata in female and male loreyi leafworms, *Leucania loreyi*. *Experientia* **1995**, *51*, 601–605. [CrossRef]

7. Hsieh, Y.C.; Hsu, E.L.; Chow, Y.S.; Kou, R. Effects of calcium channel antagonists on the corpora allata of adult male loreyi leafworm *Mythimna loreyi*: Juvenile hormone acids release and intracellular calcium level. *Arch. Insect Biochem. Physiol.* **2001**, *48*, 89–99. [CrossRef]

8. Hsieh, Y.C.; Yang, E.C.; Hsu, E.L.; Chow, Y.S.; Kou, R. Voltage-dependent calcium channels in the corpora allata of the adult male loreyi leafworm, *Mythimna loreyi*. *Insect Biochem. Mol. Biol.* **2002**, *32*, 547–557. [CrossRef]

9. Kou, R. Cholinergic regulation of the corpora allata in adult male loreyi leafworm *Mythimna loreyi*. *Arch. Insect Biochem. Physiol.* **2002**, *49*, 215–224. [CrossRef]

10. Kou, R.; Chen, S.J. Allatotropic activity in the suboesophageal ganglia and corpora cardiaca of the adult male loreyi leafworm, *Mythimna loreyi*. *Arch. Insect Biochem. Physiol.* **2000**, *43*, 78–86. [CrossRef]

11. Qin, J.; Zhang, L.; Liu, Y.; Sappington, T.W.; Cheng, Y.; Luo, L.; Jiang, X. Population Projection and Development of the *Mythimna loreyi* (Lepidoptera: Noctuidae) as Affected by Temperature: Application of an Age-Stage, Two-Sex Life Table. *J. Econ. Entomol.* **2017**, *110*, 1583–1591. [CrossRef] [PubMed]

12. Qin, J.; Liu, Y.; Zhang, L.; Cheng, Y.; Sappington, T.W.; Jiang, X. Effects of Moth Age and Rearing Temperature on the Flight Performance of the Loreyi Leafworm, *Mythimna loreyi* (Lepidoptera: Noctuidae), in Tethered and Free Flight. *J. Econ. Entomol.* **2018**, *111*, 1243–1248. [CrossRef] [PubMed]

13. El-Far, M.; Li, Y.; Fediere, G.; Abol-Ela, S.; Tijssen, P. Lack of infection of vertebrate cells by the densovirus from the maize worm *Mythimna loreyi* (MIDNV). *Virus Res.* **2014**, *99*, 17–24. [CrossRef] [PubMed]
14. Fediere, G.; El-Far, M.; Li, Y.; Bergoin, M.; Tijssen, P. Expression strategy of densonucleosis virus from *Mythimna loreyi*. *Virology* **2004**, *320*, 181–189. [CrossRef]

15. Hirai, K. The influence of rearing temperature and density on the development of two Leucania species, *M. loreyi* dup. and *L. separata* walker (Lepidoptera: Noctuidae). *Appl. Ent. Zool.* **1975**, *10*, 234–237. [CrossRef]

16. Guo, S.J.; Li, S.M.; Ma, L.P.; Zhuo, X.N. Research about biological characteristics and damage laws of *Leucania loreyi*. *J. Henan Agric. Sci.* **2003**, *9*, 37–39.

17. Jindal, V. DNA barcode reveals occurrence of *Mythimna loreyi* (Duponchel) in Punjab, India. *Indian J. Biotechnol.* **2019**, *18*, 81–84.

18. Blaser, S.; Diem, H.; von Felten, A.; Gueuning, M.; Andreou, M.; Boonham, N.; Tomlinson, J.; Müller, P.; Utzinger, J.; Frey, J.E.; et al. From laboratory to point of entry: Development and implementation of a loop-mediated isothermal amplification (LAMP)-based genetic identification system to prevent introduction of quarantine insect species. *Pest Manag. Sci.* **2018**, *74*, 1504–1512. [CrossRef]

19. Hsieh, C.H.; Wang, H.Y.; Chen, Y.F.; Ko, C.C. Loop-mediated isothermal amplification for rapid identification of biotypes B and Q of the globally invasive pest *Bemisia tabaci*, and studying population dynamics. *Pest. Manag. Sci.* **2012**, *68*, 1206–1213. [CrossRef]

20. Kim, Y.H.; Hur, J.H.; Lee, G.S.; Choi, M.Y.; Koh, Y.H. Rapid and highly accurate detection of *Drosophila suzukii*, spotted wing Drosophila (Diptera: Drosophilidae) by loop-mediated isothermal amplification assays. *J. Asia Pac. Entomol.* **2016**, *19*, 1211–1216. [CrossRef]

21. Notomi, T.; Mori, Y.; Tomita, N.; Kanda, H. Loop-mediated isothermal amplification (LAMP): Principle, features, and future prospects. *J. Microbiol.* **2015**, *53*, 1–5. [CrossRef] [PubMed]

22. Nagamine, K.; Hase, T.; Notomi, T. Accelerated reaction by loop-mediated isothermal amplification using loop primers. *Mol. Cell. Probes* **2002**, *16*, 223–229. [CrossRef] [PubMed]

23. Folmer, O.; Black, M.; Hoeh, W.; Lutz, R.; Vrijenhoek, R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. *Mol. Mar. Biol. Biotechnol.* **1994**, *3*, 294–299. [PubMed]

24. Schwartz, S.; Kent, W.J.; Smit, A.; Zhang, Z.; Baertsch, R.; Hardison, R.C.; Haussler, D.; Miller, W. Human-mouse alignments with BLASTZ. *Genome Res.* **2003**, *13*, 103–107. [CrossRef]

25. Tillich, M.; Lehwark, P.; Pellizzer, T.; Ulbricht-Jones, E.S.; Fischer, A.; Bock, R.; Greiner, S. GeSeq—versatile and accurate annotation of organelle genomes. *Nucleic Acids Res.* **2017**, *45*, W6–W11. [CrossRef]

26. Kumar, S.; Stecher, G.; Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for Bigger Datasets. *Mol. Biol. Evol.* **2016**, *33*, 1870–1874. [CrossRef]

27. Sanderson, M.J.; Wojciechowski, M.F. Improved bootstrap confidence limits in large-scale phylogenies, with an example from Neo-Astragalus (Leguminosae). *Syst. Biol.* **2000**, *49*, 671–685. [CrossRef]

28. Frazer, K.A.; Pachter, L.; Poliakov, A.; Rubin, E.M.; Dubchak, I. VISTA: Computational tools for comparative genomics. *Nucleic Acids Res.* **2004**, *32*, W273–W279. [CrossRef]

29. Mayor, C.; Brudno, M.; Schwartz, J.R.; Poliakov, A.; Rubin, E.M.; Frazer, K.A.; Patcher, L.S.; Dubchak, I. VISTA: Visualizing global DNA sequence alignments of arbitrary length. *Bioinformatics* **2000**, *16*, 1046–1047. [CrossRef]

30. Kim, J.; Kwon, M.; Park, K.J.; Maharjanm, R. Monitoring of four major lepidopteran pests in Korean cornfields and management of *Helicoverpa armigera*. *Entomol. Res.* **2018**, *48*, 308–316. [CrossRef]

31. CABI Invasive Species Compendium. Available online: https://www.cabiorg/isc/search/index (accessed on 14 August 2020).

32. Erasmus, A.; Van Rensburg, J.B.J.; Van den Berg, J. Effects of Bt maize on *Agrotis segetum* (Lepidoptera: Noctuidae): A pest of maize seedlings. *Environ. Entomol.* **2010**, *39*, 702–706. [CrossRef] [PubMed]

33. Kim, J.; Nam, H.Y.; Kwon, M.; Kim, H.; Yi, H.J.; Hänniger, S.; Unbehend, M.; Heckel, D.G. Development of a simple and accurate molecular tool for *Spodoptera frugiperda* species identification using LAMP. *bioRxiv* **2020**. [CrossRef]

34. Li, X.J.; Wu, M.-F.; Ma, J.; Gao, B.-Y.; Wu, Q.-L.; Chen, A.-D.; Liu, J.; Jiang, Y.-Y.; Zhai, B.-P.; Early, R.; et al. Prediction of migratory routes of the invasive fall armyworm in eastern China using a trajectory analytical approach. *Pest Manag. Sci.* **2020**, *76*, 454–463. [CrossRef] [PubMed]

35. Goergen, G.; Kumar, P.L.; Sankung, S.B.; Togola, A.; Tamo, M. First Report of Outbreaks of the Fall Armyworm *Spodoptera frugiperda* (J E Smith) (Lepidoptera, Noctuidae), a New Alien Invasive Pest in West and Central Africa. *PLoS ONE* **2010**, *11*, e0165632. [CrossRef]
36. Mori, Y.; Notomi, T. Loop-mediated isothermal amplification (LAMP): A rapid, accurate, and cost-effective diagnostic method for infectious diseases. *J. Infect. Chemother.* **2009**, *15*, 62–69. [CrossRef]

37. Zhang, X.Z.; Lowe, S.B.; Gooding, J.J. Brief review of monitoring methods for loop-mediated isothermal amplification (LAMP). *Biosens. Bioelectron.* **2014**, *61*, 491–499. [CrossRef]

38. Lee, P.L. DNA amplification in the field: Move over PCR, here comes LAMP. *Mol. Ecol. Resour.* **2017**, *17*, 138–141. [CrossRef]

39. Choi, B.H.; Hur, J.H.; Heckel, D.G.; Kim, J.; Koh, Y.H. Development of a highly accurate and sensitive diagnostic tool for pyrethroid-resistant chimeric P450 CYP337B3 of *Helicoverpa armigera* using loop-mediated isothermal amplification. *Arch. Insect Biochem. Physiol.* **2018**, *99*, e21504. [CrossRef]

40. Parida, M.; Sannarangaiah, S.; Dash, P.K.; Rao, P.V.L.; Morita, K. Loop mediated isothermal amplification (LAMP): A new generation of innovative gene amplification technique; perspectives in clinical diagnosis of infectious diseases. *Rev. Med. Virol.* **2008**, *18*, 407–421. [CrossRef]

41. Kim, J.; Cha, D.J.; Kwon, M.; Maharjan, R. Potato virus Y (PVY) detection in a single aphid by one-step RT-PCR with boiling technique. *Entomol. Res.* **2016**, *46*, 278–285. [CrossRef]

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_44.png" /></figure>
<figure><img src="image_45.png" /></figure>

Article

# Screening of *Helicoverpa armigera* Mobilome Revealed Transposable Element Insertions in Insecticide Resistance Genes

Khouloud KLAI<fnref n="1" />,<fnref n="2" />, Benoît CHÉNAIS<fnref n="2" />, Marwa ZIDI<fnref n="1" /> ID, Salma DJEBBI<fnref n="1" />, Aurore CARUSO<fnref n="2" />, Françoise DENIS<fnref n="2" />, Johann CONFAIS<fnref n="3" />,<fnref n="4" /> ID, Myriam BADAWI<fnref n="2" /> ID, Nathalie CASSE<fnref n="2" />,* and Maha MEZGHANI KHEMAKHEM<fnref n="1" />,* ID

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

check for
updates

Received: 29 October 2020; Accepted: 8 December 2020; Published: 11 December 2020

**Simple Summary:** Transposable elements (TEs) are mobile DNA sequences that can copy themselves within a host genome. TE-mediated changes in regulation can lead to massive and rapid changes in expression, responses that are potentially highly adaptive when an organism is faced with a mortality agent in the environment, such as an insecticide. *Helicoverpa armigera* shows a high number of reported cases of insecticide resistance worldwide, having evolved resistance against pyrethroids, organophosphates, carbamates, organochlorines, and recently to macrocyclic lactone spinosad and several *Bacillus thuringiensis* toxins. In the present study, we conducted a TE annotation using combined approaches, and the results revealed a total of 8521 TEs, representing 236,132 copies, covering 12.86% of the *H. armigera* genome. In addition, we underlined TE insertions in defensome genes and we successfully identified nine TE insertions belonging to the RTE, R2, CACTA, Mariner and hAT superfamilies.

**Abstract:** The cotton bollworm *Helicoverpa armigera* Hübner (Lepidoptera: Noctuidae) is an important pest of many crops that has developed resistance to almost all groups of insecticides used for its management. Insecticide resistance was often related to Transposable Element (TE) insertions near specific genes. In the present study, we deeply retrieve and annotate TEs in the *H. armigera* genome using the Pipeline to Retrieve and Annotate Transposable Elements, PiRATE. The results have shown that the TE library consists of 8521 sequences representing 236,132 TE copies, including 3133 Full-Length Copies (FLC), covering 12.86% of the *H. armigera* genome. These TEs were classified as 46.71% Class I and 53.29% Class II elements. Among Class I elements, Short and Long Interspersed Nuclear Elements (SINEs and LINEs) are the main families, representing 21.13% and 19.49% of the total TEs, respectively. Long Terminal Repeat (LTR) and Dictyostelium transposable element (DIRS) are less represented, with 5.55% and 0.53%, respectively. Class II elements are mainly Miniature Inverted Transposable Elements (MITEs) (49.11%), then Terminal Inverted Repeats (TIRs) (4.09%). Superfamilies of Class II elements, i.e., Transib, P elements, CACTA, Mutator, PIF-harbinger, Helitron, Maverick, Crypton and Merlin, were less represented, accounting for only 1.96% of total TEs. In addition, we highlighted TE insertions in insecticide resistance genes and we successfully identified nine TE insertions belonging to RTE, R2, CACTA, Mariner and hAT superfamilies. These insertions are hosted in genes encoding cytochrome P450 (CyP450), glutathione S-transferase (GST), and ATP-binding cassette

<footnote n="1">¹ Laboratory of Biochemistry and Biotechnology (LR01ES05), Faculty of Sciences of Tunis, University of Tunis El Manar, Tunis 1068, Tunisia; khouloud.klai@fst.utm.tn (K.K.); marwa.zidi@fst.utm.tn (M.Z.); salma.djebbi@fst.utm.tn (S.D.)</footnote>

<footnote n="2">² EA2160 Mer Molécules Santé, Le Mans Université, 72085 Le Mans, France; Benoit.Chenais@univ-lemans.fr (B.C.); Aurore.Caruso@univ-lemans.fr (A.C.); fdenis@univ-lemans.fr (F.D.); Myriam.Badawi@univ-lemans.fr (M.B.)</footnote>

<footnote n="3">³ URGI, INRAE, Université Paris-Saclay, 78026 Versailles, France; johann.confais@inrae.fr</footnote>

<footnote n="4">⁴ BioinfOmics, Plant Bioinformatics Facility, INRAE, Université Paris-Saclay, 78026 Versailles, France</footnote>

<footnote>* Correspondence: Nathalie.Casse@univ-lemans.fr (N.C.); maha.mezghani@fst.utm.tn (M.M.K.)</footnote>
(ABC) transporter belonging to the G and C1 family members. These insertions could therefore be involved in insecticide resistance observed in this pest.

**Keywords:** *Helicoverpa armigera*; transposable elements; insertion sites; insecticide resistance genes

# 1. Introduction

The cotton bollworm, *Helicoverpa armigera* (Hübner) (*Lepidoptera, Noctuidae*), is a serious crop pest having a worldwide distribution [1]. This polyphagous insect causes substantial damages to a wide range of hosts, including cotton, maize, sorghum, and tomato [2]. The biological and ecological traits of *H. armigera*, such as high reproduction rate, polyphagy, high mobility and facultative diapause, make it difficult to control [3]. Management of *H. armigera* attacks rely heavily on the use of chemicals [4]. However, this practice is harmful to the environment and has caused a rapid buildup of insecticide resistance in *H. armigera* populations [5].

Insecticide resistance in *H. armigera* is widespread and has evolved against most of commonly used insecticides [6]. To survive, pests have developed various mechanisms to resist against toxic compounds. These mechanisms include point mutations resulting in target-site resistance such as knockdown resistance (kdr), acetylcholinesterase (Ace-1) and receptor sub-unit termed (RDL) mutations, and also metabolic resistance with involvement of several detoxification enzymes [7].

Metabolic detoxification of toxins is the primary strategy occurring in three phases, each with its own set of enzymes or transporters. Cytochrome P450 monooxygenases (P450s) and carboxylesterases (CarE) carry out phase I, glutathione S-transferases (GSTs) and UDP-glycosyltransferases (UGTs) are phase II enzymes, and ATP-binding cassette transporters (ABC) ensure phase III [8–10]. The understanding of resistance mechanisms remains a challenge that next generation sequencing technologies and the increasing number of sequenced genomes can help to address [11].

Transposable Elements (TEs) are ubiquitous components of eukaryotic genomes that are strongly regulated and inactivated by mutations, which keep transposition events relatively rare [12,13]. However, because of their ability to replicate, TEs may accumulate in host genomes and generate abundant sites for chromosomal rearrangements, which may have deleterious or beneficial consequences [14]. In addition, TEs can provide a selective advantage through their insertion sites, which can enhance or repress gene expression or can be domesticated as new host gene [14–16]. Thus, TEs are an important source of variability for the genomes of their hosts and are therefore key to understanding their evolution. Indeed, TEs may be involved in the genetic adaptation of organisms such as insects to stressful environments, among which is the acquisition of insecticide resistance [17].

Several studies have shown that insecticide resistance can be associated with TE insertions in specific genes. For example, dichlorodiphenyltrichloroethane (DTT) resistance in *Drosophila melanogaster* was correlated with the insertion of a Long Terminal Repeat (LTR)-gypsy retroelement into the 5' region of the cytochrome P450 gene [18,19]. In *Helicoverpa zea*, several TE insertions in regulatory regions, exons and introns of cytochrome P450 genes were related with pyrethroid resistance and xenobiotic metabolism [20].

TEs are classified into two major classes depending on the transposition intermediate. Class I, or retroelements, replicate and transpose via an RNA intermediate; while Class II elements, or DNA transposons, are mobilized via a DNA intermediate [21]. According to the classification of Wicker et al. [22], each class is subdivided into orders and superfamilies. Class I elements are further subdivided into long terminal repeat (LTR) retrotransposons and non-LTR retrotransposons. The non-LTR retrotransposons include the long interspersed nuclear elements (LINEs) and the short interspersed nuclear elements (SINEs) as well as the Penelope-Like Elements (PLEs). Class II elements are subdivided into two subclasses. Subclass 1 includes the terminal inverted repeat (TIR) transposons and Crypton-like elements, which cleave both DNA intermediate strands, while subclass 2 elements including Mavericks
and Helitrons with a single-strand DNA intermediate have a replicative mode of transposition [22–24]. Class II TIRs transposons also include Miniature Inverted-repeat Transposable Elements (MITEs), which are short (~100 to 800 bp) non-autonomous truncated versions of autonomous transposable elements. MITEs possess conserved terminal inverted repeats (TIRs ≥ 10 bp) and a target site duplication (TSDs = 2~10 bp) [25,26].

Annotation of TEs is a challenging task because of their diversity, their repetitive nature and the complexity of their structures, and numerous tools have been designed to identify TEs [27]. In this study, we used the Pipeline to Retrieve and Annotate Transposable Elements (PiRATE) [28] to annotate the mobilome of *H. armigera* and pinpoint TEs inserted in defensome genes.

## 2. Material and Methods

### 2.1. Mobilome Annotation

The *H. armigera* genome available in GenBank-NCBI (BioProject PRJNA378437) is 337,087 Mb. This genome is assembled in 24,552 contigs and 998 scaffolds corresponding to 23.5 kb and 1000 kb N50 length respectively [29].

*H. armigera* genome assembly and the corresponding raw Illumina data were both submitted to the PiRATE pipeline to search for TEs following three steps [28]. In the first step, putative TE sequences were detected using four approaches. The first approach is a similarity-based detection of TEs using RepeatMasker [30] and TE-HMMER [31]. The second approach is a structure-based detection, using LTRharvest [32], MGEScan-nonLTR [33] and SINE-Finder [34]. The repetitiveness-based detection is the third approach using TE de novo [35] and Repeat Scout [36]. The last approach is a de novo approach using the dnaPipeTE tool [37].

After TE detection, a second step was performed to eliminate redundant sequences and classify the remaining sequences using the PASTEC tool following the Wicker's 80-80-80 rules corresponding to sequences longer than 80 bp, sharing more than 80% sequence identity and over 80% of their length [22,38]. Two libraries were generated: a "total TEs library" containing the potentially autonomous TEs and the non-autonomous TEs, and a "repeated elements library" containing the uncategorized repeated sequences and the non-TE sequences. Subsequently, two runs of TEannot [35] were performed for each library to generate the final libraries of total TEs and total repeats. To refine the annotation of TE copies in the whole genome, we used the TEannot pipeline from the REPET package v3.0 with TEs sequences of PiRATE step 1 that align at least with one Full-Length Copy (FLC) on the genome assembly [35,39,40].

Finally, a manual curation was released for all annotated TEs to find corresponding families. This analysis was performed by nucleotide Basic Local Alignment Search Tool (BLAST) against Repbase (48,225 TE sequences) and Dfam (6959 sequences) databases using a threshold value of 80% (Figure 1).

To identify putative MITE sequences, the *H. armigera* genome assembly was submitted to the MITE Tracker tool [41] (Figure 1). This tool searches for putative inverted repeat sequences ranging from 50 to 800 bp. Subsequently, putative MITEs were aligned and clustered into families by Vsearch [42] based on target sites duplication (TSD) and Terminal Inverted Repeat (TIR) sequences.

### 2.2. Search for TE Insertions in Defensome Genes

Annotated TEs from the *H. armigera* genome have been extended by 50 kb both upstream and downstream of their DNA sequences. The nucleotide BLAST was used to find defensome genes in the extended regions using 80% similarity and 80% query coverage threshold (Figure 1).
<figure><img src="image_47.png" /><figcaption>**Figure 1.** Flowchart for transposable elements annotation and insertion sites identification in the *H. armigera* genome.</figcaption></figure>

# 3. Results

## 3.1. TEs Annotation of the *H. armigera* Genome

The screening of TEs in the *H. armigera* genome using different detection tools led to the identification of 100,184 TE candidates after redundancy elimination. The classification step generated two libraries: the “total TE library” and the “total repeats library” containing 4336 sequences (11.36%) and 5201 sequences (13.63%), respectively. Among the TE library, 3133 sequences were identified as FLC belonging to Class I (2720 FLC), Class II (413 FLC) and 61 sequences were undefined (Supplementary File 1). A total of 70,030 sequences was classified as non-TE (41,411 sequences) and unclassified TEs (28,619 sequences). The results revealed that TE sequences cover 12.86% (43,349,853 bp) of the *H. armigera* genome and most of the TE sequences belong to Class II elements, accounting for 53.29% of the total TE content, while Class I elements account for 46.71%.

Among Class I elements, SINEs and LINEs were the main families, representing 21.13% and 19.49% of the total TEs, respectively. LTR elements were represented with 5.55% and Dictyostelium transposable element (DIRS) with only 0.53%. The Class II elements were represented mainly by MITEs and TIRs with 49.11% and 4.09%, respectively.

To investigate the evolutionary history of TEs in the *H. armigera* genome, we plotted the distribution of identity values between copies and their representative sequences. Distributions of TE classes showed a peak at 80% identity for Class I elements, while, for Class II elements, the distribution was linear with a recent burst at 98% identity (Figure 2).

This analysis revealed that the *H. armigera* genome has undergone a multitude of ancient and recent bursts of different TE superfamilies showing its fluidity. Distributions of TE copies showed three peaks of transposition activity (Figure 3). The first peak is at 65% divergence involving a burst of the DIRS order. We also noted a second TE burst, particularly for LINE, SINE, LTR, TIR, and Helitron orders at 80% identity. In addition, the distribution also showed the appearance of MITE elements at 95% identity, suggesting a recent invasion of the *H. armigera* genome by these TEs (Figure 3).

### 3.1.1. Class I Retrotransposons

The annotation of Class I retrotransposons in the *H. armigera* genome allowed for the identification of 3980 sequences representing 186,645 copies belonging to 10 superfamilies (Table 1).
<figure><img src="image_48.png" /><figcaption>**Figure 2.** Distribution of sequence identity values between TE copies and TE sequences with at least one full-length copy for Class I and Class II elements. The relative frequencies per percentage of identity of Class I and Class II are represented in different colors.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>**Figure 3.** Distribution of sequence identity values between TE copies and TE sequences with at least one FLC. The relative frequencies per percentage of identity of Dictyostelium transposable element (DIRS), Helitron, Long and Short Interspersed Nuclear Elements (LINE and SINE), Long Terminal Repeat (LTR), Miniature Inverted Transposable Element (MITE) and Terminal Inverted Repeat (TIR) orders are represented in different colors. Only the main orders (in terms of copy number) are represented.</figcaption></figure>

From the LTR retrotransposons, 345 TEs were full-length copies (FLC). Gypsy was the most abundant LTR superfamily with 241 sequences and 14,876 copies followed by Bel-Pao (155 sequences) and Copia (77 sequences) with 97 and 67 FLC, respectively. BLAST searches against Repbase and Dfam databases showed that, for all identified LTR sequences, no similarity was found with TEs in databases.
<table><caption>**Table 1.** Summary of the identified and annotated TEs in the *H. armigera* genome.</caption><thead><tr><th>Class</th><th>Order</th><th>Superfamily</th><th>Total Sequences <sup>1</sup></th><th>Full-Length Copies <sup>2</sup></th><th>TE Percentage</th><th>Copy Number <sup>3</sup></th></tr></thead><tbody><tr><td rowspan="12">Class I</td><td rowspan="3">LTR</td><td>Gypsy</td><td>241</td><td>181</td><td></td><td>14,876</td></tr><tr><td>Bel pao</td><td>155</td><td>97</td><td></td><td>1792</td></tr><tr><td>Copia</td><td>77</td><td>67</td><td></td><td>223</td></tr><tr><td>Total LTR</td><td></td><td>473</td><td>345</td><td>5.55%</td><td>16,891</td></tr><tr><td>DIRS</td><td>DIRS</td><td>45</td><td>38</td><td></td><td>149</td></tr><tr><td>Total DIRS</td><td></td><td>45</td><td>38</td><td>0.53%</td><td>149</td></tr><tr><td rowspan="5">LINE</td><td>Jockey</td><td>716</td><td>480</td><td></td><td>44,432</td></tr><tr><td>RTE</td><td>604</td><td>543</td><td></td><td>46,458</td></tr><tr><td>I</td><td>172</td><td>93</td><td></td><td>16,520</td></tr><tr><td>R2</td><td>169</td><td>101</td><td></td><td>10,614</td></tr><tr><td>Total LINE</td><td></td><td>1661</td><td>1217</td><td>19.49%</td><td>118,024</td></tr><tr><td rowspan="3">SINE</td><td>tRNA</td><td>1797</td><td>1120</td><td></td><td>51,581</td></tr><tr><td>5S</td><td>4</td><td>-</td><td></td><td>0</td></tr><tr><td>Total SINE</td><td></td><td>1801</td><td>1120</td><td>21.13%</td><td>51,581</td></tr><tr><td>Total Class I</td><td></td><td></td><td>3980</td><td>2720</td><td>46.71%</td><td>186,645</td></tr><tr><td rowspan="12">Class II<br>Subclass I</td><td rowspan="2">TIR</td><td>hAT</td><td>120</td><td>105</td><td></td><td>6267</td></tr><tr><td>Mariner</td><td>106</td><td>74</td><td></td><td>2284</td></tr><tr><td>Piggybac</td><td>41</td><td>32</td><td></td><td>6670</td></tr><tr><td>Transib</td><td>40</td><td>32</td><td></td><td>188</td></tr><tr><td rowspan="5">P</td><td></td><td>15</td><td>12</td><td></td><td>6103</td></tr><tr><td>CACTA</td><td>14</td><td>11</td><td></td><td>201</td></tr><tr><td>Mutator</td><td>6</td><td>4</td><td></td><td>477</td></tr><tr><td>PIF_harbinger</td><td>6</td><td>5</td><td></td><td>90</td></tr><tr><td>Merlin</td><td>1</td><td>1</td><td></td><td>30</td></tr><tr><td>Total TIR</td><td></td><td>349</td><td>276</td><td>4.09%</td><td>22,310</td></tr><tr><td rowspan="4">Class II<br>Subclass II</td><td>Crypton</td><td>Crypton</td><td>1</td><td>1</td><td></td><td>3</td></tr><tr><td>Helitron</td><td>Helitron</td><td>4</td><td>4</td><td></td><td>595</td></tr><tr><td>Maverick</td><td>Maverick</td><td>2</td><td>2</td><td></td><td>14</td></tr><tr><td>MITEs</td><td>MITEs</td><td>4185</td><td>130</td><td>49.11%</td><td>26,565</td></tr><tr><td>Total classII</td><td></td><td></td><td>4541</td><td>413</td><td>53.29%</td><td>49,487</td></tr><tr><td>Total TEs</td><td></td><td></td><td>8521</td><td>3133</td><td>100</td><td>236,132</td></tr></tbody></table>

Regarding LINEs retrotransposons, the most abundant elements belong to the Jockey and RTE superfamilies, with 716 and 604 sequences, respectively, representing 44,432 and 46,458 copies, respectively. The RTE superfamily contains the highest number of FLCs among all annotated TEs (Table 1). According to the BLAST searches, no similarity was found for Jockey elements in databases while 17 RTE sequences showed similarity ranging from 86% to 100% with RTE-1_Avan, RTE-2_Hmel_C, RTE-3_DPl, RTE-4_DPl, RTE-5_DPl and Proto2-1_BM families (Table S1).

Concerning SINEs retrotransposons, the tRNA-derived SINE superfamily corresponds to 45.15% of all Class I TEs with 1797 sequences representing 51,581 copies and 1120 FLCs. Among these sequences, 578 SINEs belong to the HaSE1 family while 85 sequences fell into the HaSE3 family with a similarity ranging from 80% to 100% (Table S1).

<footnote n="1">¹ Representative sequence identified with PiRATE Step1 with an identity ≤ 80% ² TEs sequences of PiRATE step 1 that align at least with one Full-length Copy (FLC) on the genome assembly ³ Copies annotated with TEannot REPET package v3.0. Number of undefined transposable elements are not shown in the table.</footnote>
### 3.1.2. Class II Transposons

Our results revealed that DNA transposons in the *H. armigera* genome are represented by a total of 4541 sequences representing 49,487 copies belonging to 13 superfamilies (Table 1).

Regarding TIR elements, the hAT superfamily was the most abundant, including 120 sequences, among which 105 elements were FLC. BLAST searches against Repbase and Dfam databases revealed high similarity of three hAT elements with the hAT-1_DAN family (Table S1).

The Tc1/mariner superfamily was represented by 106 elements (2284 copies), among which 74 FLC. Research using BLAST has shown that 12 of these elements are distributed among five families—ANM4, nMar-2_Avan, nMar-18_Hmel, Mariner-1_AMel and Mariner-3_BM—exhibiting similarity ranging from 80.25% to 95.34% (Table S1).

A total of 41 sequences (6670 copies) belonging to the PiggyBac superfamily was identified in the *H. armigera* genome. Further analysis revealed that only four among the identified sequences showed similarity ranging from 81% to 99% with already identified PiggyBac transposons, npiggyBac-8 and piggyBac-2.

The following superfamilies of Class II elements, i.e., Transib, P elements, CACTA, Mutator, PIF-harbinger, Helitron, Maverick, Crypt-on and Merlin, were less represented in the *H. armigera* genome corresponding to a total of 7701 copies.

In addition, the MITE tracker allowed for the identification of 4185 putative MITEs corresponding to 43.11% of all annotated TEs. The analysis of terminal TIRs and TSD sequences allowed for the classification of 3570 MITE sequences (26,565 copies) into seven superfamilies (Table 2). The Tc1\mariner and PIF-Harbinger superfamilies were the most represented, with 1817 and 1368 MITEs, respectively, followed by the CACTA superfamily with 250 MITEs, then PiggyBac with 93 elements. The hAT, Transib and Maverick superfamilies were represented by only 20, 16 and six MITEs, respectively.

<table><caption>**Table 2.** Distribution of MITEs identified in the *H. armigera* genome.</caption><thead><tr><th>Superfamily</th><th>TSD</th><th>Number of MITEs</th><th>MITE Length (bp)</th><th>TIR Length (bp)</th></tr></thead><tbody><tr><td>Tc1/mariner</td><td>TA</td><td>1817</td><td>50–360</td><td>10–21</td></tr><tr><td>PIF-Harbinger</td><td>TWA</td><td>1368</td><td>55–685</td><td>15–32</td></tr><tr><td>CACTA</td><td>2–3 bp</td><td>250</td><td>78–775</td><td>10–26</td></tr><tr><td>Piggybac</td><td>TTAA</td><td>93</td><td>50–800</td><td>15–31</td></tr><tr><td>hAT</td><td>T****A</td><td>20</td><td>56–260</td><td>17–29</td></tr><tr><td>Transib</td><td>C***G</td><td>16</td><td>83–386</td><td>13–27</td></tr><tr><td>Maverick</td><td>6 bp</td><td>6</td><td>50–800</td><td>10–24</td></tr><tr><td>Other</td><td>-</td><td>615</td><td>50–800</td><td>11–35</td></tr></tbody></table>

### 3.2. TE Insertions Scanning in Defensome Genes

Nucleotide BLAST searches for defensome genes in the regions framing identified TE sequences led to the identification of nine TE insertion sites in seven genes encoding for detoxifying enzymes (Table 3). The involved TEs are members of RTE, R2, CACTA, DIRS, Mariner, and hAT superfamilies (Table 3). Further analysis of TE insertions has shown that five TEs were inserted in four cytochrome P450 genes, an element was retrieved in a GST gene, two were hosted by ABC-G transporter gene and one was inserted in an ABC-C1 transporter gene.
<table><caption>**Table 3.** TE insertions in genes encoding for detoxifying enzymes in *H. armigera*.</caption><thead><tr><th>Gene Family</th><th>Gene Name</th><th>Gene Length (bp)</th><th>Inserted Element Name</th><th>TE Length (bp)</th><th>Insertion Position</th></tr></thead><tbody><tr><td rowspan="4">Cytochrome P450 (CYP450)</td><td>4g15-like (LOC110375617)</td><td>14,974</td><td>LINE_R2_526</td><td>4083</td><td>6483–10,566</td></tr><tr><td rowspan="2">4C1-like (LOC113006340)</td><td rowspan="2">39,707</td><td>LINE_RTE17512</td><td>620</td><td>5655–6274</td></tr><tr><td>LINE_RTE84107</td><td>740</td><td>7021–7763</td></tr><tr><td>4C1-like (LOC110381376)</td><td>8332</td><td>LTR_CACTA3565</td><td>1179</td><td>6899–8066</td></tr><tr><td>6k1-like (LOC110381042)</td><td>15,958</td><td>LINE_jockey2199</td><td>1629</td><td>9431–10,438</td></tr><tr><td>Glutathione S-transferase (GST)</td><td>GST 1-like (LOC110371343)</td><td>4001</td><td>TIR_Mariner2770</td><td>1304</td><td>1793–3097</td></tr><tr><td rowspan="2">ATP binding cassette (ABC) transporter</td><td rowspan="2">ABC-G member 20 (LOC110376033)</td><td rowspan="2">96,146</td><td>TIR_Mariner419</td><td>1085</td><td>91,998–93,071</td></tr><tr><td>TIR_hAT2824</td><td>951</td><td>1486–2436</td></tr><tr><td>ABC-C1 (Multidrug resistance protein homolog 49-like MRP1) (LOC110377844)</td><td>48,302</td><td>LINE_RTE63004</td><td>894</td><td>38,590–39,483</td></tr></tbody></table>

Six of the inserted TEs have intronic insertion sites and one TE insertion occurred in the first exon of the ABCG transporter member 20 gene (Figure 4 and Figures S1–S7). It should be noted that, for the cytochrome P450 4C1-like (LOC11300634) gene, harboring two LINE insertions, no exon or intron information was retrieved in GenBank.

<figure><img src="image_50.png" /><figcaption>**Figure 4.** Schematic representation of TE insertion sites in genes encoding for detoxifying enzymes in *H. armigera*.</figcaption></figure>
## 4. Discussion

The present study identified TEs present in the genome of *H. armigera* and searched for their occurrence among the defensome genes of this pest. Characterizing TEs is an important task for non-model organisms and several TE annotation tools have recently been developed for TE characterization in these organisms [43]. In this work, we used the PiRATE pipeline [28], which combines different TE identification tools to detect, classify and annotate TEs into known superfamilies.

The results revealed a total TE content of 4336 sequences, covering 12.86% of the *H. armigera* assembly which is much higher than the previous data from Pearce et al. [29], in which Repeat Masker and Repeat Modeller tools allowed for the identification of a TE content representing only 0.88% of the genome. This confirms the interest of using a pipeline like PiRATE to increase the detection of TEs in a genome.

In other lepidopteran genomes, such as *Bombyx mori* and *Spodoptera frugiperda*, about 53% of TEs were previously identified as Class I and half of these belong to the *LINE* order [44]. Consistent with these results, 46.71% of the TEs identified in the *H. armigera* genome were Class I elements and the two most abundant were the *SINE* and *LINE* orders. These results suggest that using a combined approach is more specific for short TEs than the use of a single tool.

Some TEs were previously identified in *H. armigera* by molecular biology methods. In 2008, Sun et al. identified two *PiggyBac* elements (HaPLE1 and HaPLE2) in the cotton bollworm genome [45]. A few years later, multiple copies of two distinct mariner elements, *Hamar1* and *Hamar2*, were isolated by Wang et al. [46]. All these sequences were retrieved in the current study with a minimum of 88% identity.

To pinpoint TE insertions in defensome genes, the upstream and downstream regions of TE sequences were scanned and nine insertions have been successfully identified. Among them, five were retrieved in cytochrome P450 genes introns. Cytochrome P450 are among the three main groups of detoxification enzymes used by insects that play crucial roles in the detoxification of most pesticides [47]. In *H. armigera*, two *LINEs* (RTE and R2) and a *CACTA* element were inserted in the P450 4C1-like gene, the expression of which is responsible for the insecticide detoxification [48]. In *Aphis gossypii*, another LINE, the *Jockey* element was also identified in cytP450 6K1-like gene intron, putting it as the strongest candidate conferring resistance to thiamethoxam [49]. Two other *LINE* elements were inserted into the P450 4g15-like gene, but no information for exon or intron regions was retrieved in GenBank for this gene, suggesting that the *H. armigera* genome needs more annotation. However, the expression of this gene in *Aphis gossypii* was correlated with imidacloprid resistance [50].

Arthropod ABC transporters also play an important role in metabolite detoxification [10]. Three TEs, which belong to *RTE*, *mariner* and *hAT* superfamilies, were inserted in an ABCG transporter and the gene in the *H. armigera* genome. TE insertions into non-coding regions may be less subject to selection and could be influenced by other mechanisms of control. These insertions could be successfully spliced out during mRNA processing and thus may have no obvious effects on the function of the corresponding defensome gene [51]. Alternatively, they could result in exon skipping, alternative splicing, or alterations in expression profiles if the corresponding IR gene introns contain regulatory sequences, as exemplified by the insertion of the *Mu* element into an intron of the *Knotted* locus in Maize [51]. Analyses of the transcript and the expression levels of the identified IR genes as well as toxicity bioassay will be necessary to determine the exact effects of the nine identified TE insertions on the expression and function of the identified IR genes, as well as on the fitness of this species in the presence of a stressful environment such as insecticides.

## 5. Conclusions

Genome-wide TE annotation has rarely been performed in *H. armigera*. This study opens the way to new searches about the role of TEs in the genome evolution of *H. armigera* and their contribution to the pest adaptation such as insecticide resistance. In the present study, we conducted an accurate TE annotation, and the results reveal a total of 8521 TEs covering 12.86% (43349853 bp) of the *H. armigera*
genome. The annotation of TEs was crossed with insertion sites search in defensome genes. Nine TEs belonging to the RTE, R2, LTR and TIR superfamilies were found to be inserted in CYP450, GST and ABC transporter genes and their insertion sites were mostly in intronic regions except for a hAT element inserted in the exon region. These results present for the first-time evidence that TEs are present in IR genes of *H. armigera*. However, further studies are necessary to elucidate the functional relationship of TEs with IR genes in *H. armigera*.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-4450/11/12/879/s1.
File 1. Full-length TEs annotated in the *Helicoverpa armigera* genome, Table S1: BLAST results of annotated TEs in *H. armigera* against reference TEs in databases (Repbase and Dfam), Figure S1. Alignment of LINE_R2526 element with Cytochrome P450 4g15-like gene, Figure S2. Alignment of LINE_RTE84107 and LINE_RTE17512 elements with Cytochrome P450 4C-like gene, Figure S3. Alignment of TIR_CACTA3565 element with Cytochrome P450 4C-like gene, Figure S4. Alignment of LINE_Jockey2199 element with Cytochrome P450 6k1-like gene, Figure S5. Alignment of TIR.mariner2770 element with glutathione S-transferase 1-like gene, Figure S6. Alignment of LINE_RTE63004 element with ABC-C1 homolog 49-like gene, Figure S7. Alignment of TIR.mariner419 and TIR_hAT824 elements with ABCG member 20 gene.

**Author Contributions:** Conceptualization, K.K., M.M.K., A.C., B.C. and N.C.; methodology, K.K., M.Z., M.B., J.C. and F.D.; investigation, data curation, K.K.; writing—original draft preparation, K.K., M.M.K., B.C. and N.C.; writing—review and editing, M.M.K., B.C., N.C., A.C., S.D., F.D.; supervision, M.M.K., B.C. and N.C.; project administration, M.M.K. and N.C.; funding acquisition, M.M.K. and N.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by PHC-Utique "Hubert Curien Franco-Tunisian partnership" (41703ZJ) and the Tunisian Ministry of Higher Education and Scientific Research and the University of Tunis El Manar.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Burgio, G.; Ravaglia, F.; Maini, S.; Bazzocchi, G.G.; Masetti, A.; Lanzoni, A. Mating Disruption of *Helicoverpa armigera* (Lepidoptera: Noctuidae) on Processing Tomato: First Applications in Northern Italy. *Insects* **2020**, *11*, 206. [CrossRef] [PubMed]

2. Tembrock, L.R.; Timm, A.E.; Zink, F.A.; Gilligan, T.M. Phylogeography of the Recent Expansion of *Helicoverpa armigera* (Lepidoptera: Noctuidae) in South America and the Caribbean Basin. *Ann. Entomol. Soc. Am.* **2019**, *112*, 388–401. [CrossRef]

3. Sarate, P.J.; Tamhane, V.A.; Kotkar, H.M.; Ratnakaran, N.; Susan, N.; Gupta, V.S.; Giri, A.P. Developmental and digestive flexibilities in the midgut of a polyphagous pest, the cotton bollworm, *Helicoverpa armigera*. *J. Insect Sci.* **2012**, *12*, 42. [CrossRef] [PubMed]

4. Pereira, F.P.; Reigada, C.; Diniz, A.J.F.; Parra, J.R.P. Potential of Two Trichogrammatidae species for *Helicoverpa armigera* control. *Neotrop. Entomol.* **2019**, *48*, 966–973. [CrossRef]

5. Ansari, M.S.; Moraiet, M.A.; Ahmad, S. Insecticides: Impact on the Environment and Human Health. In *Environmental Deterioration and Human Health: Natural and Anthropogenic Determinants*; Malik, A., Grohmann, E., Akhtar, R., Eds.; Springer: Dordrecht, The Netherlands, 2014; pp. 99–123. [CrossRef]

6. Srinivas, R.; Udikeri, S.S.; Jayalakshmi, S.K.; Sreeramulu, K. Identification of factors responsible for insecticide resistance in *Helicoverpa armigera*. *Comp. Biochem. Physiol.* **2004**, *137*, 261–269. [CrossRef]

7. Sarwar, M.; Salman, M. Insecticides Resistance in Insect Pests or Vectors and Development of Novel Strategies to Combat Its Evolution. *IJBBE* **2015**, *1*, 344–351.

8. Sampath, P.; Murukarthick, J.; Izzah, N.K.; Lee, J.; Choi, H.I.; Shirasawa, K.; Choi, B.S.; Liu, S.; Nou, I.S.; Yang, T.J. Genome-Wide Comparative Analysis of 20 Miniature Inverted-Repeat Transposable Element Families in *Brassica rapa* and *B. oleracea*. *PLoS ONE* **2014**, *9*, e94499. [CrossRef]

9. De Marco, L.; Sassera, D.; Epis, S.; Mastrantonio, V.; Ferrari, M.; Ricci, I.; Comandatore, F.; Bandi, C.; Porretta, D.; Urbanelli, S. The choreography of the chemical defensome response to insecticide stress: Insights into the *Anopheles stephensi* transcriptome using RNA-Seq. *Sci. Rep.* **2017**, *7*, [CrossRef]

10. Dermauw, W.; van Leeuwen, T. The ABC gene family in arthropods: Comparative genomics and role in insecticide transport and resistance. *Insect Biochem. Mol. Biol.* **2014**, *45*, 89–110. [CrossRef]
11. Su, H.; Gao, Y.; Liu, Y.; Li, X.; Liang, Y.; Dai, X.; Xu, Y.; Zhou, Y.; Wang, H. Comparative transcriptome profiling reveals candidate genes related to insecticide resistance of *Glyphodes pyloalis*. *Bull. Entomol. Res.* **2020**, *110*, 57–67. [CrossRef]

12. Mérel, V.; Boulesteix, M.; Fablet, M.; Vieira, C. Transposable elements in *Drosophila*. *Mob. DNA* **2020**, *11*, 23. [CrossRef] [PubMed]

13. Lavoie, C.A. Transposable Element Content in Non-Model Insect Genomes. Master's Thesis, Mississippi State University, Starkwell, MS, USA, 2020.

14. Chénais, B.; Caruso, A.; Hiard, S.; Casse, N. The impact of transposable elements on eukaryotic genomes: From genome size increase to genetic adaptation to stressful environments. *Gene* **2012**, *509*, 7–15. [CrossRef] [PubMed]

15. Bourgeois, Y.; Boissinot, S. On the Population Dynamics of Junk: A Review on the Population Genomics of Transposable Elements. *Genes* **2019**, *10*, 419. [CrossRef] [PubMed]

16. Steinberg, C.E.W. *Arms Race between Plants and Animals: Biotransformation System in Stress Ecology: Environmental Stress as Ecological Driving Force and Key Player in Evolution*; Springer: New York, NY, USA, 2012; pp. 61–105.

17. Le Goff, G.; Hilliou, F. Resistance evolution in *Drosophila*: The case of CYP6G1. *Pest Manag. Sci.* **2017**, *73*, 493–499. [CrossRef]

18. Catania, F.; Kauer, M.O.; Daborn, P.J.; Yen, J.L.; French-Constant, R.H.; Schlötterer, C. World-wide survey of an Accord insertion and its association with DDT resistance in *Drosophila melanogaster*. *Mol. Ecol.* **2004**, *13*, 2491–2504. [CrossRef]

19. Schlenke, T.A.; Begun, D.J. Strong selective sweep associated with a transposon insertion in *Drosophila simulans*. *Proc. Natl. Acad. Sci. USA* **2004**, *101*, 1626–1631. [CrossRef]

20. Chen, S.; Li, X. Transposable elements are enriched within or in close proximity to xenobiotic-metabolizing cytochrome P450 genes. *BCM Evol. Biol.* **2007**, *7*, 46. [CrossRef]

21. Finnegan, D.J. Eukaryotic transposable elements and genome evolution. *Trends Genet.* **1989**, *5*, 103–107. [CrossRef]

22. Wicker, T.; Sabot, F.; Hua-Van, A.; Bennetzen, J.L.; Capy, P.; Chalhoub, B.; Flavell, A.; Leroy, P.; Morgante, M.; Panaud, O.; et al. An unified classification system for eukaryotic transposable elements. *Nat. Rev. Genet.* **2007**, *8*, 973–982. [CrossRef]

23. Kapitonov, V.V.; Jurka, J. Rolling-Circle Transposons in Eukaryotes. *Proc. Natl. Acad. Sci. USA* **2001**, *98*, 8714–8719. [CrossRef]

24. Pritham, E.J.; Putliwala, T.; Feschotte, C. Mavericks, a novel class of giant transposable elements widespread in eukaryotes and related to DNA viruses. *Gene* **2007**, *390*, 3–17. [CrossRef] [PubMed]

25. Kuang, H.; Padmanabhan, C.; Li, F.; Kamei, A.; Bhaskar, P.B.; Ouyang, S.; Jiang, J.; Robin Buell, C.; Baker, B. Identification of miniature inverted-repeat transposable elements (MITEs) and biogenesis of their SiRNAs in the Solanaceae: New functional implications for MITEs. *Genome Res.* **2009**, *19*, 42–56. [CrossRef] [PubMed]

26. Schrader, L.; Schmitz, J. The impact of transposable elements in adaptive evolution. *Mol. Ecol.* **2019**, *28*, 1537–1549. [CrossRef] [PubMed]

27. Yan, H.; Torchiana, F.; Bombarely, A. Guideline for genome transposon annotation derived from evaluation of popular TE identification tools. *Preprints* **2020**, 2020080275. [CrossRef]

28. Berthelier, J.; Casse, N.; Daccord, N.; Jamilloux, V.; Saint-Jean, B.; Carrier, G. A transposable element annotation pipeline and expression analysis reveal potentially active elements in the microalga *Tisochrysis lutea*. *BMC Genom.* **2018**, *19*, 1–14. [CrossRef]

29. Pearce, S.L.; Clarke, D.F.; East, P.D.; Elfekih, S.; Gordon, K.H.J.; Jermiin, L.S.; McGaughran, A.; Oakeshott, J.G.; Papanikolaou, A.; Perera, O.P.; 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–30. [CrossRef]

30. Smit, A.F.A.; Hubley, R.; Green, P. Repeat Masker. Unpublished Data. Current Version: Open-4.09 (Dfam:3.0 only) 1996. Available online: http://www.repeatmasker.org (accessed on 15 April 2019).

31. Eddy, S.R. Multiple alignment using hidden Markov models. *ISMB* **1996**, *3*, 114–120.

32. Ellinghaus, D.; Kurtz, S.; Willhoefft, U. LTR harvest, an efficient and flexible software for de novo detection of LTR retrotransposons. *BMC Bioinform.* **2008**, *9*, 18. [CrossRef]

33. Rho, M.; Tang, H. MGEScan-Non-LTR: Computational identification and classification of autonomous non-LTR retrotransposons in eukaryotic genomes. *Nucleic Acids Res.* **2009**, *37*, e143. [CrossRef]
34. Wenke, T.; Döbel, T.; Sörensen, T.R.; Junghans, H.; Weisshaar, B.; Schmidta, T. Targeted Identification of Short Interspersed Nuclear Element Families Shows Their Widespread Existence and Extreme Heterogeneity in Plant Genomes. *Plant Cell* **2011**, *23*, 3117–3128. [CrossRef]

35. Flutre, T.; Duprat, E.; Feuillet, C.; Quesneville, H. Considering Transposable Element Diversification in de Novo Annotation Approaches. *PLoS ONE* **2011**, *6*, e16526. [CrossRef] [PubMed]

36. Price, A.L.; Jones, N.C.; Pevzner, P.A. De novo identification of repeat families in large genomes. *Bioinformatics* **2005**, *21* (Suppl. 1), i351–i358. [CrossRef] [PubMed]

37. Goubert, C.; Modolo, L.; Vieira, C.; Moro, C.V.; Mavingui, P.; Boulesteix, M. De Novo Assembly and Annotation of the Asian Tiger Mosquito (*Aedes albopictus*) Repeatome with DnaPipeTE from Raw Genomic Reads and Comparative Analysis with the Yellow Fever Mosquito (*Aedes aegypti*). *GBE* **2015**, *7*, 1192–1205. [CrossRef] [PubMed]

38. Hoede, C.; Arnoux, S.; Moisset, M.; Chaumier, T.; Inizan, O.; Jamilloux, V.; Quesneville, H. PASTEC: An Automatic Transposable Element Classification Tool. *PLoS ONE* **2014**, *9*, e91929. [CrossRef]

39. Quesneville, H.; Bergman, C.M.; Andrieu, O.; Autard, D.; Nouaud, D.; Ashburner, M.; Anxolabehere, D. Combined Evidence Annotation of Transposable Elements in Genome Sequences. *PLoS Comput. Biol.* **2005**, *1*, 0166–0175. [CrossRef]

40. Daccord, N.; Celton, J.M.; Linsmith, G.; Becker, C.; Choisne, N.; Schijlen, E.; van de Geest, H.; Bianco, L.; Micheletti, D.; Velasco, R.; et al. High-quality de novo assembly of the apple genome and methylome dynamics of early fruit development. *Nat. Genet.* **2017**, *49*, 1099–1106. [CrossRef]

41. Crescente, J.M.; Zavallo, D.; Helguera, M.; Vanzetti, L.S. MITE Tracker: An accurate approach to identify miniature inverted-repeat transposable elements in large genomes. *BMC Bioinform.* **2018**, *19*, 348. [CrossRef]

42. Rognes, T.; Flouri, T.; Nichols, B.; Quince, C.; Mahé, F. VSEARCH: A versatile open source tool for metagenomics. *PeerJ* **2016**, *4*, 2584. [CrossRef]

43. Goerner-Potvin, P.; Bourque, G. Computational tools to unmask transposable elements. *Nat. Rev. Genet.* **2018**, *19*, 688–704. [CrossRef]

44. D'Alençon, E.; Sezutsu, H.; Legeai, F.; Permal, E.; Bernard-Samain, S.; Gimenez, S.; Gagneur, C.; Cousserans, F.; Shimomura, M.; Brun-Barale, A.; et al. Extensive synteny conservation of holocentric chromosomes in lepidoptera despite high rates of local genome rearrangements. *Proc. Natl. Acad. Sci. USA* **2010**, *107*, 7680–7685. [CrossRef]

45. Sun, Z.C.; Wu, M.; Miller, T.A.; Han, Z.J. PiggyBac-like elements in cotton bollworm, *Helicoverpa armigera* (Hübner). *Insect Mol. Biol.* **2008**, *17*, 9–18. [CrossRef] [PubMed]

46. Wang, J.; Miller, T.A.; Park, Y. Identification of mariner-like elements belonging to the Cecropia subfamily in two closely related *Helicoverpa* species. *Insect Sci.* **2011**, *18*, 619–628. [CrossRef]

47. Zhu, W.; Yu, R.; Wu, H.; Zhang, X.; Liu, Y.; Zhu, K.Y.; Zhang, J.; Ma, E. Identification and characterization of two CYP9A genes associated with pyrethroid detoxification in *Locusta migratoria*. *Pestic. Biochem. Physiol.* **2016**, *132*, 65–71. [CrossRef] [PubMed]

48. Gong, J.; Cheng, T.; Wu, Y.; Yang, X.; Feng, Q.; Mita, K. Genome-wide patterns of copy number variations in *Spodoptera litura*. *Genomics* **2019**, *111*, 1231–1238. [CrossRef]

49. Langfield, K.L. Characterization of Neonicotinoid Resistance in the Cotton Aphid, *Aphis gossypii* from Australian Cotton. Ph.D. Thesis, University of Technology Sydney, Sydney, Australia, April 2017.

50. Kim, J.I.; Kwon, M.; Kim, G.H.; Kim, S.Y.; Lee, S.H. Two mutations in nAChR beta subunit is associated with imidacloprid resistance in the *Aphis gossypii*. *J. Asia Pac. Entomol.* **2015**, *18*, 291–296. [CrossRef]

51. Greene, B.; Wako', R.; Hake, S. Mutator Insertions in an Intron of the Maize Knotted1 Gene Result in Dominant Suppressible Mutations. *Genetics* **1994**, *138*, 1275–1285.

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
<figure><img src="image_52.png" /></figure>

Article

# A Dynamic Energy Budget Approach for the Prediction of Development Times and Variability in *Spodoptera frugiperda* Rearing

Andre Gergs * and Christian U. Baden

Bayer AG, Crop Science Division, Alfred-Nobel Straße 50, 40789 Monheim, Germany; christian.baden@bayer.com

* Correspondence: andre.gergs@bayer.com

**Simple Summary:** The fall armyworm *Spodoptera frugiperda* is a moth that is active during the night. Its larvae cause extensive damage to many crops. Laboratory experiments are conducted to find effective ways to control this pest. One important aspect in research, generally, is that experimental results are reproducible. Reproducibility directly depends on the homogeneity of the test material— the fall armyworm larvae, in our case. The more variable the conditions of the larvae in terms of larval stages or sizes, the more variable and less reliable the research results will be. We used a mathematical model to explore the causes for increased variability in the larval development of the fall armyworm. We found that low air temperatures and poor nutrition increase development times and variability compared to higher air temperature settings and good-quality food. This finding helps researchers to adjust rearing temperatures in a way that allows starting experiments with specific larval stages and low variability on time as planned for their high-quality research.

**Abstract:** A major challenge in insect rearing is the need to provide certain life cycle stages at a given time for the initiation of experimental trials. The timing of delivery, organism quality, and variability directly affect the outcome of such trials. Development times and intraspecific variability are directly linked to the availability of food and to the ambient temperature. Varying temperature regimes is an approach to adapt development times to fulfill experimental needs without impairment of larval quality. However, current practices of temperature setting may lead to increased variability in terms of development times and the frequency of particular life stages at a given point in time. In this study, we analyzed how resource availability and ambient temperature may affect the larval development of the economically important noctuid species *Spodoptera frugiperda* by means of dynamic energy budget modeling. More specifically, we analyzed how rearing practices such as raising of temperatures may affect the variability in larval development. Overall, the presented modeling approach provides a support system for decisions that must be made for the timely delivery of larvae and reduction of variability.

**Keywords:** insect rearing; dynamic energy budget (DEB) theory; *Spodoptera frugiperda*; development; temperature; variability

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

check for
updates

**Citation:** Gergs, A.; Baden, C.U. A Dynamic Energy Budget Approach for the Prediction of Development Times and Variability in *Spodoptera frugiperda* Rearing. *Insects* **2021**, *12*, 300. https://doi.org/10.3390/insects12040300

Academic Editor: Gaelle Le Goff

Received: 23 February 2021

Accepted: 25 March 2021

Published: 29 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Laboratory experiments are usually designed in a way that allows the study of biological behaviors in response to the variation of one factor only, while keeping the others as constant as possible. Attempts are made to reduce the variability in the initial conditions of the experiments to increase their statistical power. As such, the experimental outcome is directly linked to the rearing conditions of the biological test organisms and their homogenous state.

Due to its widespread and expanding dispersal and its economic and socioeconomic relevance [1], the fall armyworm, *Spodoptera frugiperda* (J.E. Smith) (Lepidoptera, Noctuidae), is an increasingly important laboratory and field model organism in biological and
agricultural research. The fall armyworm is native to the Americas, but in the last few
years, it has expanded its distribution over large parts of Africa, Asia, and Australia [2].
This polyphagous noctuid species causes major economic damage, as the 350 known host
plants species includes many important crops as well, e.g., corn, rice, sorghum, cotton, and
many varieties of vegetables [3]. In addition to its polyphagous abilities, the fall armyworm
has other properties that make it a severe pest. These include the absence of a diapause, a
high reproductive rate, and high migratory abilities [4–6].

The artificial rearing of this important pest species is crucial to obtain more information about its biology, behavior, metabolic rate, and all the other puzzle pieces needed for an integrated pest management program [7]. This becomes even more urgent because of the expanding distribution range of *S. frugiperda*; therefore, rearing conditions and artificial diets are the topics of numerous current publications [7–9]. The important influence of the temperature regime on the development times of different larval stages was recently shown by Du Plessis et al. [8], as well as Lopez et al. [10] and Montezano et al. [3].

In this paper, we illustrate how resource availability and ambient temperature may affect larval development and variability by means of dynamic energy budget (DEB) modeling. DEB theory is based on the idea that the principle of energy metabolism is largely conserved across species. DEB models can, thus, be applied to simulate the entire life cycle of a wide range of animals [11], including insects [12–14]. First proposed by Kooijman [15], DEB theory has previously been applied to analyze the effects of environmental factors on life history processes [16,17], including temperature [18] and food availability in a population context [19,20].

In this paper, we aimed to analyze how resource availability in laboratory cultures and rearing practices, such as growth-controlling temperatures, affect larval development and variability in *S. frugiperda*.

**2. Materials and Methods**

DEB models quantify the rates at which organisms assimilate energy from the environ-
ment and subsequently allocate this energy via a reserve compartment to structural growth,
the reproductive system, and the maintenance of bodily functions. The models describe
the entire life cycle of an organism starting with the embryonic phase to quantify somatic
(structural) growth, maturation, and reproduction. Maturity thresholds for birth and pu-
berty mark the onset of feeding and investment of energy into reproduction, respectively.
As all modeled rates are dependent on temperature and energy acquisition as a function
of food availability, both environmental factors are key drivers in DEB model systems.
Many decades of research in DEB theory have led to the development of the Add-my-Pet
(AmP) database, where data, models, and sets of parameters for more than 2800 species
have currently been collected. Within the collection, the parameter set for a given species is
available and can be used in combination with the well-tested model code that is used for
the parameter estimation. All of the data used for estimation, as well as the assumptions
used to link data and the model, are well documented and openly accessible. Herein, we
used a DEB model for *S. frugiperda* previously been developed and published via the AmP
collection [21]. The model makes use of the “hex” variant of the DEB system to cover
the life cycles of holometabolic insects and some other hexapods. The model covers the
morphological life stages of the egg, larva, pupa, and the imago, as well as the functional
stages, including the embryo and reproductive stages. As a deviation from the standard
model [22], the larval stage accelerates in terms of growth [23] and allocates energy to
reproduction. Larval stage transitions are marked by instar-specific stresses, depending on
their structural length. For further details of the model implementation and parameters,
see [21]. Exemplified code is available from the Supplementary Material S2.

For all scenario analyses and the simulation of developmental variability, we included
stochasticity in the model: We multiplied the maximum specific assimilation flux (pAm),
which is one of the primary DEB parameters, with a random factor drawn from a normal
distribution with a mean of one. The range of variability for a given setting was explored
by means of Monte-Carlo simulations, i.e., repeated simulation runs for different random
factors. Note that by posing variability on the parameter *p*<sub>Am</sub> and employing a constant
value for the scaled functional response *f*, we affectively assumed that the assimilation
efficiency differs among individual larvae.

For model validation, we mimicked the variable laboratory conditions as employed
in an independent experiment and compared the model output and data. The culture
conditions and experimental set-up were as follows. In *S. frugiperda* rearing, we used
cylindrical flight cages (diameter, 24 cm; height, 26 cm) containing 200 adults (male/female
ratio, approximately 1:1) fed with a sugar solution. The temperature was 25 °C, humidity
50%, and the flight cages were transparent, allowing exposure to a dimmed natural light
regime. Females deposited their eggs in batches on a piece of paper kitchen towel, forming
the lid of the cylindrical flight cages. Eggs were allowed to mature and hatch at 26 °C and
in darkness in plastic tubs (18 cm length × 13.5 cm width × 6 cm height) containing a
layer of artificial diet. At the second larval stage, the larvae were separated into individual
Petri dishes (diameter, 5.5 cm) and cultured at 25 °C, 60% humidity, and in darkness until
pupation. The recipe for 1 L of the artificial diet was: 0.43 L of hot water to solve 14.5 g
of agar and 18 g of alfalfa powder, and then after adding 0.325 L of cold water, the other
ingredients, 0.5 mL of rapeseed oil, 40 g of baker's yeast, 0.7 g of Wesson Salt Mix, 1 g of
β-Sistosterol, 0.5 g of L-Leucin, 4 g of ascorbic acid, 3.7 g of vitamin mix, 3.6 g of sorbic
acid, 133 g of bean flour, 0.5 g of antibiotics, and 3 mL of 4- hydroxybenzene S were added.
The validation experiment was performed using a variable-temperature regime and was
started with one of the mass rearing plastic containers derived from the laboratory culture.
To prevent cannibalism [24], 36 larvae were separated during the second larval instar into
small Petri dishes (5.5 cm in diameter) filled with an ad libitum amount of artificial diet.
These Petri dishes were shifted between different temperature regimes using a series of
climate cabinets maintained at 22, 26, and 30 °C. The larval stage attained was determined
individually once to twice a day.

**3. Results**

We used a DEB model for *S. frugiperda* to simulate larval development times under different rearing conditions. The simulated larval development times, on average, ranged from 9.0 to 31.2 days at 32 °C and 18 °C, respectively (Figure 1), assuming ad libitum food conditions. In order to account for the variability in development times that were usually observed in laboratory experiments, we estimated a standard deviation in the maximum specific assimilation flux, one of the primary DEB parameters, of 0.087 J d⁻¹ cm⁻², based on the empirical range observed by Du Plessis et al. [8]. The resulting variability in larval development times is shown in Figure 1.

For validation, model predictions were compared to experimental data (Supplemen-
tary Material S1) that had not been used for model development. In this laboratory
experiment, ambient temperature was varied between 22 and 30 °C. Starting from the first
instar, under these conditions, all larvae reached the sixth instar within 16 days. The DEB
model accurately predicted the larval stage transitions over time—most of the measured
data points were within the prediction interval (Figure 2), indicative of a good model
performance. Interestingly, the modeled variability increased within the fourth instar,
coinciding with a drop in ambient temperature at day 10 of the experiment.
<figure><img src="image_55.png" /><figcaption>**Figure 1.** Simulated (dots) and measured range (dashed lines; Du Plessis et al. [8]) of the individual larval development times in days (d) in *S. frugiperda*.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>**Figure 2.** Frequency of the larval instars measured over time in days (d) for an exemplified variable temperature scenario (red line). The median and 95% confidence limits of the model prediction are represented by the blue lines and shade areas, respectively. Dots indicate measured data. The instar frequencies were recorded twice per day from instar 3.</figcaption></figure>

In the scenario analyses, we subsequently analyzed the variability in larval development times with respect to food availability and fluctuating temperatures in more detail. We varied the value for the scaled functional response in the DEB model to simulate different experimental food conditions and observed the model outputs in terms of development times, starting from the first instar larvae. The scaled functional response can have values between f = 1 (ad libitum) and f = 0 (absence of food) and might represent different food
qualities or quantities. The model simulations revealed that variability in larval develop-
ment times may increase with decreasing food quality or quantity, as illustrated for the
fifth instar larvae in Figure 3. Moreover, the simulated development times increased with
increasing values for the scaled functional response, as did the standard error (Figure 5).

<figure><img src="image_57.png" /><figcaption>**Figure 3.** Simulated development times in days (d) (n = 1000) for the example of the fifth larval instar reared at 26 °C and three different feeding regimes: Ad libitum (f = 1), intermediate (f = 0.6), and low (f = 0.4) food availability. Blue columns represent the frequencies of the development times of larvae (n = 1000), and the red vertical line indicates the mean development time per instar.</figcaption></figure>

The simulated larval development times were also affected by up-/downregulation
of ambient temperatures (Figure 4). We decreased the simulated temperatures between 30
and 22 °C just before the transition into the third larval instar, and vice versa. This model
analysis revealed that the downregulation of ambient temperatures may result in increased
variability (Figure 4, left), while the opposite was the case for increasing temperature
regimes (Figure 4, right).

**4. Discussion**

The DEB model provides an excellent basis for the prediction of the developmental be-
havior in *S. frugiperda* rearing, particularly regarding changes in temperature and provision
of nutrients. The presented model adequately predicts development times in a changing
rearing environment and allows for the analysis of variability among individuals.

Increased variability, with longer development times at lower ambient temperatures, were reported for *S. frugiperda* by Du Plessis et al. [8] (Figure 1) and Garcia et al. [25]. The authors of these studies fitted linear regression models to temperature data for different stages to evaluate the development times and to estimate the temperature thresholds and the number of degree-days. However, Du Plessis et al. [8] argued that the developmental rates become non-linear at unfavorable temperatures. Deviations of the linear temperature model were thus employed, e.g., by Garcia et al. [26], to describe the development rates in *S. frugiperda* beyond empirical observations.

DEB models, by contrast, are based on first principles of bioenergetics. They capture
basic life history processes, such as feeding, development, growth, maintenance of bodily
functions, reproduction, and senescence under environmental fluctuations in tempera-
ture and food availability in one coherent framework and a relatively low number of
parameters [11]. The core concept of DEB theory is consistent with the principles of thermo-
dynamics [27,28] and general trends in evolutionary history [29]. The modular structure of
DEB models allows specific environmental attributes or stressors to be accounted for with-
out changing the core of the theory. In DEB models, changes in life history trajectories, such
as growth and development in response to environmental temperatures, are considered to
be a result of changes in bioenergetic rate constants. It has been demonstrated by means
of DEB modeling that different life history processes, such as food ingestion, growth, and
reproduction, share the same Arrhenius temperature [18].
<figure><img src="image_58.png" /><figcaption>**Figure 4.** Simulated larval development for decelerating (**left**) and accelerating (**right**) temperature regimes. Blue columns represent the frequencies of the development times of larvae (*n* = 1000), and the red vertical line in the histograms indicates the mean development time per instar.</figcaption></figure>

In the DEB model, the three-parameter Arrhenius model that was parameterized for *S. frugiperda* [21] describes the non-linear temperature response in larval development times well (Figure 1). Moreover, the same Arrhenius parameter described the instar-specific developmental times well, except for the first instar, where a deviating parameter value was needed [21]. The DEB model approach thus requires fewer parameters and less experimental testing compared to stage-specific linear regression approaches. As regards food, the modeled increase in variability with decreasing values for the scaled functional response f (Figure 3) captures the general pattern of standard errors being increased if development times are prolonged—see, e.g., daSilva et al. [30] for different food types such as oats, cotton, or an artificial diet (Figure 5). Note that simulated development times represent the time span from egg hatching until pupation, while the data [30] also cover
the pupa stage, which explains the difference in the development times and the standard errors between the simulation and data. No extra parameters were needed to simulate the different food levels, as the energy acquisition from the environment and subsequent bioenergetic fluxes are inherent properties of the DEB model.

<figure><img src="image_59.png" /><figcaption>**Figure 5.** Standard error as function of development time in days (d). The data cover development times (larvae to adults), as measured by daSilva et al. [30], for different food types (soybean, cotton, maize, oats, wheat, and an artificial diet) at an ambient temperature of 25 ± 2 °C. Simulated mean development times (first instar to pupation) and standard errors (n = 10) were calculated 1000 times for random values of the scaled functional response between f = 0.4 and f = 1 and temperatures in the range of 25 ± 2 °C.</figcaption></figure>

The DEB model revealed that lower temperatures and food sources of low nutritional value lead to more heterogenous development times (Figures 3 and 4), which subsequently result in more heterogenous colonies or cohorts (e.g., different larval stages). We gained these results both via the model and our experimental data independently of one another (Figure 4). As these effects lead to more homogenous or heterogenous colonies even within a rather short time (Figure 4), every change in a temperature regime should be conducted thoughtfully. For example, the lowering of temperatures overnight (e.g., for slowing down the growth rate to ensure the desired larval stage for a given experiment day) would lead to a higher developmental diversity within the cohort, whereas a higher temperature would have the opposite effect. Therefore, for more homogenous rearing, conditions should include a rather high temperature, between 26 and 30 °C [8], with a good quality food source (Figure 3) [7] to ensure reproducible experimental data with less variability in *S. frugiperda* research.

If the temperature conditions should change during rearing, the DEB model for *S. frugiperda* is a useful tool to predict and understand the resulting changes in development times and variability. On the contrary, with the help of the model, rearing conditions can be calculated and predicted beforehand, decreasing the need to interfere in running larval rearing. The DEB model can therefore enable researchers to obtain the desired larval stage at a specific date.

A further potential application of the model is for estimating the development rates in field populations of *S. frugiperda*, providing an alternative calculation type besides the degree-days model [8], if there is sufficient information about the temperature and host plant nutrients available (e.g., *S. frugiperda* larval development times on cotton are much longer than on maize [31]). Therefore, DEB models provide important tools for future insect pest management research in laboratory (rearing) and field population biology. They could enhance risk management tools, e.g., for seasonal multigenerational migration, to further improve existing population dynamics and migration models [25,26]. For instance, the NicheMapR package can be used to make field calculations as a function of microclimate using the DEB model [32].
**5. Conclusions**

The dynamic energy budget model is a powerful tool to predict the development
times in *S. frugiperda* rearing and can provide guidance on how to adjust the development
times to individual needs. Furthermore, it helps to reduce possible side effects by avoiding
unnecessary developmental variabilities within a larval cohort caused by unfavorable tem-
perature or nutritional regimes. Many potential interactions with other models predicting
population dynamics in the field can be envisaged.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/insects12040300/s1: S1: Validation data and S2: model code.

**Author Contributions:** Conceptualization, A.G. and C.U.B.; methodology, A.G. and C.U.B.; software, A.G.; validation, A.G. and C.U.B.; formal analysis, A.G.; investigation, A.G. and C.U.B.; data curation, A.G.; writing—original draft preparation, A.G. and C.U.B.; writing—review and editing, A.G. and C.U.B.; visualization, A.G. All authors read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Validation data (Figure 2) are available in the Supplementary Materials.

**Acknowledgments:** The authors thank Petra Papachristou, Dorina Pöttgen, Tina Teiwes, and Markus Schaeben for the maintenance of the *Spodoptera frugiperda* rearing and their technical support.

**Conflicts of Interest:** The authors are employed by Bayer AG. Bayer products were used in the fall armyworm control. The authors declare no additional conflicts of interest.

**References**

1. EPPO. Available online: https://gd.eppo.int/taxon/LAPHFR/distribution (accessed on 8 December 2020).

2. FAO. The Global Action for Fall Armyworm Control: Action framework 2020–2022. In *Working Together to Tame the Global Threat*; Food and Agriculture Organization of the United Nations: Rome, Italy, 2020. [CrossRef]

3. Montezano, D.G.; Specht, A.; Sosa-Gómez, D.R.; Roque-Specht, V.F.; Sousa-silva, J.C.; Paula-Moraes, S.V.; Peterson, J.A.; Hunt, T.E. Host plants of *Spodoptera frugiperda* (Lepidoptera: Noctuidae) in the Americas. *Afr. Entomol.* **2018**, *26*, 286–300. [CrossRef]

4. Pair, S.D.; Raulston, J.R.; Sparks, A.N.; Westbrook, J.K.; Douce, G.K. Fall armyworm distribution and population dynamics in the southeastern states. *Fla. Entomol.* **1986**, *69*, 468–487. [CrossRef]

5. Johnson, S.J. Migration and the life history strategy of the fall armyworm, *Spodoptera frugiperda* in the western hemisphere. *Int. J. Trop. Insect Sci.* **1987**, *8*, 543–549. [CrossRef]

6. Westbrook, J.; Fleischer, S.; Jairam, S.; Meagher, R.; Nagoshi, R. Multigenerational migration of fall armyworm, a pest insect. *Ecosphere* **2019**, *10*, 1–12. [CrossRef]

7. Pinto, J.R.L.; Torres, A.F.; Truzi, C.C.; Vieira, N.F.; Vacari, A.M.; De Bortoli, S.A. Artificial Corn-Based Diet for Rearing *Spodoptera frugiperda* (Lepidoptera: Noctuidae). *J. Insect Sci.* **2019**, *19*, 1–8. [CrossRef]

8. Du Plessis, H.; Schlemmer, M.-L.; Van den Berg, J. The Effect of Temperature on the Development of *Spodoptera frugiperda* (Lepidoptera: Noctuidae). *Insects* **2020**, *11*, 228. [CrossRef]

9. Jin, T.; Lin, Y.; Chi, H.; Xiang, K.; Ma, G.; Peng, Z.; Yi, K. Comparative Performance of the Fall Armyworm (Lepidoptera: Noctuidae) Reared on Various Cereal-Based Artificial Diets. *J. Econ. Entomol.* **2020**, *116*, 2986–2996. [CrossRef] [PubMed]

10. Lopez, R.Y.; Ortega, A.V.; Centeno, J.H.A.; Ruiz, J.S.; Carranza, J.A.Q. Sistema de alerta contra el gusano cogollero *Spodoptera frugiperda* (J.E. Smith) (Insecta: Lepidoptera: Noctuidae). *Rev. Mex. Cienc. Agríc.* **2019**, *10*, 405–416.

11. Kearney, M.R.; Domingos, T.; Nisbet, R. Dynamic Energy Budget Theory: An Efficient and General Theory for Ecology. *BioScience* **2015**, *65*, 34. [CrossRef]

12. Gergs, A.; Jager, T. Body size-mediated starvation resistance in an insect predator. *J. Anim. Ecol.* **2013**, *83*, 758–768. [CrossRef] [PubMed]

13. Maino, J.; Kearney, M. Testing mechanistic models of growth in insects. *Proc. R. Soc. B* **2015**, *282*, 1–9. [CrossRef]

14. Llandres, A.; Marques, G.; Maino, J.; Kooijman, S.A.L.M.; Kearney, M.; Casas, J. A Dynamic Energy Budget for the whole life-cycle of holometabolous insects. *Ecol. Monogr.* **2015**, *85*, 353–371. [CrossRef]

15. Kooijman, S.A.L.M. Energy budgets can explain body size relations. *J. Theor. Biol.* **1986**, *121*, 269–282. [CrossRef]

16. Pieters, B.J.; Jager, T.; Kraak, M.H.S.; Admiraal, W. Modeling responses of Daphnia magna to pesticide pulse exposure under varying food conditions: Intrinsic versus apparent sensitivity. *Ecotoxicology* **2006**, *15*, 601–608. [CrossRef] [PubMed]

17. Muller, E.B.; Nisbet, R.M. Dynamic energy budget modeling reveals the potential of future growth and calcification for the coccolithophore Emiliania huxleyi in an acidified ocean. *Glob. Chang. Biol.* **2004**, *20*, 2031–2038. [CrossRef]

18. Kooijman, S.A.L.M. *Dynamic Energy Budget Theory for Metabolic Organization*; Cambridge University Press: Cambridge, UK, 2010; p. 514.
19. Martin, B.T.; Jager, T.; Preuss, T.G.; Nisbet, R.; Grimm, V. Predicting population dynamics from the properties of individuals: A cross-level test of Dynamic Energy Budget theory. *Am. Nat.* **2013**, *181*, 506–519. [CrossRef]

20. Gergs, A.; Preuss, T.G.; Palmqvist, A. Double Trouble at High Density: Cross-Level Test of Resource-Related Adaptive Plasticity and Crowding-Related Fitness. *PLoS ONE* **2014**, *9*, e91503. [CrossRef] [PubMed]

21. Gergs, A.; Baden, C.U. AmP Spodoptera frugiperda, version 2020/07/21. 2020. Available online: https://www.bio.vu.nl/thb/deb/delab/add_my.pet/entries_web/Spodoptera-frugiperda/Spodoptera-frugiperda_res.html (accessed on 8 December 2020).

22. Kooijman, S.A.L.M.; Sousa, T.; Pecquerie, L.; Van der Meer, J.; Jager, T. From food dependent statistics to metabolic parameters, a practical guide to the use of dynamic energy budget theory. *Biol. Rev.* **2008**, *83*, 533–552. [CrossRef] [PubMed]

23. Kooijman, S.A.L.M. Metabolic acceleration in animal ontogeny: An evolutionary perspective. *J. Sea Res.* **2014**, *94*, 128–137. [CrossRef]

24. Chapman, J.W.; Williams, T.; Escribano, A.; Caballero, P.; Cave, R.D.; Goulson, D. Fitness consequences of cannibalism in the fall armyworm, *Spodoptera frugiperda*. *Behav. Ecol.* **1999**, *10*, 298–303. [CrossRef]

25. Garcia, A.; Godoy, W.A.C.; Thomas, J.M.C.; Nagoshi, R.N.; Meagher, R. Delimiting Strategic Zones for the Development of Fall Armyworm (Lepidoptera: Noctuidae) on Corn in the State of Florida. *J. Econ. Entomol.* **2017**, *111*, 120–126. [CrossRef] [PubMed]

26. Garcia, A.G.; Ferreira, C.P.; Godoy, W.A.C.; Meagher, R. A computational model to predict the population dynamics of *Spodoptera frugiperda*. *J. Pest. Sci.* **2019**, *92*, 429–441. [CrossRef]

27. Jusup, M.; Sousa, T.; Domingos, T.; Labinac, V.; Marn, N.; Wang, Z.; Klanjscek, T. Physics of metabolic organization. *Phys. Life Rev.* **2017**, *20*, 1–39. [CrossRef] [PubMed]

28. Sousa, T.; Domingos, T.; Poggiale, J.C.; Kooijman, S. Dynamic energy budget theory restores coherence in biology Introduction. *Philos. Trans. R. Soc. Lond. B Biol. Sci.* **2010**, *365*, 3413–3428. [CrossRef] [PubMed]

29. Kooijman, S.A.L.M.; Troost, T.A. Quantitative steps in the evolution of metabolic organisation as specified by the Dynamic Energy Budget theory. *Biol. Rev.* **2007**, *82*, 113–142. [CrossRef] [PubMed]

30. Da Silva, D.M.; de Freitas Bueno, A.; Andrade, K.; dos Santos Stecca, C.; Oliveira, P.M.; Neves, J.; de Oliveira, M.C.N. Biology and nutrition of *Spodoptera frugiperda* (Lepidoptera: Noctuidae) fed on different food sources. *Sci. Agric.* **2017**, *74*, 18–31. [CrossRef]

31. Ali, A.; Luttrell, R.G.; Schneider, J.C. Effects of temperature and larval diet on development of the fall armyworm (Lepidoptera: Noctuidae). *Ann. Entomol. Soc. Am.* **1990**, *83*, 725–733. [CrossRef]

32. Kearney, M.R.; Porter, W.P. NicheMapR—An R package for biophysical modelling: The ectotherm and Dynamic Energy Budget models. *Ecography* **2020**, *43*, 85–96. [CrossRef]
<figure><img src="image_60.png" /></figure>
Article

# Geographic Monitoring of Insecticide Resistance Mutations in Native and Invasive Populations of the Fall Armyworm

Sudeeptha Yainna<fnref n="1" />,<fnref n="2" /> ID, Nicolas Nègre<fnref n="1" /> ID, Pierre J. Silvie<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> ID, Thierry Brévault<fnref n="3" />,<fnref n="4" /> ID, Wee Tek Tay<fnref n="5" />, Karl Gordon<fnref n="5" />, Emmanuelle dAlençon<fnref n="1" /> ID, Thomas Walsh<fnref n="5" />,* and Kiwoong Nam<fnref n="1" />,* ID

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

check for
updates

**Citation:** Yainna, S.; Nègre, N.; Silvie, P.J.; Brévault, T.; Tay, W.T.; Gordon, K.; dAlençon, E.; Walsh, T.; Nam, K. Geographic Monitoring of Insecticide Resistance Mutations in Native and Invasive Populations of the Fall Armyworm. *Insects* **2021**, *12*, 468. https://doi.org/10.3390/insects12050468

Academic Editor: Béla Darvas

Received: 31 March 2021

Accepted: 12 May 2021

Published: 18 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Simple Summary:** The moth fall armyworm (*Spodoptera frugiperda*) is a major agricultural pest insect damaging a wide range of crops, especially corn. Field evolved resistance against *Bacillus thuringiensis* (Bt) toxins and synthetic insecticides has been repeatedly reported. While the fall armyworm is native to the Americas, its biological invasion was first reported from West Africa in 2016. Since then, this pest has been detected across sub-Saharan and North Africa, Asia, and Oceania. Here, we examine the geographical distribution of mutations causing resistance against Bt or synthetic insecticides to test if the invasion was accompanied by the spread of resistance mutations using 177 individuals collected from 12 geographic populations including North and South America, West and East Africa, India, and China. We observed that Bt resistance mutations generated in Puerto Rico or Brazil were found only from their native populations, while invasive populations had higher copy numbers of cytochrome P450 genes and higher proportions of resistance mutations at AChE, which are known to cause resistance against synthetic insecticides. This result explains the susceptibility to Bt insecticides and the resistance against synthetic insecticides in invasive Chinese populations. This information will be helpful in investigating the cause and consequence associated with insecticide resistance.

**Abstract:** Field evolved resistance to insecticides is one of the main challenges in pest control. The fall armyworm (FAW) is a lepidopteran pest species causing severe crop losses, especially corn. While native to the Americas, the presence of FAW was confirmed in West Africa in 2016. Since then, the FAW has been detected in over 70 countries covering sub-Saharan Africa, the Middle East, North Africa, South Asia, Southeast Asia, and Oceania. In this study, we tested whether this invasion was accompanied by the spread of resistance mutations from native to invasive areas. We observed that mutations causing Bt resistance at ABCC2 genes were observed only in native populations where the mutations were initially reported. Invasive populations were found to have higher gene numbers of cytochrome P450 genes than native populations and a higher proportion of multiple resistance mutations at acetylcholinesterase genes, supporting strong selective pressure for resistance against synthetic insecticides. This result explains the susceptibility to Bt insecticides and resistance to various synthetic insecticides in Chinese populations. These results highlight the necessity of regular and standardized monitoring of insecticide resistance in invasive populations using both genomic approaches and bioassay experiments.

**Keywords:** ABCC2; *Bacillus thuringiensis*; biological invasion; Cytochrome P450; Fall armyworm; insecticide resistance; *Spodoptera frugiperda*

<footnote n="1">¹ DGIMI, University of Montpellier, INRAE, F-34095 Montpellier, France; sudeepha.yainna-kumarihami@inrae.fr (S.Y.); nicolas.negre@umontpellier.fr (N.N.); emmanuelle.d-alencon@inrae.fr (E.d.)</footnote>

<footnote n="2">² CIRAD, UPR AIDA, F-34398 Montpellier, France; pierre.silvie@cirad.fr</footnote>

<footnote n="3">³ PHIM Plant Health Institute, University of Montpellier, IRD, CIRAD, INRAE, Institut Agro, F-34398 Montpellier, France; thierry.brevault@cirad.fr</footnote>

<footnote n="4">⁴ AIDA, University of Montpellier, CIRAD, F-34398 Montpellier, France</footnote>

<footnote n="5">⁵ Black Mountain Laboratories, CSIRO, Canberra, ACT 2601, Australia; weetek.tay@csiro.au (W.T.T.); karl.gordon@csiro.au (K.G.)</footnote>

<footnote>* Correspondence: tom.walsh@csiro.au (T.W.); ki-woong.nam@inrae.fr (K.N.)</footnote>
# 1. Introduction

Insecticide resistance is one of the main challenges for the control of insect pests. Commonly used insecticides to control pest insects can be classified into two main types. The first group is Bt (*Bacillus thuringiensis*) toxins, which are generally produced by transgenic crops. The global areas of planted Bt crops is positively correlated with the number of Bt resistance pest species, implying that field-evolved Bt resistance is prevalent [1]. Field-evolved Bt resistance is a particularly serious issue in corn and cotton because the majority of those plants are Bt crops in the USA and because only GM-cotton expressing Bt toxins are grown in Australia. Field-evolved Bt resistance has been observed in several major insect Lepidoptera pest species including maize stalk borer (Busseola fusca), western corn rootworm (Diabrotica virgifera virgifera), cotton bollworm (Helicoverpa zea), Old World cotton bollworm (Helicoverpa armigera), Helicoverpa punctigera, fall armyworm (Spodoptera frugiperda), and pink bollworm (Pectinophora gossypiella) [1]. The other group is synthetic insecticides such as organophosphates, pyrethroids, neonicotinoids, organochlorides, carbamates, ryanoids, and spinosyns. By 2021, more than 600 arthropod species have been included in the arthropod pesticide resistance database (https://www.pesticideresistance.org/) (accessed on 21 January 2021).

Several genetic mechanisms of Bt resistance have been reported (reviewed in [2]). If a mutation in a midgut receptor causes reduced physical binding with Bt proteins, the toxicity of Bt can be decreased. These receptors include cadherin, aminopeptidase-N, alkaline phosphatases, and ATP-binding cassette transporters. Alternatively, altered processing of Bt prototoxin, increased immune status, sequestration of Bt protein, and accelerated recovery of epithelial integrity can also increase Bt resistance. The resistance against synthetic insecticides can be caused by detoxification genes [3] such as cytochrome P450 gene, esterase, and glutathione S transferase. In addition, mutations in calcium channels, acetylcholinesterase, and nicotinic acetylcholine receptors may also cause resistance [4–7].

Fall armyworm (FAW, *Spodoptera frugiperda*, J.E.Smith) (Lepidoptera: Noctuidae: Noctuinae) is one of the most damaging pest insects, partly due to the extreme polyphagy by feeding on at least 353 species in 76 plant families, which include several economically important crops such as corn, rice, sorghum, sugarcane, cotton, and soybean [8]. Genomic analyses have demonstrated a rapid expansion of detoxification genes [9–11], potentially associated with this extreme polyphagy by overcoming plant defense toxins from diverse plants.

Field-evolved resistance to Bt insecticide in the FAW has been reported from Puerto Rico [12,13], Brazil [14,15], Argentina [16], and the USA [17]. Five mutations causing Bt resistance have been reported. Interestingly, all these mutations were observed in the ATP Binding Cassette Subfamily C Member 2 (ABCC2) gene, which encodes ATP-binding cassette transporter proteins, implying parallel evolution. In a population from Puerto Rico, 2bp GC insertion at ABCC2 causes a frameshift mutation and a premature stop codon, which results in improper binding between ABCC2 and Bt Cry1F and Cry1A proteins [18,19]. Boaventura et al. identified three causal mutations of Cry1F resistance at ABCC2 including GY deletion, P799K/R substitution, and G1088D substitution from a Brazilian population [20]. Guan et al. also identified a 12bp insertion, which causes a frameshift mutation and a premature stop codon, from another Brazilian population [21]. As well as the Bt resistance, field-evolved resistance to synthetic insecticides has been widely reported. For example, a population in Puerto Rico was reported to have resistance against a wide range of synthetic insecticides including flubendiamide, chlorantraniliprole, methomyl, thiodicarb, permethrin, chlorpyriphos, zeta-cypermethrin, deltamethrin, triflumuron, and spinetoram [22]. We observed that P450 genes cause resistance against deltamethrin [23]. Boaventura et al. reported that I4734M and G4946E substitutions at the ryanodine receptor (RyR) cause resistance against diamides [24]. Carvalho et al. reported that resistance against organophosphate and carbamates is caused by A314S, G340A, and F402V mutations at acetylcholinesterase (AChE) and that pyrethroid resistance is caused by T929I, L932F, and L1014F mutations at the voltage-gated sodium channel (VGSC) gene [4].
FAW is native to North and South America, and its presence was first confirmed in West Africa in 2016 [25]. Since then, FAW has been detected across almost the entire Sub-Saharan Africa, in northern Africa/the Near East (Egypt), Middle East Asia, South Asia, South East Asia, East Asia, Australia, and more recently, the Canary Islands located off the coast between Western Africa and Morocco, and New Caledonia in the South Pacific (https://www.cabi.org/ISC/fallarmyworm) (accessed on 21 January 2021). The invasion caused a serious reduction in corn production, between 21–53% in Africa [26]. While native populations are composed of two sympatric strains, corn (sfC) and rice strains (sfR), named after their preferred host-plants [27–29], the invasive populations are observed as hybrids between corn and rice strains [30–32], but are predominantly observed in corn fields. Genomic analyses support the multiple introductions of FAW from multiple native populations including Florida and Brazil to multiple invasive areas including China and Africa [31,33]. However, the invasive populations exhibited much more homogeneous DNA sequences than native populations, implying that these invasive populations had experienced extensive admixture since their establishment across the Old World [31].

The invasion process opens up the possibility that resistance mutations have been
spread as well in the invasive area. Recently, Boaventura et al. observed that resistance
mutations at AChE were frequently observed from two invasive populations in Kenya and
Indonesia and two populations in Brazil and Puerto Rico, while GY deletion of ABCC2
genes was observed only in the Brazilian population [34]. Boaventura et al. also reported
the presence of the VGSC L1014F allele with low frequency in an Indonesian population
that has, to-date, yet to be reported in other invasive populations [34]. Zhang et al. reported
that the 2 bp GC insertion was not observed from Chinese populations [32], while Guan
et al. reported that the 12 bp insertion in exon 15 of the ABCC2 gene was not detected
in FAW populations from Uganda, Malawi, and from six Chinese provinces [21]. How-
ever, we still lack a comprehensive understanding of the geographic distribution of these
resistance mutations because global profiling of the causal resistance mutations has not
yet been performed. We performed whole-genome resequencing from 177 individuals
from both invasive and native populations sampled from 12 geographic populations in
North America, South America, Africa, India, and China to investigate adaptive evolution
associated with the invasive success in our previous study [31]. In this study, we undertook
reanalysis of the resequencing data of the FAW populations to test the presence of causal
resistance mutations against insecticides in the FAWs.

2. Methods

2.1. *Resequencing Data*

Resequencing data generated in our previous study were reused [31,33]. More specifi-
cally, this dataset included 29 sfC and 49 sfR individuals from invasive populations (Benin,
Malawi, Uganda, China, and India) and 70 sfC and 29 sfR individuals from native popula-
tions (Mississippi, Florida, Puerto Rico, Brazil, French Guinea, Guadeloupe, and Mexico)
(Figure 1) (please see Table S1 for more information). This dataset includes four Bt-resistant
and six susceptible Brazilian individuals. Paired-end whole genome resequencing was
performed with approximately 20X coverage with 300 bp insert size, and 150 bp read length
using Illumina (San Diego, CA, USA) HiSeq 2500, HiSeq 4000, and Novaseq 6000. Variant
calling was performed using the GATK-4.0.11.0 package [35] for single nucleotide variant
(SNV) and CNVCaller [36] for copy number variation (CNV). The numbers of variations
were 27,117,672 and 22,916 for SNV and CNV, respectively.
<figure><img src="image_63.png" /><figcaption>**Figure 1.** Countries where the analyzed individuals were sampled. The green color indicates native populations. Purple (2016), orange (2017), and yellow colors (2018) indicate the reported years of detection in countries from where individuals were collected. The map was generated using mapchart [37].</figcaption></figure>

## 2.2. Gene Annotation and Statistical Analysis

ABCC2, RyR, AChE, and VGSC were annotated by aligning protein sequences obtained from NCBI and reference genomes using exonerate 2.2.0 [38] with the protein2genome model. The accession numbers of these genes were QGS83596, XP_022819835, MK226188, and KC435025, for ABCC2, AChE, RyR, and VGSC, respectively. Genetic variations within these genes were identified from the resequencing data [31] using IGV tools [39]. P450 genes were identified from the gene annotation files generated in our previous study (OGS7.0) [31]. The identification of the P450 clan was performed by blasting against OGS2.2 [9], which includes manually annotated P450 genes, with 98% cutoff using ncbi-blast-2.10.0+ [40]. The Fisher's exact test was performed using the R package to test a statistical significance of increased gene duplication of P450 genes. Weir and Cockerham's FST [41] between Bt susceptible and resistant individuals was calculated at the ABCC2 gene using vcftools v0.1.13 [42]. Increased FST at ABCC2 was tested by calculating the proportion of FST at randomly chosen loci with the same length of ABCC2 from whole genome sequences with 100,000 replications.

# 3. Results

## 3.1. *Bt Insecticide Resistance-ABCC2 Gene*

We tested the presence of five reported mutations at the ABCC2 gene causing the resistance against Bt insecticide. These mutations include (i) 2 bp GC insertion leading to a premature stop codon due to the frameshift (originally identified from a population in Puerto Rico) [18,19]; (ii) GY deletion (Brazil) [20]; (iii) 12 bp insertion causing premature stop codon due to frameshift (Brazil) [21]; (iv) P799K/R substitution (Brazil) [20]; and (v) G1088D substitution (Brazil) [20] (Figure 2). We identified the 2 bp GC insertion from 11 individuals from Puerto Rico in the resequencing data (Table 1). As this mutation was originally reported from a population in Puerto Rico, this result supports that this mutation remained in this original population in Puerto Rico and did not spread to the Old World by invasion.
<figure><img src="image_64.png" /><figcaption>**Figure 2.** Full amino acid sequence of the ABCC2 gene with information on all mutations studied. The green and orange arrows indicate the identified and unidentified mutations from the resequencing data, respectively.</figcaption></figure>

<table><caption>**Table 1.** Individuals with observed resistant mutations (* indicates individuals heterozygous for the mutations). PR: Puerto Rico, CC: susceptible Brazilian individuals, rCC: resistant Brazilian individuals.</caption><thead><tr><th>Mutations</th><th>Individuals with Resistance Mutations</th></tr></thead><tbody><tr><td>2 bp insertion and a frameshift mutation</td><td>* PR1, PR5, PR12, PR14, PR16, * PR18, * PR19, * PR27, PR30, * PR31, * PR33</td></tr><tr><td>GY deletion</td><td>CC44, CC69</td></tr><tr><td>P799K/R</td><td>CC44, CC69</td></tr><tr><td>12 bp insertion and a frameshift mutation</td><td>* rCC25, * rCC5</td></tr></tbody></table>

GY deletion was identified from two individuals in the Brazilian population, and P799K/R was observed from the same individuals. This result implies that this mutation also remained in the original population. Intriguingly, these mutations exist at the homozygous state in these two Bt susceptible individuals, suggesting a possibility that these mutations alone are not sufficient for Bt resistance.

We identified two individuals with 12 bp insertion only from Brazil, also supporting that this mutation remains in the original population and confirms the analysis of Guan et al. [21]. G1088D was not observed from the resequencing data.

As unidentified mutations causing Bt resistance might exist, we investigated sequences with complete genetic differentiation between resistant and susceptible Brazilian individuals. In total, 49 SNPs had complete genetic differentiation (i.e., F<sub>ST</sub> = 1), and all these positions were found to be intronic. This result shows that other causal mutations are not likely to exist in the tested individuals.

We tested the possibility that natural selection on the ABCC2 gene is responsible for the field-evolved Bt resistance in the Brazilian population. In this case, ABCC2 is expected to have increased genetic differentiation between resistant and susceptible groups. F<sub>ST</sub> calculated between these two groups was 0.2922 at ABCC2 gene. In total, 7.7% of 100,000 randomly chosen loci with the same length of the ABCC2 gene had F<sub>ST</sub> higher than 0.2922 (equivalent to *p*-value = 0.077) (Figure 3), indicating marginally significant
statistical support for increased genetic differentiation. This result implies that natural selection might have increased the Bt resistance of the Brazilian population.

<figure><img src="image_65.png" /><figcaption>**Figure 3.** *F*<sub>ST</sub> calculated between CC and rCC (red vertical line) and from random grouping among CC and rCC.</figcaption></figure>

## 3.2. *Synthetic Insecticide Resistance*

P450 genes increase the level of synthetic insecticide resistance in FAW, and P450 gene duplication has been suggested as a possible cause underpinning resistance [23]. In total, 99 P450 genes were identified from the reference genome (Table S2). Copy number variation compared with the reference genome was observed from 47 P450 genes, which include4, 12, and 22 genes belonging to clan2, clan3, and clan4, respectively. In all clades, invasive populations have higher numbers of duplicated P450 genes per individual than native populations for all clades (Figure 4). On the other hand, native populations have higher numbers of deleted P450 genes per individuals than invasive populations for all clades. Consequently, invasive populations have overrepresented numbers of duplicated P450 genes compared with deleted ones than the native populations (p-value = 6.018 × 10<sup>-21</sup>, 0.0004553, 0.04414 for clan2, clan3, and clan4, respectively; two-tailed Fisher's exact test). This result supports the hypothesis that invasive populations have higher P450 gene numbers than native populations through CNVs.
<figure><img src="image_66.png" /><figcaption>**Figure 4.** Average numbers of deleted or duplicated P450 genes per individual in invasive and native populations for *clan2*, *clan3*, and *clan4*.</figcaption></figure>

The resequencing data show all three reported resistance mutations at the AChE (Figure 5), as shown by Boaventura et al. [34]. F209V mutations were observed from populations in Uganda, Malawi, Brazil, and China, as reported from Guan et al. [21] as well as India, Benin, Puerto Rico, Mexico, Mississippi, and Florida. Invasive populations have a higher proportion of individuals with resistance F209V mutations (80.26%) than native populations (72.37%) (p = 0.001097, two-tailed Fisher's exact test). A201S mutation was found from populations in Uganda, Malawi, and Brazil, as shown by Guan et al., but also from Benin, Mississippi, and Florida. Invasive populations have a higher proportion of individuals with A201S resistance mutations (20.0%) than native populations (7.14%) (p = 0.01966). Interestingly, two codons were found for the susceptible 'A' in the population from Benin (GCG and GCA). In Benin, the GCG codon was observed from 28 individuals (as well as other populations). The GCA codon was observed from three individuals in Benin. The most parsimonious explanation is that the GCA codon was newly generated from GCG by a point mutation at the third codon position. The G227A mutation was found in Brazil, as per Guan et al. [21] as well as in Puerto Rico and Florida. The proportion of individuals with this resistance mutation was 17.65% in native populations. Resistance mutations at RyR (I4734M and G4946E substitutions) or the VGSC (T929I, L932F, and L1014F substitutions) were not observed in the analyzed populations of FAW.

<figure><img src="image_67.png" /><figcaption>**Figure 5.** The proportion of individuals with susceptibility and resistance at AChE for (A) F209V, (B) A201S, and (C) G227A in invasive and native populations.</figcaption></figure>
**4. Discussion**

The global spread of FAW through invasion may be accompanied by the spread of resistance mutations to Bt or synthetic insecticides to invasive areas. This is known to occur in other species. The population of the cotton bollworm *Helicoverpa armigera* is thought to have invaded South America with a pyrethroid resistance allele, which was then shared with a sister species, *H. zea* [43,44]. Therefore, monitoring the geographical distribution of resistance mutations in the context of invasion can provide useful information for pest FAW control and insecticide resistance management. We identified known resistance mutations from the whole genome sequences of 177 individuals collected from 12 geographic populations in North America, South America, Africa, India, and China. We showed that resistance Bt mutations at ABCC2 remained at their origin in Puerto Rico and Brazil, meaning that these mutations do not spread to invasive areas. However, invasive populations have increased copy numbers of P450 genes compared with native populations through copy number variations. We also observed that invasive populations had increased proportions of resistance mutations at AChE. These results can be best explained by selective pressure causing increased resistance against synthetic insecticides in invasive populations.

Bioassay experiments demonstrated that Chinese populations are susceptible to Bt Cry1Fa toxins [31]. The analyzed Chinese samples in this study did not have mutations causing Cry1Fa resistance, explaining the susceptibility to this Bt toxin. Interestingly, a South African population showed moderate resistance against Bt Cry1A [32], while GC insertion causing Cry1A resistance [18,19] was not identified from invasive populations in the resequencing data. Therefore, the South African population could potentially have a resistance mutation that is yet to be characterized in native populations, and we cannot exclude the possibility of on-going gene flow of unidentified resistance mutations either from native to invasive populations or from other invasive populations that represent independent introduction from South African FAWs. More studies are necessary to identify causal resistance mutation in the South African populations and to identify the origin of this resistance mutation.

If a resistance Bt mutation exists in invasive populations with low frequency, individuals with the resistance mutations might not be included in the analyzed resequencing data. If this is true, existing resistance mutations in invasive populations might be undetected. The probability of not finding a resistance mutation is $(1 - p)^{2N}$, where *p* is the allele frequency of resistance mutations and *N* is the number of diploid organisms. If the allele frequency of resistance Bt mutation is 0.01, representing a rare allele, the probability of missing this mutation in our resequencing data is only $(1 - 0.01)^{2 \times 99} = 13.67\%$. Therefore, we do not believe that a resistance Bt mutation is likely to exist. Even in the case that this possibility of 13.67% is true, the resistance mutation is likely to be eliminated in a population by genetic drift due to the low allele frequency (0.01).

A whole-genome analysis is necessary to contrast the abundance of P450 genes between invasive and native populations, because P450 genes exist in multi-copies scattered across the genome [9–11]. Our whole genome analysis at a population level shows an increased copy number of P450 genes in the invasive FAW populations. Bioassay experiments demonstrated that Chinese populations have resistance to pyrethroid insecticides [32,45]. As P450 genes are known to cause resistance against deltamethrin [23], a type of pyrethroid in FAWs, the increased P450 gene numbers might be responsible for this resistance. Chinese populations also have resistance against organophosphate [32,45] as well as oxadiazine and diamide [32]. As P450 proteins are probably able to detoxify other types of insecticides than pyrethroid, the increased P450 genes may be responsible for the resistance against a wider range of insecticides.

In addition, the increased proportion of resistance mutations at AChE can cause the increased resistance against synthetic insecticides in invasive populations. Multiple resistance mutations with the increased proportion imply strong selective pressure for resistance against synthetic insecticide. The vast majority of beneficial mutations are quickly removed in a population due to stochastic fluctuation of allele frequency by the genetic
bottleneck, and the survival probability of beneficial mutations is just two times that of the selective coefficient. Therefore, the probability of observing a single resistance mutation at AChE is proportional to the extent of increased fitness due to the resistance. The probability of observing multiple resistance mutations is a function of the multiplication of each selective coefficient, ranging between 0 to 1, at each resistance mutation. Therefore, multiple resistance mutations are difficult to observe without very high selective coefficients, and we argue that invasive populations are under strong selective pressure for the resistance against synthetic insecticides.

Glutathione S transferases or UDP glycosyltransferase might also cause resistance against synthetic insecticides. However, the direct relationship between these two gene families and resistance is yet to be tested in FAW. Genomic analyses and bioassay experiments will generate a more comprehensive list of causal resistance mutations against synthetic insecticides.

The differential usage of synthetic insecticides may explain the observations that invasive populations have increased allele frequency of resistance mutations at AchE genes and increased number of P450 genes by CNVs. The average annual insecticide use is particularly high in East Asia including China [46]. We speculate that initially introduced FAW insects in East Asia might have survived through adaptive evolution by the increased P450 gene number or the increased allele frequency of resistance mutations at AchE genes, while these FAWs remained undetected due to the low population sizes before explosive population expansion in West Africa, where genetic admixture occurred between populations with different invasive origins [31]. More studies are necessary to test the role of adaptive evolution associated with synthetic insecticide resistance in invasive success.

The FST result supports natural selection for Bt resistance in the Brazilian population. If selective pressure with the same directionality exists in invasive populations, newly generated Bt resistance mutations can be selected in the near future. Theoretically, species with large population sizes have higher adaptive potential than species with small population sizes because of a higher population-scaled rate of mutation [47,48]. Since the estimated population size of FAWs is as high as 3.75 million [49], the population-scaled mutation rate can be sufficient to generate a genetic variation for drastically rapid Darwinian positive selection. We observed two synonymous codons at a resistance position in AChE from a population Benin, implying that mutational influx could not be a limiting factor of adaptive evolution through positive selection. Therefore, invasive populations might experience adaptive evolution, causing resistance against Bt mutations in the near future.

**5. Conclusions**

In this study, we showed that Bt resistance mutations did not spread during this period to the Old World by the surveyed invasive FAW populations and that invasive populations have increased copy numbers of P450 genes and increased proportions of resistance mutations at AchE compared with native populations. These results imply that invasive populations might be under selective pressure for the resistance against synthetic insecticides.

The samples used in this study were obtained between 2017 to 2018. Since the first reported FAWs in Western Africa in 2016 [31], the genetic information in this study represents snapshots at the early phase of invasions. As Bt crops and synthetic insecticides are commonly used in the invasive area, the level of insecticide resistance can be changed together with underlying genotypes. This change is likely because our result is in line with selective pressure both for Bt and, in particular, synthetic insecticides. Therefore, monitoring changes in field susceptibility to various insecticides will be helpful in FAW control, together with the analysis of changes in genomics sequences to identify known or new causal resistance mutations.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/insects12050468/s1 (Table S1: Sample collection site and year, Table S2: The list of analyzed P450 genes in this study).
**Author Contributions:** Conceptualization, S.Y., T.W., and K.N.; Formal analysis, S.Y. and K.N.; Resources, N.N., P.J.S., T.B., W.T.T., K.G., an E.d.; Writing—original draft preparation, S.Y. and K.N.; Writing—review and editing, N.N., P.J.S., T.B., W.T.T., K.G., E.d., and T.W.; Supervision, K.N., T.B., and E.d.; Funding acquisition, K.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work (ID 1702-018, given to K.N.) was publicly funded through ANR (the French National Research Agency) under the "Investissements d'avenir" program with the reference ANR-10-LABX-001-01 Labex Agro and coordinated by Agropolis Fondation under the frame of I-SITE MUSE (ANR-16-IDEX-0006). In addition, a grant from the department of Santé des Plantes et Environnement at Institut national de la recherche agronomique was awarded to K.N. (adaptivesv). This work was also financially supported by EUPHRESCO (FAW-spedcom, given to Anne-Nathalie Volkoff) and by CSIRO Health & Biosecurity (given to W.T.T., T.W., and K.G.). S.Y. was supported by a CIRAD-INRAE PhD fellowship.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** We will provide all files generated in this study upon request. These files include vcf files, gene annotation files, and computer programming scripts.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Tabashnik, B.E.; Brévault, T.; Carrière, Y. Insect Resistance to Bt Crops: Lessons from the First Billion Acres. *Nat. Biotechnol.* **2013**, *31*, 510–521. [CrossRef]

2. Jurat-Fuentes, J.L.; Heckel, D.G.; Ferré, J. Mechanisms of Resistance to Insecticidal Proteins from *Bacillus thuringiensis*. *Annu. Rev. Entomol.* **2021**, *66*, 121–140. [CrossRef]

3. Hemingway, J.; Field, L.; Vontas, J. An Overview of Insecticide Resistance. *Science* **2002**, *298*, 96–97. [CrossRef]

4. Carvalho, R.A.; Omoto, C.; Field, L.M.; Williamson, M.S.; Bass, C. Investigating the Molecular Mechanisms of Organophosphate and Pyrethroid Resistance in the Fall Armyworm *Spodoptera frugiperda*. *PLoS ONE* **2013**, *8*, e62268. [CrossRef]

5. Shan, T.; Chen, C.; Ding, Q.; Chen, X.; Zhang, H.; Chen, A.; Shi, X.; Gao, X. Molecular Characterization and Expression Profiles of Nicotinic Acetylcholine Receptors in *Bradysia odoriphaga*. *Pestic. Biochem. Physiol.* **2020**, *165*, 104563. [CrossRef] [PubMed]

6. Kim, A.-Y.; Kwon, D.H.; Jeong, I.H.; Koh, Y.H. An Investigation of the Molecular and Biochemical Basis Underlying Chlorantraniliprole-Resistant *Drosophila* Strains and Their Cross-Resistance to Other Insecticides. *Arch. Insect Biochem. Physiol.* **2018**, *99*, e21514. [CrossRef]

7. Lin, L.; Hao, Z.; Cao, P.; Yuchi, Z. Homology Modeling and Docking Study of Diamondback Moth Ryanodine Receptor Reveals the Mechanisms for Channel Activation, Insecticide Binding and Resistance. *Pest Manag. Sci.* **2020**, *76*, 1291–1303. [CrossRef]

8. Montezano, D.G.; Specht, A.; Sosa-Gómez, D.R.; Roque-Specht, V.F.; Sousa-Silva, J.C.; de Paula-Moraes, S.V.; Peterson, J.A.; Hunt, T.E. Host Plants of *Spodoptera frugiperda* (Lepidoptera: Noctuidae) in the Americas. *Afr. Entomol.* **2018**, *26*, 286–300. [CrossRef]

9. Gouin, A.; Bretaudieu, A.; Nam, K.; Gimenez, S.; Aury, J.-M.; 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*, 11816. [CrossRef]

10. Gui, F.; Lan, T.; Zhao, Y.; Guo, W.; Dong, Y.; Fang, D.; Liu, H.; Li, H.; Wang, H.; Hao, R. Genomic and Transcriptomic Analysis Unveils Population Evolution and Development of Pesticide Resistance in Fall Armyworm *Spodoptera frugiperda*. *Protein Cell* **2020**. [CrossRef]

11. Xiao, H.; Ye, X.; Xu, H.; Mei, Y.; Yang, Y.; Chen, X.; Yang, Y.; Liu, T.; Yu, Y.; Yang, W. The Genetic Adaptations of Fall Armyworm *Spodoptera frugiperda* Facilitated Its Rapid Global Dispersal and Invasion. *Mol. Ecol. Resour.* **2020**, *20*, 1050–1068. [CrossRef] [PubMed]

12. Jakka, S.R.K.; Gong, L.; Hasler, J.; Banerjee, R.; Sheets, J.J.; Narva, K.; Blanco, C.A.; Jurat-Fuentes, J.L. Field-Evolved Mode 1 Resistance of the Fall Armyworm to Transgenic Cry1Fa-Expressing Corn Associated with Reduced Cry1Fa Toxin Binding and Midgut Alkaline Phosphatase Expression. *Appl. Environ. Microbiol.* **2016**, *82*, 1023–1034. [CrossRef] [PubMed]

13. Storer, N.P.; Babcock, J.M.; Schlenz, M.; Meade, T.; Thompson, G.D.; Bing, J.W.; Huckaba, R.M. Discovery and Characterization of Field Resistance to Bt Maize: *Spodoptera frugiperda* (Lepidoptera: Noctuidae) in Puerto Rico. *J. Econ. Entomol.* **2010**, *103*, 1031–1038. [CrossRef]

14. Monnerat, R.; Martins, E.; Macedo, C.; Queiroz, P.; Praça, L.; Soares, C.M.; Moreira, H.; Grisi, I.; Silva, J.; Soberon, M.; et al. Evidence of Field-Evolved Resistance of *Spodoptera frugiperda* to Bt Corn Expressing Cry1F in Brazil That Is Still Sensitive to Modified Bt Toxins. *PLoS ONE* **2015**, *10*, e0119544. [CrossRef]
15. Omoto, C.; Bernardi, O.; Salmeron, E.; Sorgatto, R.J.; Dourado, P.M.; Crivellari, A.; Carvalho, R.A.; Willse, A.; Martinelli, S.; Head, G.P. Field-Evolved Resistance to Cry1Ab Maize by *Spodoptera frugiperda* in Brazil. *Pest Manag. Sci.* **2016**, *72*, 1727–1736. [CrossRef] [PubMed]

16. Chandrasena, D.I.; Signorini, A.M.; Abratti, G.; Storer, N.P.; Olaciregui, M.L.; Alves, A.P.; Pilcher, C.D. Characterization of Field-Evolved Resistance to *Bacillus Thuringiensis*-Derived Cry1F δ-Endotoxin in *Spodoptera frugiperda* Populations from Argentina. *Pest Manag. Sci.* **2018**, *74*, 746–754. [CrossRef]

17. Schlum, K.A.; Lamour, K.; de Bortoli, C.P.; Banerjee, R.; Meagher, R.; Pereira, E.; Murua, M.G.; Sword, G.A.; Tessnow, A.E.; Viteri Dillon, D.; et al. Whole Genome Comparisons Reveal Panmixia among Fall Armyworm (*Spodoptera frugiperda*) from Diverse Locations. *BMC Genom.* **2021**, *22*, 179. [CrossRef]

18. Banerjee, R.; Hasler, J.; Meagher, R.; Nagoshi, R.; Hietala, L.; Huang, F.; Narva, K.; Jurat-Fuentes, J.L. Mechanism and DNA-Based Detection of Field-Evolved Resistance to Transgenic Bt Corn in Fall Armyworm (*Spodoptera frugiperda*). *Sci. Rep.* **2017**, *7*, 10877. [CrossRef]

19. Flagel, L.; Lee, Y.W.; Wanjugi, H.; Swarup, S.; Brown, A.; Wang, J.; Kraft, E.; Greenplate, J.; Simmons, J.; Adams, N.; et al. Mutational Disruption of the ABCC2 Gene in Fall Armyworm, *Spodoptera frugiperda*, Confers Resistance to the Cry1Fa and Cry1A.105 Insecticidal Proteins. *Sci. Rep.* **2018**, *8*, 7255. [CrossRef]

20. Boaventura, D.; Ulrich, J.; Lueke, B.; Bolzan, A.; Okuma, D.; Gutbrod, O.; Geibel, S.; Zeng, Q.; Dourado, P.M.; Martinelli, S.; et al. Molecular Characterization of Cry1F Resistance in Fall Armyworm, *Spodoptera frugiperda* from Brazil. *Insect Biochem. Mol. Biol.* **2020**, *116*, 103280. [CrossRef]

21. Guan, F.; Zhang, J.; Shen, H.; Wang, X.; Padovan, A.; Walsh, T.K.; Tay, W.T.; Gordon, K.H.J.; James, W.; Czepak, C.; et al. Whole-Genome Sequencing to Detect Mutations Associated with Resistance to Insecticides and Bt Proteins in *Spodoptera frugiperda*. *Insect Sci.* **2020**. [CrossRef] [PubMed]

22. Gutiérrez-Moreno, R.; Mota-Sanchez, D.; Blanco, C.A.; Whalon, M.E.; Terán-Santofimio, H.; Rodriguez-Maciel, J.C.; DiFonzo, C. Field-Evolved Resistance of the Fall Armyworm (Lepidoptera: Noctuidae) to Synthetic Insecticides in Puerto Rico and Mexico. *J. Econ. Entomol.* **2019**, *112*, 792–802. [CrossRef] [PubMed]

23. Gimenez, S.; Abdelgaffar, H.; Goff, G.L.; Hilliou, F.; Blanco, C.A.; Hänniger, S.; Bretaudieu, A.; Legeai, F.; Nègre, N.; Jurat-Fuentes, J.L.; et al. Adaptation by Copy Number Variation Increases Insecticide Resistance in the Fall Armyworm. *Commun. Biol.* **2020**, *3*, 664. [CrossRef] [PubMed]

24. Boaventura, D.; Bolzan, A.; Padovez, F.E.; Okuma, D.M.; Omoto, C.; Nauen, R. Detection of a Ryanodine Receptor Target-Site Mutation in Diamide Insecticide Resistant Fall Armyworm, *Spodoptera frugiperda*. *Pest Manag. Sci.* **2020**, *76*, 47–54. [CrossRef]

25. Goergen, G.; Kumar, P.L.; Sankung, S.B.; Togola, A.; Tamò, M. First Report of Outbreaks of the Fall Armyworm *Spodoptera frugiperda* (J E Smith) (Lepidoptera: Noctuidae), a New Alien Invasive Pest in West and Central Africa. *PLoS ONE* **2016**, *11*, e0165632. [CrossRef]

26. Day, R.; Abrahams, P.; Bateman, M.; Beale, T.; Clottey, V.; Cock, M.; Colmenarez, Y.; Corniani, N.; Early, R.; Godwin, J. Fall Armyworm: Impacts and Implications for Africa. *Outlooks Pest Manag.* **2017**, *28*, 196–201. [CrossRef]

27. Pashley, D.P. Host-Associated Genetic Differentiation in Fall Armyworm (Lepidoptera: Noctuidae): A Sibling Species Complex? *Ann. Entomol. Soc. Am.* **1986**, *79*, 898–904. [CrossRef]

28. Pashley, D.P. Host-associated differentiation in armyworms (Lepidoptera: Noctuidae): An allozymic and mitochondrial DNA perspective. In *Electrophoretic Studies on Agricultural Pests*; Loxdale, H.D., Hollander, J.D., Eds.; Clarendon Press: Oxford, UK, 1989; Volume 39, pp. 103–144.

29. Dumas, P.; Legeai, F.; Lemaitre, C.; Scaon, E.; Orsucci, M.; Labadie, K.; Gimenez, S.; Clamens, A.-L.; Henri, H.; Vavre, F.; et al. *Spodoptera frugiperda* (Lepidoptera: Noctuidae) Host-Plant Variants: Two Host Strains or Two Distinct Species? *Genetica* **2015**, *143*, 305–316. [CrossRef]

30. Nagoshi, R.N.; Goergen, G.; Plessis, H.D.; van den Berg, J.; Meagher, R. Genetic Comparisons of Fall Armyworm Populations from 11 Countries Spanning Sub-Saharan Africa Provide Insights into Strain Composition and Migratory Behaviors. *Sci. Rep.* **2019**, *9*, 8311. [CrossRef]

31. Yainna, S.; Tay, W.T.; Fiteni, E.; Legeai, F.; Clamens, A.-L.; Gimenez, S.; Frayssinet, M.; Asokan, R.; Kalleshwaraswamy, C.M.; Deshmukh, S.; et al. Genomic Balancing Selection Is Key to the Invasive Success of the Fall Armyworm. *bioRxiv* **2020**. [CrossRef]

32. Zhang, L.; Liu, B.; Zheng, W.; Liu, C.; Zhang, D.; Zhao, S.; Li, Z.; Xu, P.; Wilson, K.; Withers, A.; et al. Genetic Structure and Insecticide Resistance Characteristics of Fall Armyworm Populations Invading China. *Mol. Ecol. Resour.* **2020**, *20*, 1682–1696. [CrossRef] [PubMed]

33. Tay, W.T.; Rane, R.; Padovan, A.; Walsh, T.; Elfekih, S.; Downes, S.; Nam, K.; d'Alençon, E.; Zhang, J.; Wu, Y.; et al. Whole Genome Sequencing of Global *Spodoptera frugiperda* Populations: Evidence for Complex, Multiple Introductions across the Old World. *bioRxiv* **2020**. [CrossRef]

34. Boaventura, D.; Martin, M.; Pozzebon, A.; Mota-Sanchez, D.; Nauen, R. Monitoring of Target-Site Mutations Conferring Insecticide Resistance in *Spodoptera frugiperda*. *Insects* **2020**, *11*, 545. [CrossRef]

35. McKenna, A.; Hanna, M.; Banks, E.; Sivachenko, A.; Cibulskis, K.; Kernytsky, A.; Garimella, K.; Altshuler, D.; Gabriel, S.; Daly, M.; et al. The Genome Analysis Toolkit: A MapReduce Framework for Analyzing next-Generation DNA Sequencing Data. *Genome Res.* **2010**, *20*, 1297–1303. [CrossRef] [PubMed]
36. Wang, X.; Zheng, Z.; Cai, Y.; Chen, T.; Li, C.; Fu, W.; Jiang, Y. CNVcaller: Highly Efficient and Widely Applicable Software for Detecting Copy Number Variations in Large Populations. *GigaScience* **2017**, *6*. [CrossRef]

37. Create Your Own Custom Map. Available online: https://mapchart.net/index.html (accessed on 17 May 2021).

38. Slater, G.S.C.; Birney, E. Automated Generation of Heuristics for Biological Sequence Comparison. BMC Bioinform. **2005**, *6*, 31. [CrossRef] [PubMed]

39. Robinson, J.T.; Thorvaldsdóttir, H.; Winckler, W.; Guttman, M.; Lander, E.S.; Getz, G.; Mesirov, J.P. Integrative Genomics Viewer. Nat. Biotechnol. **2011**, *29*, 24–26. [CrossRef]

40. Camacho, C.; Coulouris, G.; Avagyan, V.; Ma, N.; Papadopoulos, J.; Bealer, K.; Madden, T.L. BLAST+: Architecture and Applications. BMC Bioinform. **2009**, *10*, 421. [CrossRef]

41. Weir, B.S.; Cockerham, C.C. Estimating F-Statistics for the Analysis of Population Structure. *Evolution* **1984**, *38*, 1358–1370.

42. Danecek, P.; Auton, A.; Abecasis, G.; Albers, C.A.; Banks, E.; DePristo, M.A.; Handsaker, R.E.; Lunter, G.; Marth, G.T.; Sherry, S.T.; et al. The Variant Call Format and VCFtools. *Bioinformatics* **2011**, *27*, 2156–2158. [CrossRef]

43. Walsh, T.K.; Joussen, N.; Tian, K.; McGaughran, A.; Anderson, C.J.; Qiu, X.; Ahn, S.-J.; Bird, L.; Pavlidi, N.; Vontas, J.; et al. Multiple Recombination Events between Two Cytochrome P450 Loci Contribute to Global Pyrethroid Resistance in *Helicoverpa armigera*. *PLoS ONE* **2018**, *13*, e0197760. [CrossRef] [PubMed]

44. Valencia-Montoya, W.A.; Elfekih, S.; North, H.L.; Meier, J.I.; Warren, I.A.; Tay, W.T.; Gordon, K.H.; Specht, A.; Paula-Moraes, S.V.; Rane, R. Adaptive Introgression across Semipermeable Species Boundaries between Local *Helicoverpa zea* and Invasive *Helicoverpa armigera* Moths. *Mol. Biol. Evol.* **2020**, *37*, 2568–2583. [CrossRef] [PubMed]

45. Zhang, D.; Xiao, Y.; Xu, P.; Yang, X.; Wu, Q.; Wu, K. Insecticide Resistance Monitoring for the Invasive Populations of Fall Armyworm, *Spodoptera frugiperda* in China. *J. Integr. Agric.* **2021**, *20*, 783–791. [CrossRef]

46. Zhang, W. Global Pesticide Use: Profile, Trend, Cost/Benefit and More. Proc. Int. Acad. Ecol. Environ. Sci. **2018**, *8*, 1–27.

47. Lanfear, R.; Kokko, H.; Eyre-Walker, A. Population Size and the Rate of Evolution. *Trends Ecol. Evol.* **2014**, *29*, 33–41. [CrossRef] [PubMed]

48. Nam, K.; Munch, K.; Mailund, T.; Nater, A.; Greminger, M.P.; Krützen, M.; Marquès-Bonet, T.; Schierup, M.H. Evidence That the Rate of Strong Selective Sweeps Increases with Population Size in the Great Apes. Proc. Natl. Acad. Sci. USA **2017**, *114*, 1613–1618. [CrossRef] [PubMed]

49. Nam, K.; Nhim, S.; Robin, S.; Bretaudieu, A.; Nègre, N.; d'Alençon, E. Positive Selection Alone Is Sufficient for Whole Genome Differentiation at the Early Stage of Speciation Process in the Fall Armyworm. BMC Evol. Biol. **2020**, *20*, 152. [CrossRef]
<figure><img src="image_68.png" /></figure>

# Evaluation of Reference Genes and Expression Level of Genes Potentially Involved in the Mode of Action of Cry1Ac and Cry1F in a Susceptible Reference Strain of *Chrysodeixis includens*

Macarena Martin<fnref n="1" />, Debora Boaventura<fnref n="2" /> and Ralf Nauen<fnref n="2" />,*

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

check for
updates

**Citation:** Martin, M.; Boaventura, D.; Nauen, R. Evaluation of Reference Genes and Expression Level of Genes Potentially Involved in the Mode of Action of Cry1Ac and Cry1F in a Susceptible Reference Strain of *Chrysodeixis includens*. *Insects* **2021**, *12*, 598. https://doi.org/10.3390/insects12070598

Academic Editor: Yulin Gao

Received: 2 June 2021

Accepted: 29 June 2021

Published: 30 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Simple Summary:** Soybean looper (a moth species) is a major pest of soybean plants in the American continent and its larvae need to be kept under economic damage thresholds to guarantee sustainable yields. Soybean looper control relies mostly on the use of insecticides and genetically modified crops expressing *Bacillus thuringiensis* (Bt) insecticidal proteins. Due to the high selection pressure exerted by these control measures, resistance has developed to different insecticides and Bt proteins. Here, we tested the basal sensitivity of a soybean looper laboratory reference strain against two insecticidal proteins and determined the level of expression of potential receptors of these proteins across all (six) larval stages. Furthermore, we identified stable reference genes across all larval stages to normalize gene expression data obtained by quantitative polymerase chain reaction (qPCR). The results presented in this communication are useful to support future studies on insecticide and insecticidal protein resistance in soybean looper.

**Abstract:** Soybean looper (SBL), *Chrysodeixis includens* (Walker), is one of the major lepidopteran pests of soybean in the American continent. SBL control relies mostly on the use of insecticides and genetically modified crops expressing *Bacillus thuringiensis* (Bt) insecticidal Cry proteins. Due to the high selection pressure exerted by these control measures, resistance has developed to different insecticides and Bt proteins. Nevertheless, studies on the mechanistic background are still scarce. Here, the susceptibility of the laboratory SBL-Benzon strain to the Bt proteins Cry1Ac and Cry1F was determined in diet overlay assays and revealed a greater activity of Cry1Ac than Cry1F, thus confirming results obtained for other sensitive SBL strains. A reference gene study across larval stages with four candidate genes revealed that *RPL10* and *EF1* were the most stable genes for normalization of gene expression data obtained by RT-qPCR. Finally, the basal expression levels of eight potential Bt protein receptor genes in six larval instars were analyzed, including ATP-binding cassette (ABC) transporters, alkaline phosphatase, aminopeptidases, and cadherin. The results presented here provide fundamental knowledge to support future SBL resistance studies.

**Keywords:** soybean looper; *Bacillus thuringiensis*; resistance; reference genes; ABCC2

## 1. Introduction

Soybean looper (SBL), *Chrysodeixis includens* (Walker) (Lepidoptera: Noctuidae), is a distinctive species from the Western hemisphere and a key pest of soybean in many countries of the American continent [1]. The control of SBL relies on genetically modified crops expressing *Bacillus thuringiensis* (Bt) insecticidal proteins and foliar insecticide applications [2].

<footnote n="1">¹ Department of Agronomy, Food, Natural Resources, Animals and Environment, University of Padova, 35020 Padova, Italy; macamar_20@hotmail.com</footnote>

<footnote n="2">² Bayer AG, Crop Science Division, R&D Pest Control, Alfred Nobel Str. 50, 40789 Monheim, Germany; debora.duarteboaventura@bayer.com</footnote>

<footnote>* Correspondence: ralf.nauen@bayer.com</footnote>
Frequent insecticide applications have led to insecticide resistance against some major chemical classes, such as pyrethroids, carbamates, organophosphates, benzoylureas, and diamides [3–5]. Moreover, the high selection pressure exerted by millions of hectares of soybean expressing Bt proteins has resulted in resistance towards Cry1Ac [6].

A reduction in Cry protein binding due to changes in the expression and/or muta-
tions in midgut receptors involved in Cry protein-mediated intoxication was the most
frequent mechanism of resistance [7]. Different proteins have been identified as receptors
for Cry proteins in lepidopteran species, including aminopeptidases (APN), membrane-
bound alkaline phosphatases (mALP), cadherins (CAD), and ATP-binding cassette (ABC)
transporters [8].

Brazilian field populations of SBL selected against Cry1Ac in the laboratory showed up to 127-fold resistance [6]. However, the molecular mechanism involved in the resistance of SBL towards Cry1Ac is still unknown. In other lepidopteran pest species, resistance to Cry1Ac has been associated with mutations and/or differences in the expression levels of ABCC2, APN, mALP, and CAD [9–13].

Here, we have tested the basal sensitivity of the commonly known reference SBL-Benzon strain towards two major Bt proteins, Cry1Ac and Cry1F. Moreover, we have investigated, for the first time, the stability of potential reference genes with low expression variance across larval instars of SBL for normalization of RT-qPCR expression data. As resistance towards Bt proteins in other species has been linked to differences in the expression of potential Bt receptor genes, we have analyzed the basal expression of eight potential Bt protein receptor genes in six larval instars to support the characterization of resistance mechanisms to Cry toxins in SBL in future studies.

**2. Material and Methods**

**2.1. Insects**

SBL-Benzon strain from Benzon Research Inc., USA was reared in the laboratory under controlled conditions (25 ± 1°C, 55 ± 5% relative humidity, and a photoperiod of L16:D8). Briefly, adults were kept in cages covered with filter paper for egg laying and fed with honey water pads every second day. Larvae were mass-kept in a plastic box (18 × 13 × 5 cm) containing Soybean-Looper-Diet (Southland Products Inc., Lake Village, AR, USA) until the third instar. Afterwards, two third-instar larvae were confined to petri dishes (diameter 6 cm) containing the same diet until pupation.

2.2. *Cry1Ac and Cry1F bioassays*

The bioassays were performed according to Boaventura et al. 2020 [14]. Briefly, seven concentrations of Cry1Ac (450-0.62 ng cm⁻²) and eight concentrations of Cry1F (12,150-5.56 ng cm⁻²) were tested. Dilutions were prepared from Cry1Ac lyophilized material (28.2% (w/w) purity) and Cry1F (TIC-842) crystal spore solution [15] in 50 mM sodium carbonate buffer (pH 10.4) and 0.1% (v/v) Triton X-100. Twenty-five µL of Bt solution was applied to the surface of artificial diet in a 48-well plate (Greiner CELLSTAR®, Merck, Darmstadt, Germany) and a single SBL neonate larva (<24 h old) was placed into each well and kept at 20 ± 2 °C.

A total number of 12 neonates were used per protein concentration and control (buffer and 0.1% (v/v) Triton X-100). For the highest concentration, only six neonates were used and the other six neonates were placed in wells containing only artificial diet. The bioassay was replicated three times and the mortality was scored after five days; larvae that did not move or reach the third instar were considered dead. The results were fitted by a logistic regression model to calculate LC₅₀/LC₉₅ values and 95% confidence intervals (GraphPad Software Inc. v.5, San Diego, CA, USA).

2.3. RNA Extraction and cDNA Synthesis

Total RNA from each instar was extracted from 30 pooled (first instar) or 10 pooled
larvae (second to sixth instar) of SBL in four biological replicates. The RNA extraction
from first and second instars was performed with the PicoPure™ RNA Isolation Kit (Thermo Fisher, Vilnius, Lithuania) according to the manufacturer's instructions. For the third to sixth instar, the larval tissue was ground with a Retsch TissueLyser (QIAGEN, Hilden, Germany) or mortar and pestle. The total RNA was extracted from powdered samples using TRIzol® reagent (Invitrogen, Waltham, MA, USA) followed by RNA purification according to RNeasy® Plus Universal Mini Kit (QIAGEN, Hilden, Germany) recommendations including gDNA elimination with DNase I. The RNA concentration was normalized to 100 ng µL⁻¹ and 1 µg total RNA was used in 20 µL reactions for cDNA synthesis using iScript™ cDNA synthesis (Bio-Rad, Laboratories, CA, USA) according to the manufacturer's instructions.

## 2.4. Screening for Reference Genes

The expression of four candidate reference genes was tested: ribosomal protein L10 (RPL10), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), elongation factor-1 alpha (EF1), and beta actin 1 (ACT1). The reference gene sequences were obtained from a custom-made SBL transcriptome generated at Bayer AG through queries using local BLAST searches and submitted to GenBank (Table S1). The primers were designed using Geneious software v. 10.2.3 (Biomatters Ltd, Auckland, New Zealand) (Table S1). The RT-qPCR reactions were performed as recently described [14] using reverse/forward primers shown in Table S1. Amplification efficiencies were determined by a series of 5-fold dilutions to create standard curves by a linear regression model [16] (Table S1).

The cycle threshold values (Ct values) were obtained from CFX Maestro 1.0 v4.0 (Bio-Rad) software. The expression stability of four candidate reference genes was analyzed using geNorm qBase Plus v3.1 software (Biogazelle, Belgium) and NormFinder (moma.dk/normfinder-software) using stability value (M) below 0.5 and pairwise variation (V) value below 0.15 [16]. The Excel-based tool Normfinder was used according to Andersen et al. (2004) [17].

## 2.5. Expression Analysis of Possible Bt Protein Target Genes

The expression level of eight potential Bt receptor genes (ABCC2, ABCC3, CAD, mALP, and APN (APN1-4)) was determined by RT-qPCR in six larval instars. The genes EF1 and RPL10 were the most stable genes across SBL larval instars and were used for normalization. In total, four biological replicates containing cDNA from a pool of 10-30 larvae, as described in Section 2.3, were used for analyzing the normalized expression levels. Primer pairs are described in Table S1 and the RT-qPCR conditions were chosen as recently described [14]. The normalized expression values obtained for the different developmental stages were compared to the Ct values from the first instar and statistically tested by one-way ANOVA, followed by Tukey-Kramer multiple-comparison test with qbase + software (v. 3.2; Biogazelle, Belgium).

# 3. Results

## 3.1. Bioassays

For Cry1Ac, the LC₅₀ and LC₉₅ values were 10.41 ng cm⁻² and 44.56 ng cm⁻², respectively. Cry1F showed an LC₅₀ of 228.1 ng cm⁻² and an LC₉₅ of 4,967 ng cm⁻² (Figure 1, Table S2).

## 3.2. Screening for Reference Genes

According to geNorm, RPL10 was identified as the most stable gene across larval stages of SBL (M = 0.458), followed by EF1, GAPDH, and ACT1 with M values of 0.519, 0.559, and 0.637, respectively (Table 1). When using NormFinder, EF1 (stability value = 0.242) and RPL10 (0.244) were the most stable genes, followed by GAPDH (0.309) and ACT1 (0.351) (Table 1).
<figure><img src="image_71.png" /><figcaption>**Figure 1.** Toxicity of Bt proteins against neonates of soybean looper strain SBL-Benzon. (A) Toxicity of Cry1Ac (N = 288) and (B) toxicity of Cry1F (N = 324). Dashed lines indicate the 95% confidence band.</figcaption></figure>

<table><caption>**Table 1.** Stability analysis of the candidate reference genes *ribosomal protein L10* (RPL10), *glyceraldehyde-3-phosphate dehydrogenase* (GAPDH), *elongation factor-1 alpha* (EF1), and *beta actin 1* (ACT1) analyzed by geNorm (average expression stability or M value) and NormFinder (stability value) in different larval stages of *Chrysodeixis includens*.</caption><thead><tr><th>Method</th><th>EF1</th><th>RPL10</th><th>GAPDH</th><th>ACT1</th></tr></thead><tbody><tr><td>geNorm</td><td>0.519</td><td>0.458</td><td>0.559</td><td>0.637</td></tr><tr><td>NormFinder</td><td>0.242</td><td>0.244</td><td>0.309</td><td>0.351</td></tr></tbody></table>

### 3.3. Expression Analysis of Potential Bt Protein Target Genes

The basal expression levels of the genes *ABCC2*, *ABCC3*, *CAD*, *mALP*, and *APN* (*APN1-4*), reported to be potential receptors of Cry1Ac and Cry1F in other lepidopteran species, were analyzed in the different larval stages of the Bt-susceptible SBL-Benzon strain. For *ABCC3*, no significant difference in the expression level across larval stages was observed, with average normalized values between 1 and 2.02. In contrast, the expression level of *ABCC2* was significantly lower in second-instar larvae (*p*-value < 0.05) (Figure 2).

<figure><img src="image_72.png" /><figcaption>**Figure 2.** *Cont.*</figcaption></figure>
<figure><img src="image_73.png" /><figcaption>**Figure 2.** RT-qPCR expression analysis of genes possibly involved in Bt toxicity, analyzed by geNorm across different larval stages of *Chrysodeixis includens*. Average normalized quantities of (A) ABCC2, (B) ABCC3, (C) APN1, (D) APN2, (E) APN3, (F) APN4, (G) ALP, and (H) CAD. CAD: cadherin; ABCC2/ABCC3: ATP-binding cassette subfamily C2/C3 transporter; ALP: alkaline phosphatase; APN: amino peptidase class 1, 2, 3, and 4. Different letters denote a significant difference (One-way ANOVA; post hoc Tukey comparison, *p* < 0.05).</figcaption></figure>

The expression level of the individual *APN* genes was not significantly different between the different larval stages. Furthermore, no significant differences in the expression level of *mALP* were detected. However, the expression level of *CAD* in second-instar larvae was significantly lower than in fifth-instar larvae (Figure 2).

**4. Discussion**

The use of insecticide-susceptible reference strains such as SBL-Benzon is a key aspect in resistance monitoring and mechanism studies. SBL is a notorious pest of soybean crops, but it also feeds on cotton [18]. Cry1Ac is utilized in both crops and a relatively high survival rate of SBL in Cry1Ac Bt cotton has been reported [19]. The primary Insecticide Resistance Management (IRM) strategy for Bt crops has been the use of refuges (presence of non-Bt host plants near Bt crops). The success of refuge-based IRM depends on a low frequency of potential resistance alleles, high susceptibility to the expressed Bt protein(s), and a recessive mode of inheritance [20]. Utilizing Bt crops expressing more than one Bt protein (pyramid strategy), such as Bt soybean expressing Cry1Ac and Cry1F, is another valuable IRM strategy [21]. The combination of Cry1F and Cry1Ac in transgenic cotton was shown to be highly effective against SBL [22]. However, soybean and cotton are usually planted in succession in Brazil, thus increasing the resistance risk of SBL by continuous selection pressure, particularly against Cry1F, which exhibited significantly lower activity than Cry1Ac against the susceptible reference strain SBL-Benzon. The data obtained here were close to dose-response data published recently and showing the superior efficacy of Cry1Ac compared to Cry1F against neonates of a susceptible Brazilian SBL field strain [6]. This is in contrast to toxicity data previously published in another study, where Cry1F was more active than Cry1Ac against neonates of a sensitive SBL laboratory colony, whereas brush-border membrane vesicle (BBMV) binding studies revealed the opposite, i.e., a higher affinity of Cry1Ac than Cry1F [23].
The mechanisms driving Bt resistance are quite diverse and often associated with the differential expression of genes of potential Bt receptor proteins, particularly in noctuid pests reviewed in [24]. Bt protein resistance in SBL is not yet widespread, and mechanistic studies of factors driving resistance in SBL are scarce. However, the elucidation of mechanisms conferring Bt resistance based on differential gene expression analysis requires stable reference genes for normalization of RT-qPCR data [25]. Such studies have not yet been conducted in SBL, so we investigated the stability of four candidate reference genes across six larval instars to allow the selection of at least two appropriate reference genes typically considered sufficient for the normalization of RT-qPCR data to reduce any systematic bias [26]. Based on the results obtained in our study, *RPL10* and *EF1* were most stably expressed across larval stages of SBL and deemed suitable to be used in future gene expression studies. However, *EF1* was recently already employed for the normalization of the expression level of detoxification genes in a pyrethroid-resistant SBL strain [27].

Finally, we utilized the identified reference genes to normalize RT-qPCR data obtained for the basal expression level of potential Bt receptor genes in strain SBL-Benzon across six larval instars. Such data are a useful resource for the comparative analysis of gene expression levels in future studies with SBL strains resistant to Bt proteins. *ABCC2* has been described as an important receptor for Cry1Ac [10] and Cry1F [28] in lepidopteran pests, and mutations in *ABCC2* were shown to confer resistance to Cry1 proteins in *H. armigera* [10] and other noctuid pests such as fall armyworm [14,15,24]. Meanwhile, in *Plutella xylostella*, Cry1Ac resistance was linked to the downregulation of *ABCC2*, *ABCC3*, and *mALP* [13]. However, in a *Helicoverpa zea* strain selected with Cry1Ac, resistance was associated with an increase in the expression of *ALP* [29]. Target-site mutations in *APN1* and differences in its expression level have been reported to cause resistance to Cry1Ac in *H. armigera* and *Trichoplusia ni*, respectively [9,11]. Cadherins were also shown to be receptors of Cry1Ac in different lepidopteran species [24,30]. Differences in cadherin levels and/or mutations in this gene were, for example, associated with Cry1Ac resistance in *Ostrinia furnacalis* [12]. Our results reveled no differences in *APN* and *ALP* expression levels across larval instars of Bt-susceptible SBL, but significant differences in *CAD* expression levels between the second and fifth larval instar; however, additional work will be necessary to investigate if such differences have implications for Bt protein toxicity.

In summary, we have shown the stability of the genes *RPL10* and *EF1* across six different larval stages of SBL and recommend their use for the normalization of gene expression data in future studies. In addition, we identified and described potential Bt target genes in SBL and investigated their basal expression profiles across larval instars in the reference SBL-Benzon strain to support future resistance studies in SBL.

**5. Conclusions**

Resistance towards Bt proteins has already been reported in SBL and the results presented in this study will help to support future research on SBL resistance towards both Bt proteins and synthetic insecticides. Besides the data presented on the susceptibility towards Cry1Ac and Cry1F proteins, this study provides useful insights into the stability of reference genes, with *RPL10* and *EF1* as the most stable ones. Moreover, for the first time, potential Bt target genes (*ABCC2*, *ABCC3*, *cadherin*, *mALP*, and several *APNs*) and their respective expression profiles were investigated for a susceptible SBL reference strain, SBL-Benzon.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10.3390/insects12070598/s1, Table S1: Gene names of orthologous candidate reference genes (*RPL10*, *GADPH*, *EF1*, and *ACT1*), potential orthologous Cry toxin target genes, species, GenBank accession numbers and primer sequences for the amplification of the respective Chrysodeixis includens genes, Table S2: Toxicity (LC50-and LC95-values) of Cry1Ac and Cry1F to neonates of soybean looper Benzon strain and the respective goodness of fit calculated for the non-linear regressions used.
**Author Contributions:** R.N. and D.B. conceived the study. R.N. supervised the project. D.B. and M.M. performed the research and analyzed data. D.B., M.M. and R.N. wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** M.M. thanks the Erasmus+—KA1 Erasmus Mundus Joint Master’s Degree Programme of the European Commission as this paper has been developed as a result of a mobility stay funded by the PLANT HEALTH Project.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in the supplementary materials or on request from the corresponding author.

**Conflicts of Interest:** R.N. and D.B. are employed by Bayer AG, a manufacturer of insecticides and transgenic soybean. This work was partially financed by Bayer AG. The authors declare no additional conflict of interest.

## References

1. Herzog, D.C. Sampling Soybean Looper on Soybean. In *Sampling Methods in Soybean Entomology*; Kogan, M., Herzog, D.C., Eds.; Springer: New York, NY, USA, 1980; pp. 141–168. [CrossRef]

2. Di Oliveira, J.R.G.; Ferreira, M.D.C.; Román, R.A.A. Diferentes diâmetros de gotas e equipamentos para aplicação de inseticida no controle de Pseudoplosia includens. *Eng. Agrícola* **2010**, *30*, 92–99. [CrossRef]

3. Boernel, D.J.; Mink, J.S.; Wier, A.T.; Thomas, J.D.; Leonard, B.R.; Gallardo, F. Management of Insecticide Resistant Soybean Loopers (Pseudoplosia Includens) in the Southern United States. In *Pest Management in Soybean*; Copping, L.G., Green, M.B.M., Rees, R.T., Eds.; Springer: Dordrecht, The Netherlands, 1992; pp. 66–87.

4. Owen, L.N.; Catchot, A.; Musser, F.; Gore, J.; Cook, D.; Jackson, R. Susceptibility of Chrysodeixis includens (Lepidoptera: Noctuidae) to Reduced-Risk Insecticides. *Fla. Éntomol.* **2013**, *96*, 554–559. [CrossRef]

5. Stacke, R.F.; Godoy, D.; Pretto, V.E.; Führ, F.M.; Gubiani, P.D.S.; Hettwer, B.L.; Garlet, C.G.; Somavilla, J.C.; Muraro, D.S.; Bernardi, O. Field-evolved resistance to chitin synthesis inhibitor insecticides by soybean looper, Chrysodeixis includens (Lepidoptera: Noctuidae), in Brazil. *Chemosphere* **2020**, *259*, 127499. [CrossRef]

6. Rodrigues-Silva, N.; Canuto, A.F.; Oliveira, D.F.; Teixeira, A.F.; Amaya, O.F.S.; Picanço, M.C.; Pereira, E.J.G. Negative cross-resistance between structurally different Bacillus thuringiensis toxins may favor resistance management of soybean looper in transgenic Bt cultivars. *Sci. Rep.* **2019**, *9*, 199. [CrossRef] [PubMed]

7. Heckel, D.G.; Gahan, L.J.; Baxter, S.W.; Zhao, J.; Shelton, A.M.; Gould, F.; Tabashnik, B.E. The diversity of Bt resistance genes in species of Lepidoptera. *J. Invertebr. Pathol.* **2007**, *95*, 192–197. [CrossRef]

8. Bravo, A.; Gill, S.S.; Soberón, M. Mode of action of Bacillus thuringiensis Cry and Cyt toxins and their potential for insect control. *Toxicon* **2007**, *49*, 423–435. [CrossRef] [PubMed]

9. Tiewsiri, K.; Wang, P. Differential alteration of two aminopeptidases N associated with resistance to Bacillus thuringiensis toxin Cry1Ac in cabbage looper. *Proc. Natl. Acad. Sci. USA* **2011**, *108*, 14037–14042. [CrossRef] [PubMed]

10. Xiao, Y.; Zhang, T.; Liu, C.; Heckel, D.; Li, X.; Tabashnik, B.E.; Wu, K. Mis-splicing of the ABCC2 gene linked with Bt toxin resistance in Helicoverpa armigera. *Sci. Rep.* **2014**, *4*, 6184. [CrossRef]

11. Zhang, S.; Cheng, H.; Gao, Y.; Wang, G.; Liang, G.; Wu, K. Mutation of an aminopeptidase N gene is associated with Helicoverpa armigera resistance to Bacillus thuringiensis Cry1Ac toxin. *Insect Biochem. Mol. Biol.* **2009**, *39*, 421–429. [CrossRef]

12. Jin, T.; Chang, X.; Gatehouse, A.M.R.; Wang, Z.; Edwards, M.G.; He, K. Downregulation and Mutation of a Cadherin Gene Associated with Cry1Ac Resistance in the Asian Corn Borer, Ostrinia furnacalis (Guenée). *Toxins* **2014**, *6*, 2676–2693. [CrossRef]

13. Guo, Z.; Kang, S.; Chen, D.; Wu, Q.; Wang, S.; Xie, W.; Zhu, X.; Baxter, S.W.; Zhou, X.; Jurat-Fuentes, J.L.; et al. MAPK Signaling Pathway Alters Expression of Midgut ALP and ABCC Genes and Causes Resistance to Bacillus thuringiensis Cry1Ac Toxin in Diamondback Moth. *PLoS Genet.* **2015**, *11*, e1005124. [CrossRef] [PubMed]

14. Boaventura, D.; Ulrich, J.; Lueke, B.; Bolzan, A.; Okuma, D.; Gutbrod, O.; Geibel, S.; Zeng, Q.; Dourado, P.M.; Martinelli, S.; et al. Molecular characterization of Cry1F resistance in fall armyworm, Spodoptera frugiperda from Brazil. *Insect Biochem. Mol. Biol.* **2020**, *116*, 103280. [CrossRef]

15. Flagel, L.; Lee, Y.W.; Wanjugi, H.; Swarup, S.; Brown, A.; Wang, J.; Kraft, E.; Greenplate, J.; Simmons, J.; Adams, N.; et al. Mutational disruption of the ABCC2 gene in fall armyworm, Spodoptera frugiperda, confers resistance to the Cry1Fa and Cry1A.105 insecticidal proteins. *Sci. Rep.* **2018**, *8*, 1–11. [CrossRef]

16. Pfaffl, M.W.; Tichopad, A.; Prgomet, C.; Neuvians, T.P. Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper—Excel-based tool using pair-wise correlations. *Biotechnol. Lett.* **2004**, *26*, 509–515. [CrossRef]
17. Andersen, C.L.; Jensen, J.L.; Ørntoft, T.F. Normalization of Real-Time Quantitative Reverse Transcription-PCR Data: A Model-Based Variance Estimation Approach to Identify Genes Suited for Normalization, Applied to Bladder and Colon Cancer Data Sets. *Cancer Res.* **2004**, *64*, 5245–5250. [CrossRef]

18. Panizzi, A.R. History and Contemporary Perspectives of the Integrated Pest Management of Soybean in Brazil. *Neotrop. Entomol.* **2013**, *42*, 119–127. [CrossRef]

19. Funichello, M.; Fern, J.; Grigolli, o.J.; de Souza, B.H.S.; Boiccedil, A.L.; Junior, A.; Busoli, A.C. Effect of transgenic and non-transgenic cotton cultivars on the development and survival of Pseudoplusia includens (Walker) (Lepidoptera: Noctuidae). *AJAR 2013*, *8*, 5424–5428. [CrossRef]

20. Tabashnik, B.E.; Gould, F.; Carriere, Y. Delaying evolution of insect resistance to transgenic crops by decreasing dominance and heritability. *J. Evol. Biol.* **2004**, *17*, 904–912. [CrossRef] [PubMed]

21. Marques, L.H.; Castro, B.A.; Rossetto, J.; Silva, O.A.B.N.; Moscardini, V.F.; Zobiole, L.H.S.; Santos, A.C.; Valverde-Garcia, P.; Babcock, J.M.; Rule, D.M.; et al. Efficacy of Soybean's Event DAS-81419-2 Expressing Cry1F and Cry1Ac to Manage Key Tropical Lepidopteran Pests Under Field Conditions in Brazil. *J. Econ. Entomol.* **2016**, *109*, 1922–1928. [CrossRef] [PubMed]

22. Tindall, K.V.; Siebert, M.W.; Leonard, B.R.; All, J.; Haile, F.J. Efficacy of Cry1Ac:Cry1F proteins in cotton leaf tissue against fall armyworm, beet armyworm, and soybean looper (Lepidoptera: Noctuidae). *J. Econ. Entomol.* **2009**, *102*, 1497–1505. [CrossRef]

23. Bel, Y.; Sheets, J.J.; Tan, S.Y.; Narva, K.E.; Escriche, B. Toxicity and Binding Studies of Bacillus thuringiensis Cry1Ac, Cry1F, Cry1C, and Cry2A Proteins in the Soybean Pests Anticarsia gemmatalis and Chrysodeixis (Pseudoplusia) includens. *Appl. Environ. Microbiol.* **2017**, *83*, e00326-17. [CrossRef]

24. Jurat-Fuentes, J.L.; Heckel, D.G.; Ferré, J. Mechanisms of Resistance to Insecticidal Proteins from Bacillus thuringiensis. *Annu. Rev. Entomol.* **2021**, *66*, 121–140. [CrossRef] [PubMed]

25. Bustin, S.A.; Benes, V.; Garson, J.A.; Hellemans, J.; Huggett, J.; Kubista, M.; Mueller, R.; Nolan, T.; Pfaffl, M.W.; Shipley, G.L.; et al. The MIQE Guidelines: Minimum Information for Publication of Quantitative Real-Time PCR Experiments. *Clin. Chem.* **2009**, *55*, 611–622. [CrossRef] [PubMed]

26. Vandesompele, J.; De Preter, K.; Pattyn, F.; Poppe, B.; Van Roy, N.; De Paepe, A.; Speleman, F. Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. *Genome Biol.* **2002**, *3*, research0034.1. [CrossRef] [PubMed]

27. Perini, C.R.; Tabuloc, C.A.; Chiu, J.C.; Zalom, F.G.; Stacke, R.F.; Bernardi, O.; Nelson, D.R.; Guedes, J.C. Transcriptome Analysis of Pyrethroid-Resistant Chrysodeixis includens (Lepidoptera: Noctuidae) Reveals Overexpression of Metabolic Detoxification Genes. *J. Econ. Entomol.* **2021**, *114*, 274–283. [CrossRef] [PubMed]

28. Coates, B.S.; Siegfried, B.D. Linkage of an ABCC transporter to a single QTL that controls Ostrinia nubilalis larval resistance to the Bacillus thuringiensis Cry1Fa toxin. *Insect Biochem. Mol. Biol.* **2015**, *63*, 86–96. [CrossRef] [PubMed]

29. Caccia, S.; Moar, W.J.; Chandrashekhar, J.; Oppert, C.; Anilkumar, K.J.; Jurat-Fuentes, J.L.; Ferré, J. Association of Cry1Ac Toxin Resistance in Helicoverpa zea (Boddie) with Increased Alkaline Phosphatase Levels in the Midgut Lumen. *Appl. Environ. Microbiol.* **2012**, *78*, 5690–5698. [CrossRef]

30. Pigott, C.R.; Ellar, D.J. Role of Receptors in Bacillus thuringiensis Crystal Toxin Activity. *Microbiol. Mol. Biol. Rev.* **2007**, *71*, 255–281. [CrossRef]
<figure><img src="image_74.png" /></figure>

Article

# Genetic Screening to Identify Candidate Resistance Alleles to Cry1F Corn in Fall Armyworm Using Targeted Sequencing

Katrina Schlum<fnref n="1" />, Kurt Lamour<fnref n="1" />,<fnref n="2" />, Peter Tandy<fnref n="2" />, Scott J. Emrich<fnref n="1" />,<fnref n="3" />, Caroline Placidi de Bortoli<fnref n="2" />, Tejas Rao<fnref n="2" />, Diego M. Viteri Dillon<fnref n="4" />, Angela M. Linares-Ramirez<fnref n="5" /> and Juan Luis Jurat-Fuentes<fnref n="1" />,<fnref n="2" />,*

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

check for
updates

**Citation:** Schlum, K.; Lamour, K.; Tandy, P.; Emrich, S.J.; de Bortoli, C.P.; Rao, T.; Viteri Dillon, D.M.; Linares-Ramirez, A.M.; Jurat-Fuentes, J.L. Genetic Screening to Identify Candidate Resistance Alleles to Cry1F Corn in Fall Armyworm Using Targeted Sequencing. *Insects* **2021**, *12*, 618. https://doi.org/10.3390/insects12070618

Academic Editors: Ralf Nauen and
Gaelle Le Goff

Received: 26 May 2021

Accepted: 30 June 2021

Published: 8 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2021 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Simple Summary:** Monitoring of resistance alleles is critical to the sustainability of transgenic crops producing insecticidal Cry proteins. Highly sensitive and cost-effective DNA-based methods are needed to improve current bioassay-based resistance screening. Our goal was to evaluate the use of targeted sequencing in detecting known and novel candidate resistance alleles to Cry proteins. As a model, we used field-collected fall armyworm (*Spodoptera frugiperda*) from Puerto Rico, the first location reporting continued practical field-evolved *S. frugiperda* resistance to corn producing the Cry1F insecticidal protein, and sequenced the *SfABCC2* gene previously identified as critical to Cry1F toxicity. Targeted sequencing of *SfABCC2* detected a previously reported Cry1F resistance allele and mutations originally identified in populations from Brazil. Importantly, targeted sequencing also identified nonsynonymous and frameshift mutations as novel candidate resistance alleles. These results advocate for the use of targeted sequencing in screening for resistance alleles to Cry proteins and support potential gene flow, including resistance alleles, between *S. frugiperda* from Brazil and the Caribbean.

**Abstract:** Evolution of practical resistance is the main threat to the sustainability of transgenic crops producing insecticidal proteins from *Bacillus thuringiensis* (Bt crops). Monitoring of resistance to Cry and Vip3A proteins produced by Bt crops is critical to mitigate the development of resistance. Currently, Cry/Vip3A resistance allele monitoring is based on bioassays with larvae from inbreeding field-collected moths. As an alternative, DNA-based monitoring tools should increase sensitivity and reduce overall costs compared to bioassay-based screening methods. Here, we evaluated targeted sequencing as a method allowing detection of known and novel candidate resistance alleles to Cry proteins. As a model, we sequenced a Cry1F receptor gene (*SfABCC2*) in fall armyworm (*Spodoptera frugiperda*) moths from Puerto Rico, a location reporting continued practical field resistance to Cry1F-producing corn. Targeted sequencing detected a previously reported Cry1F resistance allele (*SfABCC2mut*), in addition to a resistance allele originally described in *S. frugiperda* populations from Brazil. Moreover, targeted sequencing detected mutations in *SfABCC2* as novel candidate resistance alleles. These results support further development of targeted sequencing for monitoring resistance to Bt crops and provide unexpected evidence for common resistance alleles in *S. frugiperda* from Brazil and Puerto Rico.

**Keywords:** *Spodoptera frugiperda*; fall armyworm; resistance; Cry1F; genotyping; targeted sequencing; resistance screen; ABCC2

<footnote n="1">¹ Genome Science and Technology Graduate Program, University of Tennessee, Knoxville, TN 37996, USA; kschlum@vols.utk.edu (K.S.); klamour@utk.edu (K.L.); semrich@utk.edu (S.J.E.)</footnote>

<footnote n="2">² Department of Entomology and Plant Pathology, University of Tennessee, Knoxville, TN 37996, USA; ptandy@vols.utk.edu (P.T.); cplacidi@utk.edu (C.P.d.B.); trao2@utk.edu (T.R.)</footnote>

<footnote n="3">³ Department of Electrical Engineering and Computer Science, University of Tennessee, Knoxville, TN 37996, USA</footnote>

<footnote n="4">⁴ Isabela Research Substation, Department of Agro-Environmental Sciences, University of Puerto Rico, Isabela, PR 00662, USA; diego.viteri@upr.edu</footnote>

<footnote n="5">⁵ Lajas Research Substation, Department of Agro-Environmental Sciences, University of Puerto Rico, Lajas, PR 00667, USA; angela.linares@upr.edu</footnote>

<footnote>* Correspondence: jurat@utk.edu; Tel.: +1-(865)-974-5931</footnote>
# 1. Introduction

The fall armyworm (*Spodoptera frugiperda*, J.E. Smith) (Lepidoptera: Noctuidae) is a global invasive pest affecting numerous food and fiber staple crops, although the highest damage is observed in corn (*Zea mays* L.) [1]. In North America, *S. frugiperda* moths display long-distance northward migration over several generations from overwintering sites in southern Texas and Florida, reaching Canada in the summer months [2]. Substantial gene flow also occurs between *S. frugiperda* populations from the Caribbean and Florida, but not between these populations and moths from South America [3].

Effective control of *S. frugiperda* in the Western hemisphere in addition to pesticides, has been provided by genetically modified corn and cotton producing insecticidal proteins from the bacterium *Bacillus thuringiensis* (*Bt*), namely Cry1F, Cry1Ab, Cry1A.105, Cry2Ab and Vip3Aa. In the past decade, resistance to Cry1F, Cry1Ab and Cry1A.105 was reported for *S. frugiperda* populations in Puerto Rico and the continental US (Florida and North Carolina) [4,5], Brazil [6] and Argentina [7]. Fall armyworm resistance to these Cry proteins is recessive and genetically linked to mutations in an ABC transporter superfamily C2 gene (*SfABCC2*) [8–10]. In Puerto Rico, two *SfABCC2* resistance alleles have been reported, a nonsense 2 bp insertion (*SfABCC2mut* allele) [8] and an insertion near the start of the fourth exon associated with aberrant splicing [9]. In contrast, resistance to Cry1F in Brazil is linked to missense mutations primarily localized to the fourth extracellular loop in SfABCC2, in support of this region containing a Cry1F binding domain [11] and being a hotspot for resistance mutations [10]. Genetic knockout of SfABCC2 leads to resistance against Cry1F [12,13] further supporting that truncation of SfABCC2 is predictive of a Cry1F-resistant phenotype.

Identification of resistance alleles in *S. frugiperda* is critical to developing highly sensitive DNA-based genotyping tools, which are less laborious and costly compared to F1 [14] or F2 [5] resistance screens. Successful *SfABCC2* genotyping efforts were described using TaqMan probes [8,9] and pyrosequencing [10] in populations from Puerto Rico, USA and Brazil. Results from these tests support high-resistance allele frequencies in areas where practical resistance had been reported, further associating identified *SfABCC2* alleles with field-relevant resistance. An unexpected observation, considering predicted migratory models [3], was not detecting the *SfABCC2mut* allele reported in Puerto Rico at migratory Caribbean destinations or in the Florida populations tested [8,9]. Relevant fitness costs are not associated with *SfABCC2mut* under laboratory conditions [15] and thus probably do not explain the limited distribution to Puerto Rico of *SfABCC2mut*.

Hi-Plex (HP) is a highly multiplexed PCR-based approach for targeted massive parallel sequencing with demonstrated efficacy in disease gene screening [16,17]. It involves amplifying small (e.g., 180 bp) overlapping portions of a larger target region (tiling) in one (or a few) multiplex PCR reaction(s) adding a unique sample-specific barcode sequence during the PCR reaction. The HP technology is attractive because it requires very little genomic DNA to amplify multiple target amplicons (e.g., 20–40 ng of template DNA to amplify a few hundred targets), the DNA can be highly complex (e.g., mixed with genomic DNA from many different organisms) and the overall cost in time and reagents is greatly reduced compared to F1 and F2 screens. The barcoded amplicons are pooled and sequenced on a next-generation sequencing (NGS) device, and the sample-specific amplicons mapped to a reference genome.

The goal of this project was to adapt and provide preliminary data for a new HP targeted screen for resistance alleles to *Bt* insecticidal proteins. We focused our analysis on *S. frugiperda* samples from Puerto Rico, as both *SfABCC2mut* and Cry1F resistance are highly frequent at this location [8,18]. Furthermore, mutations resulting in a truncated SfABCC2 protein result in a Cry1F-resistant phenotype [8,9,12,13], potentially allowing for phenotypic predictions. Analysis of Taqman and targeted SfABCC2 sequences detected the known SfABCC2mut allele in 34–37% of tested samples, depending on the technique used. In addition, the HP approach coupled to a bioinformatics pipeline detected one nonsense and nine frameshift SfABCC2 mutations as novel candidate Cry1F resistance