Whole-Genome Resequencing of the VGSC Reveals the Evolutionary Mechanism of Pesticide Resistance in Liriomyza trifolii in Hainan

1. Introduction

Liriomyza trifolii (Burgess), an invasive pest species in China, initiated its global dissemination from Florida, USA, in 1968 [1]. Currently, the distribution of the L. trifolii is primarily concentrated in the southeastern coastal regions of China, encompassing Hainan, Guangdong, Guangxi, Shanghai, Fujian, Zhejiang, Jiangsu, and Taiwan. L. trifolii exhibits a broad host range, impacting crops such as legumes, Solanaceae, Cucurbitaceae, and cruciferous plants [2], with documented evidence of host range expansion [3]. This pest is distinguished by its diminutive size, high reproductive capacity, brief development cycle, and rapid spread dispersal, frequently resulting in overlapping generations, which predisposes it to outbreaks and substantial agricultural damage [4]. In warmer regions, L. trifolii can cause continuous damage throughout the year due to overlapping generations.

Hainan Province, located in southern China, functions as the principal distribution hub for winter fruits and vegetables in the region, with an annual cultivation area exceeding 60,000 hectares. That’s a significant challenge for agriculture in Hainan. The year-round cowpea cultivation creates a perfect storm for L. trifolii populations, essentially providing a continuous host and breeding ground. This, combined with the extended and frequent pesticide exposure, creates a classic scenario for the rapid development of pesticide resistance. Currently, both domestically and internationally, chemical pesticides remain the primary method for controlling the L. trifolii. The widespread and excessive application of pesticides has led to a continuous rise in field population resistance.

Pyrethroids are a significant class of neurotoxic pesticides [5]. Their mode of action involves disrupting the function of the insect’s peripheral and central nervous systems [6]. Specifically, pyrethroids target voltage-gated sodium channels (VGSC) as their main site of activity [7]. Mutations at VGSC sites can lead to knockdown resistance (kdr), the primary mechanism of pyrethroid resistance, by reducing the pesticide’s affinity for its target [8,9]. This same VGSC target is also exploited by other insecticides, including indoxacarb and flufenidin [10,11,12,13]. VGSCs are comprised of four homologous internal domains (I–IV), each containing six transmembrane α-helices (S1–S6) (Figure S1). These domains are interconnected by intracellular linkers, while the helices are linked by either intracellular or extracellular loops [14]. The homologous channel model suggests that pyrethroids and VGSCs possess two distinct binding regions: PyR1 (II L45-II S5-III S6) [15] and PyR2 (I L45-I S5-I S6-IIS6) [16]. Moreover, studies on mammals have demonstrated that the drug Cenobamate binds to the canonical human Na_v1.5 ion channels by interacting primarily with S6 in domains II and III, along with a sole residue (Ile 1768) in domain IV, which differs slightly from insects [17]. Research has extensively focused on mutations within these binding regions, as such mutations are likely to diminish the sensitivity of insect VGSCs to pyrethroids, thereby contributing to increased insect resistance to these pyrethroids.

Whole-genome resequencing (WGS) is a process that swiftly and accurately identifies all the “differences” between a specific individual and the “standard map” (reference genome) based on the existing “map” information. WGS has emerged as a prominent application in Next-Generation Sequencing (NGS). Relevant studies have conducted genotyping on populations exhibiting varying levels of pyrethroids resistance in northern Brazil, focusing on genetic polymorphisms within the genomes of Aedes aegypti with distinct characteristics. The allele frequencies across the entire genome were assessed using SNP chips, leading to the identification of single nucleotide polymorphisms (SNPs) directly associated with resistance, as well as one superior SNP pair [18]. WGS of Spodoptera frugiperda from the Indochinese Peninsula to northern China (2019–2023) identified eight stable genomic markers (six SNPs, two InDels) that clearly distinguish eastern and western populations. This work established a reliable molecular toolkit for tracing the origin of individual moths [19].

This study employed WGS to investigate the underlying mechanisms of high resistance in L. trifolii populations from Hainan. By comparing the VGSC across various geographic strains and a susceptible reference, four critical VGSC mutations were identified. These mutations confer high levels of resistance and cross-resistance to pyrethroids and indoxacarb, offering crucial molecular targets for future resistance management strategies.

2. Results

2.1. Complete Genomic Information of VGSC in L. trifolii

Following the analysis of first-generation sequencing results, the complete genomic sequence of the VGSC of L. trifolii was determined to be 36,774 base pairs (bp). This sequence comprises specific lengths for Domain I (6,696 bp), Domain II (10,162 bp), Domain III (11,025 bp), and Domain IV (8,891 bp). Within the entire genomic sequence, the intron and coding sequence (CDS) lengths are 30,360 bp and 6,414 bp, respectively (Figure S2). The focused 6,414 bp CDS segment contains a total of 2,138 amino acids, including the stop codon. This segment encodes 415, 614, 531, and 578 amino acids for Domain I, Domain II, Domain III, and Domain IV, respectively (Figure 1).

2.2. Information of Exonic Synonymous sSNP in L. trifolii

Based on codon degeneracy, the 31 identified sSNPs were classified as L145L, I290I, S319S, Y354Y, A470A, L474L, E496E, E505E, E661E, A685A, L733L, L778L, T814T, Y866Y, N931N, S978S, A1112A, L1231L, D1232D, P1282P, F1372F, I1471I, Y1483Y, H1996H, Q2006Q, R2011R, G2022G, G2033G, A2038A, G2040G, and A2050A. These sSNPs are primarily distributed across the following regions: “D I -S1”, “D I -L-S56”, “L-D I - II”, “D II -S1”, “D II -L-S23”, “D II -S5”, “D II -L-S56”, “L-D II - III”, “D III -S3”, “D III -L-S56” and “L-D IV <“. The positions of these sSNPs within the VGSC, along with the nucleotide alterations and occurrence frequencies, are detailed in Table 1. Notably, 14 sSNPs were found in all five field populations (MY, TS, DA, JY, TY), specifically L145L, I290I, S319S, E496E, E505E, L733L, N931N, S978S, A1112A, P1282P, Q2006Q, R2011R, G2033G, and A2050A. Among these, Q2006Q exhibited the highest occurrence rate at 94.31%, while N931N had the lowest at 7.08%. Furthermore, both the highest and lowest frequencies were observed in the JY population. Of the 31 sSNPs, A470A, L1231L, D1232D, I1471I, Y1483Y, A2038A, and G2040G were found exclusively in a single population, with frequencies of 8.57% (TS), 11.76% (DA), 12.24% (DA), 10.10% (TS), 16.19% (TS), 7.04% (TS), and 8.33% (TS), respectively (Figure 2).

2.3. Information of Exonic nsSNP in L. trifolii

After conducting multiple comparisons and analyses of the CDS region in the genome resequencing data, our results identified 6 nsSNPs in the VGSC of L. trifolii: M922T (918), T933I, L1018F (1014), Q1285H, V1845I (1851), and D2036E. These nsSNPs are primarily situated in specific regions of the VGSC: D II -L-S45 (M922T), D II -S5 (T933I), D II -S6 (L1018F), L-D II - III (Q1285H), D IV -S6 (V1845I), and L-D IV < (D2036E) (Figure 3A). The positions, nucleotide alterations, and frequencies of these nsSNPs in the VGSC are detailed in Table 2. Notably, 3 nsSNPs, namely M922T, L1018F, and V1845I, were present across all five field populations (MY, TS, DA, JY, TY). Among these, L1018F exhibited the highest occurrence rate of 100.00%, which was consistent across all five populations. M922T exhibits the highest incidence rate in JY (90.91%) and the lowest in MY (58.46%). Similarly, V1845I shows the highest frequency in JY (96.52%) and the lowest in MY (36.54%), mirroring the pattern observed for M922T. T933I and D2036E are unique to the TS population among the 6 nsSNPs. The frequencies of Q1285H and V1845I in the DA and JY populations are identical at 80.00% and 96.52%, respectively. Figure 3B illustrates the frequency comparison of various nsSNPs across populations in a north-to-south order based on decreasing latitude.

3. Discussion

The high level of resistance and unique reproductive capabilities are the primary factors contributing to the L. trifolii ‘s status as a significant pest [20]. For decades, deltamethrin has been extensively employed to manage the L. trifolii. However, subsequent studies have demonstrated that deltamethrin is no longer effective against this pest [21]. In 1979, permethrin was introduced in California, United States, for the control of the L. trifolii. Unfortunately, its effectiveness diminished due to the pest’s increased resistance [22,23]. Additionally, over the years, the L. trifolii exhibited varying degrees of resistance to permethrin, deltamethrin, and several other insecticides [24,25]. The sustainable management of crop pests in agricultural fields has been significantly compromised by the emergence of pesticide resistance [26], a challenge that is particularly pronounced in tropical regions [27].

For the identified nsSNPs, M918T (a nsSNP of M922T in L. trifolii), conferred exceptionally high resistance to deltamethrin in rat VGSC [28]. Furthermore, in 2011 and 2012, M918T was confirmed to be linked to resistance against pyrethroid insecticides in Tetranychus evansi and Aphis gossypii [29]. The M918L variant at the same position was identified in Hyelella azteca [30]. Additionally, M918V was discovered in Bemisia tabaci [31]. When examined in Xenopus oocytes, various substitutions can confer differing levels of resistance. Distinct mutant VGSCs offer varying degrees of protection against type I and type II pyrethroids [16]. This phenomenon has been observed at the M918 (T/ L/ V) site, where the M918T substitution conferred the highest level of protection against permethrin and deltamethrin [28]. In this study, only the M918T variant was identified in the L. trifolii, and the mutation frequency exhibited a gradual increase as latitude decreased, ultimately reaching a maximum of 90.91%. This finding suggests that L. trifolii in Hainan Province has developed an exceptionally high resistance to pyrethroids. This explains the widespread emergence of the potent M918T substitution across different populations and its remarkably high frequency. Furthermore, the frequency of these mutations correlates with rising temperatures at lower latitudes. This trend is likely attributed to the shortened generation cycle of L. trifolii at higher temperatures, which accelerating the accumulation of insecticide resistance in their offspring.

Among the six nsSNPs identified in this study, mutations consistently linked to pyrethroids resistance include M918, as well as L1014 (F, C, H, S, or W) (L1018F in L. trifolii) and T933 (I, C, or V). Notably, L1014F is situated in PyR1, similar to M918, whereas T933 is found in PyR2. Therefore, M918 and L1014F within PyR1 will be analyzed in our results firstly, followed by a discussion of T933 located in PyR2.

The first mutation identified and confirmed as kdr-related was L1014F (L1018F in L. trifolii) [32]. L1014F provide variable levels of protection to Type I or Type II pyrethroids or DDT [33]. Prior research has demonstrated that the L1014F single mutation confers resistance to pyrethroids in various species, including Anopheles stephensi [34], A. gossypii [29], Blattella germanica [8,35], Ctenocephalides felis [36], Cydia pomonella [37], Frankliniella occidentalis [38], Leptinotarsa decemlineata [39], Liriomyza huidobrensis [40], Meligethes aeneus [41], Musca domestica [35,42,43], Sitobion avenae [44], Triatoma infestans [45] and Liriomyza sativae [40], etc. The variants L1014S, L1014H, L1014C, and L1014W were identified at position 1014 [40,46,47]. In contrast, only one mutant type, L1014F, was observed in the L. trifolii, exhibiting a frequency of occurrence of 100% across all five populations. This finding indicates that the resistance of the L. trifolii in Hainan Province to pyrethroids has reached an exceptionally high level.

The T933I, corresponding to the amino acid number in the L. trifolii, has been linked to resistance against pyrethroids in various pests, including Thrips palmi [48], Trialeurodes vaporariorum [49], L. decemlineata [50] and Thrips tabaci [51]. Other single mutations at this position include T933C and T933V. Additionally, T929C was identified in F. occidentalis [38], while T933V has been observed in C. felis [36] and F. occidentalis [38]. In this study, only the T933I mutation was identified, and this mutation only occurred at a relatively low frequency of 8.57% in the TS. The outcome could be attributed to TS being a classic “mountain city” enclosed by mountains from all directions, leading to reproductive isolation or minimal gene flow with other populations. In small and secluded populations, particular mutations are prone to enhancing survival chances, despite lacking an immediate advantage. The T933I evidently confers a survival advantage special to the environmental challenges encountered by this population, including the intensity and frequency of local pesticide use and the cold temperatures prevalent in mountainous regions.

Prior research has demonstrated that L. trifolii populations with avermectin tolerance display the highest resistance to cypermethrin, a type II pyrethroid [52]. Consequently, cross-resistance between pyrethroids and other insecticides has surfaced in L. trifolii. Moreover, in the Sanya City field of Hainan Province, China, in 2022, the sensitivity of second-instar L. trifolii larvae to indoxacarb decreased by 776.17 times, with resistance to avermectin exhibiting a consistent upward trajectory [53]. One of the nsSNP identified in this study, V1845I (amino acid position in the L. trifolii), has been confirmed to be associated with resistance to indoxacarb in Plutella xylostella [12]. Subsequently, V1845I resistance mutation were also detected in Tuta absoluta [11]. Furthermore, this mutation has been reported in the L. trifolii [27]. Research has confirmed, using a two-electrode voltage clamp system, that the V1845I mutant channel are more resistant to indoxacarb, DCJW and metaflumizone than wild type channels [10]. Moreover, the resistance mutation V1845I of the VGSC blocker insecticides (SCBIs) was functionally validated in Drosophila through molecular simulation and genomic engineering [54]. These findings suggest that the L. trifolii populations in Hainan Province, China, not only demonstrates exceptionally high resistance to pyrethroids but also shows a relatively elevated level of cross-resistance to indoxacarb. The findings of this study indicate that as latitude decreases, the mutation frequency of V1845I exhibits a general upward trend, peaking at 96.52%. The mutation frequency in the TS region is notably higher than that observed in high-latitude areas. This discrepancy may be attributed to the widespread application of indoxacarb in the agricultural practices of the TS region.

Moreover, research has demonstrated that double mutations, in contrast to single mutations, nearly abolish the sensitivity of Drosophila melanogaster VGSC to deltamethrin [28,39]. This finding suggests that the impact of the double mutation surpasses that of the single mutation. Notably, the mutation Q1285H+V1845I exhibited identical mutation frequencies in the DA and JY populations, at 80.00% and 96.52%, respectively. Due to the sequencing methods and characteristics of WGS, it remains uncertain whether Q1285H+V1845I represents a genuine tandem double mutation. Nevertheless, based on frequency analysis, there is a likelihood of the occurrence of the double mutant Q1285H+V1845I. The functional implications of this dual process necessitate validation through electrophysiology-related experiments in future investigations. The D2036E variant identified in this study is located in the non-domain region of the VGSC of the L. trifolii. Further investigation is necessary to determine its potential association with drug resistance in this species.

In addition to the numerous nsSNPs analyzed above, our results identified a substantial number of sSNPs in the VGSC of the L. trifolii, with their prevalence significantly exceeding that of nsSNPs. Historically referred to as “silent mutations”, sSNPs are now understood to be not entirely silent. In general, sSNPs are still considered benign or minimally harmful. Furthermore, accurately quantifying the number of harmful variations (specifically nsSNPs), and neutral variations (specifically sSNPs), present in a population can elucidate the influence of the population’s developmental history on genetic variation. This assessment can also determine whether a small-scale population is prone to the accumulation of harmful variations, potentially resulting in uncontrolled mutations or even species extinction and outbreaks [55,56]. In evolutionary research, the ratio of nsSNP to sSNP (dN/dS) serves as a crucial indicator for assessing whether a gene is subject to natural selection pressure. When natural selection facilitates alterations in protein sequences, the dN/dS ratio is anticipated to exceed 1. Conversely, when natural selection restricts changes in proteins, the ratio falls below 1 [58,59]. This interpretation of dN/dS is further substantiated by theoretical analyses of the relationship between dN/dS statistics and potential selection pressure within the Wright-Fisher model [60]. In this study, all five populations adhered to the principle of natural selection, as indicated by a dN/dS ratio of less than 1. Our results propose that this observation signifies several important aspects. Firstly, in the VGSC of the L. trifolii, the majority of deleterious mutations (nsSNPs) are purged from the population. Individuals harboring these mutations, which disrupt the original structure and function of proteins, experience significantly reduced survival and reproductive success. Consequently, natural selection persistently eliminates these harmful nsSNPs, resulting in a dN/dS ratio of less than 1. Next, the functionality of VGSC is critically important. Their protein sequences are essential for survival and reproduction, such that any non-synonymous alteration in their amino acid sequences carries substantial implications. Furthermore, the emergence of any nsSNP likely represents a significant site that aids the population in resisting natural selection pressures, such as pesticide resistance. Finally, a dN/dS ratio of less than 1 also suggests that the VGSC sequence has remained highly conserved throughout a prolonged evolutionary period, with minimal changes in its amino acid sequence across different species or generations. Key regions, including structural or functional sites of VGSCs, are subject to considerable selection pressure to preserve their existing functions, thereby limiting the potential for alteration. In Hainan Province, characterized by persistent high temperatures, humidity, and continuous pesticide use, further investigation is warranted regarding the unique evolution of VGSC in L. trifolii and associated nsSNPs.

4. Materials and Methods

4.1. Insects Collection

The sensitive strain of L. trifolii (LT-S) used in this study was provided by Professor Yuzhou Du of Yangzhou University. This strain has been maintained in the laboratory without exposure to insecticides for over two decades. Wild populations selected for WGS were gathered from five locations: Tian-Ya (TY: N 18.31, E 109.48) and Ji-Yang (JY: N 18.29, E 109.54) in Sanya, Mao-Yang (MY: N 18.94, E 109.51) and Tong-Shi (TS: N 18.79, E 109.52) in Wuzhishan, and Da-An in Ledong (DA: N 18.47, E 108.90) (Figure 4). The entire life cycle of the insects was raised at a temperature of 26 ± 1 ℃, the humidity should be maintained at 70% ± 5% and L: D = 16h: 8 h. The adult insects laid eggs and the larvae were raised using cultivated bean seedlings (Baijiali cowpea seeds, USA) that were not exposed to any pesticides.

4.2. Acquisition of Complete Genomic of VGSC in L. trifolii and Sequencing WGS

Genomic DNA extracted from LT-S was single-headed and processed using the Universal Genomic DNA Extraction Kit (Beijing Solarbio Science & Technology Co., Ltd., China). Segmented gene amplification on the complete genomic sequence of VGSC of L. trifolii was conducted utilizing the Phanta-Max Super-Fidelity DNA Polymerase kit (Vazyme Biotech Co., Ltd., China). By integrating the genomic data of the VGSC from the L. trifolii (GenBank: GCA_001014935.1), primers were formulated for cloning to acquire the complete genome sequence of the susceptible strain of the VGSC from the L. trifolii. Details of the primers and annealing temperatures employed for fragment amplification can be found in Table 3. Subsequently, PCR products were submitted to Tsingke Biotechnology Co., Ltd. (China) for first-generation sequencing, and the DNASTAR Lasergene 11 Core Suite software was utilized to assemble the complete genomic sequence of VGSC.

Large quantities of sensitive strains cultivated indoors and resistant populations gathered from the field were sent to the Beijing Genomics Institute (BGI Genomics Co., Ltd.) for WGS. Genomic data analysis was conducted utilizing the GeneAn cloud platform.

4.3. Analysis and Acquisition of sSNP and nsSNP in VGSC of the L. trifolii

The genomic fragment sequences of the VGSC from the LT-S strain were concatenated to derive the complete full-length genomic sequence of the VGSC of L. trifolii. This complete VGSC genomic sequence will be annotated according to the Coding Sequence (CDS) (GenBank: MT648287.1), and the nucleotide position information of the intron and exon segments will be scrutinized. Mutants situated in the exon region were filtered from the outcomes of WGS using the aforementioned nucleotide position information. Subsequently, by correlating the codon information of each mutation in the exon with the CDS position, our results classified the mutations as synonymous single nucleotide polymorphism (sSNP) or non-synonymous single nucleotide polymorphism (nsSNP), determined their positions on the VGSC, and computed their frequencies of occurrence. The approach for distinguishing between sSNPs and nsSNPs is illustrated in Figure S3, where the section highlighted in red elucidates the method for identifying sSNPs as discussed in this study.