Introduction

Human-induced climate change, particularly the continued change in temperature and precipitation patterns (IPCC, 2023), is affecting the survival and distribution of species (Fitzpatrick and Edelsparre, 2018), and impacts biodiversity and ecosystem services (Lawlor et al., 2024). Terrestrial insects are particularly vulnerable to temperature change because they are poikilothermic animals and small with high surface area to volume ratio (Wang et al., 2022). Over the past century, many insect populations have shown declining abundance, with climate change considered a major contributing factor after habitat loss and use of agricultural pesticides (Cardoso et al., 2020). It is predicted that a species distributed across latitudinal or elevational gradients shows a positive relationship between the thermal tolerance range and the level of climatic variability (Harvey et al., 2020). Tropical insects have a narrower thermal range, indicating that even a small increase in temperature can lead to a collapse in their populations (Harvey et al., 2022). Global climate change may also alter the geographical distribution of insect pests, allowing those previously confined by temperature barriers to spread to new areas, posing a threat to natural and managed ecosystems, crop production and food security (Deutsch et al., 2008; Outhwaite et al., 2022). Species with higher genetic diversity and adaptive potential may better cope with varying or extreme temperature environments, giving them an advantage to move to new habitats.

Adaptive evolution is a crucial mechanism for insects to tackle climate change and maintain or expand their geographical ranges (McCulloch and Waters, 2022; Burc et al., 2025). Understanding the genetic basis of thermal adaptation is essential for predicting how species will respond to rapid global climate change (Gibert et al., 2019). Genetic mutations, particularly non-synonymous mutations, are the primary sources of DNA sequence variability and can alter protein structure or function (Belfield et al., 2018). Increased thermal tolerance in organisms can originate from such mutations. For example, nonsynonymous mutations in the AcVIAAT gene of the eastern honeybee, Apis cerana (e.g., P42L substitution), are associated with enhanced thermal adaptation (Li et al., 2024c). During thermal evolution, alpine ground beetle, Nebria vandykei, can achieve survival capability in extreme thermal environments through adaptive selection of mutations in key genes (TREH, EIF3A, LRPPRC, etc.) coupled with immediate responses from heat shock proteins (Schoville et al., 2024). Although significant progress has been made in identifying genetic mutations that facilitate adaptation to extreme environments in non-model organisms, our understanding remains incomplete, especially regarding the integrated molecular, metabolic, and life history changes that collectively enable thermal adaptation.

In the case of diamondback moth (DBM), Plutella xylostella (Lepidoptera: Plutellidae), a globally serious pest of cruciferous plants, the availability of complete genome sequences and a large dataset of genome-wide single nucleotide polymorphism (SNP) sequences (You et al., 2013; You et al., 2020), allow for the investigation of its genetic plasticity, evolutionary history, and climate adaptation to changes in temperatures. Our recent research has shown that P. xylostella originated in South America and has expanded to all countries and regions where cruciferous crops are grown (Furlong et al., 2013; You et al., 2020). Based on the genomic data of field populations from 114 locations in a diverse range of biogeographical zones, we have demonstrated that P. xylostella is genetically capable of tolerating anticipated future climates in most regions of the world and will remain a global pest beyond 2050 (Chen et al., 2021). Here, we investigate the mechanisms of P. xylostella’s genetic adaptation and evolutionary responses to different thermal environments. Specifically, we employ thermal regime patterns of 12h/12h in hot (32°C/27°C) or cold (15°C/10°C) environments over the course of three years, compared to the favorable constant condition at 26°C, to investigate the adaptive evolution of P. xylostella in climate-controlled chambers.

Age-stage, two sex life tables of P. xylostella measured the demographic variation of the three P. xylostella strains evolved in the favorable (ancestral), hot and cold environments. Using metabolomic and transcriptomic analyses, we identified the key genes that could facilitate the adaptation of P. xylostella to thermal extremes. Our results showed that a large number of differentially expressed genes and metabolites were produced in populations adapted to high and low temperatures through multi-generational selection. We find the mutant of a key gene, PxSODC, which can alter the superoxide dismutase activity and increase the ability to scavenge superoxide anions, thereby significantly affecting the adaptability under high and low temperature environments. CRISPR-Cas9 was used to functionally validate the role of PxSODC in facilitating adaptive evolution and influencing regulatory networks. These findings not only improve our understanding of how novel thermal environments can drive genetically-based adaptive evolution manifested through integrated transcriptomic, metabolic, and life history changes in P. xylostella but also help to predict the distribution and dynamics of its populations in the context of global climate change.

Results

Life history trait divergence among temperature-adapted strains

Compared to the ancestral strain, both hot and cold strains had a significantly shorter preadult duration than the ancestral strain. The hot strain also have significantly shorter female longevity, and oviposition days, while the cold strain had an longer male longevity when compared to the ancestral strain (Supplemental Table S1). The female fecundity, intrinsic rate of increase (r) and finite rate of increase (λ) were all significantly higher in the hot strain than the cold strain which was not significantly different than the ancestral strain (Supplemental Table S1). Overall, following three years of evolution under contrasting thermal regimes, the hot and cold strains exhibited divergent life history traits compared to the ancestral strain. The hot strain exhibited accelerated life cycle and increased fecundity while the cold strain extended male longevity.

Based on the age-stage, two-sex life tables of the strains, we generated curves of the age-stage specific survival rate (s_xj), age-specific survival rate (l_x), female age-specific fecundity (f_x), and population age-specific fecundity (m_x) (Supplemental Figure S1). These curves present the dynamic variation in survival and fecundity of different strains and allow accurate estimation of daily survival and fecundity. The results showed that the daily survival rates at the pre-adult stages of the hot and ancestral strains were higher than that of the cold strain, with the maximum daily survival rate of HS = 78.9%, AS = 76.7% and CS = 66.7%. The daily fecundity of the hot and ancestral strains was higher than that of the cold strain, with the peak daily egg counts of HS = 53.22, AS = 53.00 and CS = 45.10 (Supplemental Table S1).

To assess the evolved thermal tolerance of the temperature-adapted strains, we further examined the stage-specific survival rates of the hot and ancestral strains under extremely high temperatures, as well as the supercooling and freezing points of the cold and ancestral strains at pupae stage. The survival rates of eggs, 3^rd-instar larvae and adults in the hot strain were higher than those of the ancestral strain at 42°C, and the survival rate of pupae in the hot strain was higher than that in the ancestral strain at 43°C and 44°C (Figure 1A). The supercooling and freezing points of pupae in the cold strain were lower than those in the ancestral strain (Figure 1B). The variation in survival rates and the supercooling/freezing points at extreme temperatures suggest that the hot and cold strains of P. xylostella have undergone profound adaptive adjustments.

Figure 1.

Omics-based evidence for adaptive evolution

Our previous studies have identified metabolites such as trehalose and very long chain fatty acids that play a role in adaptation of P. xylostella to both high and low temperatures (Zhou et al., 2022; Lei et al., 2023). Here, we performed a broad analysis of targeted metabolites of 3^rd-instar larvae of each strain using high-throughput UPLC-MS/MS. A total of 781 metabolites were identified, including 199 amino acids and their metabolites, 146 lipids, 90 organic acids and their derivatives, 78 nucleotides and their metabolites, 61 heterocyclic compounds, 45 benzene and substituted derivatives, 42 alcohols and amines, 37 carboxylic acids and derivatives, 21 coenzymes and vitamins, and 62 other metabolites (Figure 2A). Principal component analysis (PCA) and inter-sample correlation heat maps revealed significant metabolic changes in the 3^rd-instar larvae from the ancestral strain to the hot and cold strains (Figure 2B; Supplemental Figure S2A). These comprised 77 differential metabolites (34 up-regulated, 43 down-regulated compared to the ancestral strain) in the hot strain, and 37 differential metabolites (13 up-regulated, 24 down-regulated compared to the ancestral strain) in the cold strain (Figure 2C). These differential metabolic profiles are likely to play crucial roles in the adaptation of P. xylostella to high and low temperature environments.

Figure 2.

Compared to the ancestral strain, the common differential metabolites of the hot and cold strains included lipids, amino acids and their metabolites, organic acids and their derivatives, nucleotides and their metabolites, and benzene and substituted derivatives (Figure 2D; Supplemental Figure S2B). Inter-replicate analysis of differential metabolites showed a low correlation between the ancestral strain and hot/cold strains, but a high correlation between the hot and cold strains. Notably, 30 common metabolites were identified across the differential sets based on comparison of the hot/cold strains to the ancestral strain (Figure 2E; Supplemental Figure S2C). These metabolites, except for N1, N8-diacetylpiperidine, exhibited similar fold changes when comparing the hot/cold strains to the ancestral strain (Figure 2F). These results indicate that P. xylostella responds to different environmental stresses by regulating similar metabolic pathways. Further analysis revealed a reduction in most of the lipid metabolites by comparing the hot/cold strains to the ancestral strain (Figure 2F), providing evidence for adaptive evolution of P. xylostella in response to high and low temperatures, with the synthesis, degradation and metabolism of lipid playing an important role.

We then profiled and compared the transcriptomes of the three strains to identify the key genes involved in adaptation of P. xylostella to temperature extremes. This revealed significant variation in gene expression among strains (Figure 3A; Supplemental Figure S3A), with 1364 (825 up-regulated, 539 down-regulated) and 2029 (1205 up-regulated, 824 down-regulated) genes differentially expressed in the hot and cold strains, compared to the ancestral strain, respectively (Figure 3B). Pearson correlation showed, in contrast to the metabolomics data, a lack of strong correlation between the differentially expressed genes of the hot/cold strains and the ancestral strain (Figure 3C), with 498 common differentially expressed genes (Figure 3D). However, KEGG analysis revealed that the differentially expressed genes between the hot/cold strains and the ancestral strain were enriched in similar pathways, such as transport and catabolism, and the digestive system (Figure 3E), indicating that while multiple genes are involved in the adaptation of P. xylostella to high and low temperatures, a relatively limited range of biological functions might be affected.

Figure 3.

The gene expression-based clustering tree using a weighted gene co-expression network analysis (WGCNA) was divided into 29 modules as shown with arbitrary color codes (Figure 4A; Supplemental Figure S3B). The turquoise module contained the most genes (3463), while the grey module contained the fewest (31) (Figure 4B). Selecting the common differential metabolites (30 in total) shared between the hot and cold strains as compared to the ancestral strain, we performed a correlation analysis with co-expressed networks and discovered the strongest correlations between genes in the red module and 28 metabolites (Figure 4C; Supplemental Figure S3C). Further analysis revealed that 79 genes within the red module were differentially expressed in the hot and cold strains when compared with the ancestral strain (Figure 4D). These results suggest that genes in the red module are most likely to play a critical regulatory role in the adaptation of P. xylostella to extreme temperatures.

Figure 4.

Genetic basis of temperature adaptation

To further elucidate the genetic basis of P. xylostella adaptation to hot and cold environments, we analyzed 15 annotated genes with high expression levels (FPKM > 10) within the red module (Figure 4D) and identified 11 genes being successfully amplified. Comparative results revealed eight genes with nonsynonymous mutations and one with a synonymous mutation in both the hot and cold strains (Supplemental Table S2). Among these genes, we focused on the role of Px04C00666 (PxSODC) in the temperature adaptation of P. xylostella because the deletion of SOD genes can alter the response of insects to abiotic stresses including temperature (Bittner et al., 2019; Quan et al., 2024).

The NCBI database predicted that PxSODC contained three exons and two introns, with a conserved domain belonging to the copper-zinc superoxide dismutase superfamily (Supplemental Figure S4A) which plays an antioxidative role in cellular defense systems, protecting cells from damage caused by reactive oxygen species (Fridovich, 1995). The open reading frame of PxSODC was 633 bp and encodes 210 amino acids. Expasy predicted that the molecular weight of the PxSODC protein was 22,168.24 Da, with the isoelectric point being 6.29. Its secondary structure mainly consisted of strands, helices and coils (Supplemental Figure S4B). Evolutionary analysis showed PxSODC of P. xylostella clustered with that of other Lepidoptera insects such as Brenthis ino, Vanessa tameamea, and Vanessa cardui, indicating its conserved evolution within this taxonomic order, but distinct from Diptera flies and Hemiptera aphids (Supplemental Figure S4C). The coding region of PxSODC in the hot and cold strains had 23 SNP sites, including 20 synonymous and three non-synonymous mutations (Leu9-Val9, Lys25-Gln25, Leu194-Met194) (Figure 5A, 5B). Leu9-Val9 and Leu194-Met194 mutations were involved in the substitution of hydrophobic amino acids, with the Leu194-Met194 mutation exhibiting a higher frequency in both the hot and cold strains. The expression of PxSODC at different developmental stages of the hot and cold strains was significantly lower than that of the ancestral strain (Supplemental Figure S5A). After 2 h exposure of the 3^rd-instar larvae to the stress of high (32°, 34°, 36°, 38° and 40°C) or low (12°, 10°, 8°, 6° and 4°C) temperature environments, the expression of PxSODC in the hot and cold strains was significantly lower than in the ancestral strain (Supplemental Figure S5B, S5C). This suggests that the PxSODC gene can alter its expression pattern and function in response to environmental change and contribute to the temperature-adaptive evolution in P. xylostella.

Figure 5.

To explore the specific function of PxSODC in temperature adaptation of P. xylostella, we successfully knocked out PxSODC in the ancestral strain using CRISPR/Cas9 technology. Of 162 eggs treated with CRISPR/Cas9, 75 successfully developed into adults. We confirmed three mutant strains in the G0 generation of P. xylostella: +1 bp (SODC-MU1), +2 bp (SODC-MU2) and −1 bp (SODC-MU3). Self-crossing continued, and three types of homozygous mutation were obtained in the G5 generation (Supplemental Figure S6). Age-stage, two sex life tables were constructed for the ancestral strain, and the three SODC-MU mutant (MU1, MU2, MU3) strains under constant favorable (26°C), hot (32°C/27°C) and cold (15°C/10°C) environments to examine their phenotypic traits including development, survival, fecundity and fitness (Supplemental Figure S7; Supplemental Tables S3-S5). Overall, the three mutant strains had a prolonged development time (T), lower survival rates (s_xj and l_x), and reduction in the fecundity (f_x and m_x) and population fitness parameters (r, λ, R₀) compared to the ancestral strain. The observed change in these parameters indicated that the loss of PxSODC gene function significantly affected the life history and population fitness of P. xylostella, particularly in the hot/cold environments. Prolonged developmental time and reduced fecundity suggest that the physiological functions of the mutant strains may be damaged when exposed to both high and low temperature conditions.

We then carried out a detailed study on the response of ancestral and mutant strains (AS, SODC-MU1, SODC-MU2 and SODC-MU3) to extremely high temperature conditions. At 42°C, the survival rates of eggs, 3^rd-instar larvae, female adults and male adults of the mutant strains were significantly lower from those of the ancestral strain at several time points (Figure 5C). However, survival rates of the mutant pupae exposed to the high temperature (43°C) were not significantly different from those of the ancestral strain at different time points (Figure 5C). The pupal stage exhibited enhanced stability to environmental temperature fluctuations compared to other life stages, possibly due to a special heat resistance mechanism in the pupae, such as protective chrysalis and U-shaped metabolism (Kaiser et al., 2010; Barros-Cordeiro et al., 2014). In addition, supercooling and freezing points of the mutant strains (MU1: −21.32±0.41 and −9.75±0.38; MU2: −21.50±0.38 and −9.93±0.43; MU3: −21.23±0.48 and −9.94±0.41) were significantly higher than those of the ancestral strain (−23.09±0.26 and −11.58±0.52) at pupal stage (Figure 5D), indicating a key role of the PxSODC gene in the adaptability and tolerance of P. xylostella to extreme temperatures.

To further investigate the effect of PxSODC gene mutations on the temperature adaptability of P. xylostella, we identified five genes from the same SOD family in transcriptomes of the 3^rd-instar larvae from the three tested strains. We found that Px04C00505 and Px13C00423 showed SNP mutations in the hot and cold strains, whereas Px20C00248, Px15C00224 and Px15C00223 were not mutated (Supplemental Table S6). Further comparison of gene expression across different strains revealed that, relative to the ancestral strain, the expression levels of PxSODC, Px04C00505, and Px13C00423 were significantly reduced in the hot and cold strains, while the remaining genes maintained stable or increased expression levels (Figure 6A). Concurrently, SOD activity decreased in the hot and cold strains, along with a reduction in O₂^- levels (Figure 6B, 6C). When SOD expression and activity, as well as O₂^- levels, were compared under different temperature conditions between the ancestral and hot/cold strains, similar patterns were observed (Figure 6D, 6E, 6F). These results demonstrate that the hot and cold strains could increase their O₂^-scavenging capacity through non-synonymous mutations, thereby better adapting to high or low temperature environments. Compared to the ancestral strain, expression of the mutated SODC genes (Px04C00505 and Px13C00423) was increased in the male adults in SODC-MU1 and SODC-MU2. The expression levels of the non-mutated genes Px20C00248, Px15C00224 and Px15C00223 also increased (Figure 6G). Further, SOD activity decreased, while O₂^- levels increased in the two mutant strains (Figure 6H, 6I), which were unable to fully compensate for the effects caused by the deletion of the PxSODC gene, implying that the SOD protein encoded by PxSODC plays a crucial role in O₂^- scavenging.

Figure 6.

PxSODC-allied metabolic networks

To explore the metabolic networks associated with the PxSODC gene in P. xylostella and better understand the biological functions underlying their complex relationships, we performed an untargeted metabolomic analysis on the ancestral and mutant (SODC-MU1 and SODC-MU2) strains at different developmental stages. We detected 608 metabolites across the strains, including 167 organic acids and their derivatives, 102 nucleotides and their metabolites, 102 lipids, 92 amino acids and their metabolites, 33 heterocyclic compounds, 24 carbohydrates and their metabolites, and 88 other metabolites (Supplemental Figure S8A). Principal component analysis (PCA) revealed a distinction between the ancestral and mutant strains along the first component (PC1). However, a trend of clustering was observed in the 3^rd-instar larvae of ancestral and SODC-MU1 strains and the male adults of ancestral and SODC-MU2 strains, suggesting that the difference at these developmental stages may be driven by a few key metabolites (Supplemental Figure S8B). Compared to ancestral strain, different metabolites were identified in the eggs, 3^rd-instar larvae, pupae, female adults and male adults of mutant strains (Supplemental Figure S8C). The mutant strains had 68, 103, 110, 77 and 35 common differential metabolites compared to the ancestral strain at different developmental stages of P. xylostella (Supplemental Figure S8D). In addition, the lipid content significantly decreased in the 3^rd-instar larvae, pupae and female adults, while the content of nucleotides and their metabolites increased in the eggs, 3^rd-instar larvae and pupae following deletion of the PxSODC gene in P. xylostella (Supplemental Figure S9A). The up-regulation and down-regulation of various metabolites may be associated with the antioxidant function of PxSODC. We observed that the deletion of PxSODC gene could lead to elevated levels of the O₂^- within P. xylostella (Figure 6I).

We then performed enrichment analysis of differential metabolite KEGG pathways for the SODC-MU (MU1 and MU2) eggs, 3^rd-instar larvae, pupae, female adults and male adults resulting in these differential metabolites wre distributed in 23, 28, 28, 22, and 18 pathways, respectively. They were mainly involving lipid metabolism, nucleotide metabolism, carbohydrate metabolism, cofactor and vitamin metabolism, and amino acid metabolism (Supplemental Figure S9B). These pathways were related to the PxSODC gene and might have contributed to the temperature-adaptive capacity in P. xylostella through the regulation of biological functions/processes in different temperature environments.

In the metabolome, we observed that the abundance of 5-hydroxymethyluracil, 2-methylcitric acid and 5’-deoxyadenosine differed in all developmental stages of mutant strains (SODC-MU1 and SODC-MU2) compared to the ancestral strain (Figure 7A). Apart from a significant increase in 5’-deoxyadenosine in the 3^rd-instar larvae of mutant strains compared to the ancestral strain, the levels of these three metabolites were significantly reduced in other developmental stages, suggesting a possible direct association with PxSODC (Figure 7B) and representing a key biological regulatory response adopted by P. xylostella in adaptation to different environmental conditions. The significant roles of 5-hydroxymethyluracil (involved in dynamic methylation regulation) and 5’-deoxyadenosine (as a precursor to methyl donors) in DNA methylation has been previously confirmed (Bhutani et al., 2011; McKean et al., 2019). We measured the expression levels and enzymatic activities of DNA methyltransferase 1 (PxDnmt1) in the 3^rd-instar larvae of different strains. The results showed that both the expression and enzymatic activity of PxDnmt1 were significantly reduced in the hot and cold strains compared to the ancestral strain (Figure 7C, 7D). Using RNA interference (RNAi) technology, we specifically silenced the expression of PxDnmt1 in the ancestral strain of P. xylostella (Figure 7E) and observed significantly reduced levels of 5-methylcytosine (5-mC, a marker of DNA methylation) in both pupae and female adults (Figure 7F) (Ni et al., 2023). Further, we found that silencing of PxDnmt1 could significantly decrease the critical thermal maximum (CTmax) of female adults and increase the supercooling and freezing points of pupae (Figure 7G, 7H). These results confirm that DNA methyltransferase is associated with the adaptive response of P. xylostella to extreme temperatures.

Figure 7.

Discussion

Insects are valuable bio-indicators of the effects of climate change via their phenology, distribution and species interactions, population dynamics and community structure responses (Bale et al., 2002; Halsch et al., 2021). The present study demonstrates that P. xylostella can undergo rapid genetic adaptation to thermal extremes, with evolved changes spanning transcriptomic, metabolomic, and life history traits. By conducting the life history trait, demography, and fitness assessment of the ancestral, hot and cold strains of P. xylostella, we observed that the hot and cold strains had evolved significant genetic differences from the ancestral strain in multiple traits related to thermal tolerance and fitness. In the context of global warming, P. xylostella may evolve greater flexibility across its range and ecological niche leading the hot-evolved populations to be able to persist in regions with increased temperatures due to climate change (Chen et al., 2021). With climate change, cold-adapted P. xylostella may be favored during episodes of late frost in the spring or early frost in the autumn in temperate regions. The lower supercooling and freezing points of the cold strain pupae facilitate the survival of insects in cold climates, extending their ecological adaptability to low-temperature environments (Block, 1997). The demonstrated capacity of P. xylostella to adapt to extreme thermal conditions of both forms implies that P. xylostella may survive under a broader range of climatic conditions, posing new challenges for the management and control of this worldwide pest.

Our findings reveal a significant metabolic adjustment in P. xylostella during its adaptive evolution to both high and low temperatures. Lipids such as octadecenoic acid, epoxystearic acid and carnitine are reduced in the hot and cold strains compared to the ancestral strain, which indicates that after exposing and adapting to novel temperature environments, P. xylostella may maintain basic metabolic rates by lowering the level of certain metabolites, thereby allocating more energy to development and reproduction (Mallard et al., 2018). Lipids play a key role in energy storage and membrane stability in insects; adjusting the lipid metabolism may be a common mechanism for insects to adapt to temperature changes (Rommelaere et al., 2019). Carnitine has important functions in fatty acid transport and energy metabolism; change in its levels may impact energy supply (Bremer, 1983). Although gene expression between the hot and cold strains does not show a strong correlation with metabolites, KEGG analysis indicates that differentially expressed genes between the ancestral strain and the hot or cold strains have similarly enriched pathways. More differentially expressed genes, but fewer differential metabolites, were found in the cold strain compared to the ancestral strain than in the hot strain, a phenomenon different from that observed in Drosophila (Mallard et al., 2018). This might be due to the species-specific difference or the difference in selection pressure in the cold environment and requires further investigation.

Integrating information from metabolomic and transcriptomic analyses revealed extensive changes in gene expression and metabolic rewiring, suggesting that P. xylostella could adopt similar adaptive evolutionary strategies in hot and cold environments. Although the specific genes involved in the adaptation may be different, their functional pathways are similar, especially the important role of lipid metabolism in temperature adaptation, possibly affecting the environmentally sensitive phenotypes of heat tolerance, cold hardiness, and female fecundity (Tigano et al., 2020; Sherpa et al., 2022; Li et al., 2024a).

Insect genomes contain genes encoding the proteins involved in thermal adaptation (Kelley et al., 2014). Mutations in these genes during adaptive evolution in response of thermal adaptation can lead to novel phenotypic traits, such as changes in the life history and population fitness of insects (Gibson et al., 2019). The PxSODC gene encodes a superoxide dismutase that scavenges superoxide anions in cells, maintaining redox balance and protecting cellular structures (Sheng et al., 2014). Under extreme temperatures, insects can adjust their survival strategy, allocating more energy to maintain fundamental life functions (Hoffmann and Sgro, 2011). Here, the deletion of the PxSODC gene led to reduced SOD enzyme activity and increased O₂^- levels, reducing the tolerance of P. xylostella to extreme temperatures. These findings provide novel insights into how genetic variation translates into phenotypic variation (You et al., 2024), and the ways in which P. xylostella responds – and has responded – to a changing climate (Chen et al., 2021).

Expression of the PxSODC gene in both hot- and cold-evolved strains was lower than that in the ancestral strain, and the two evolved strains also had lower O₂^- levels, suggesting that they should have adapted to extreme temperature environments. This adaptation may result from non-synonymous mutations that alter the conformation of the PxSODC-encoded protein, increasing its capacity to scavenge superoxide anions and reducing the demand for SOD enzymes. A previous study has also shown that mutations in the Mn-SODC protein at positions K182R and A183P can increase its enzyme activity (Sheng et al., 2013). After long-term adaptation, insects may acquire the ability to maintain cellular homeostasis in new thermal environments by reducing basal metabolism and allocating more energy to development and reproduction (Mallard et al., 2018). At 34°C and 36°C, the trends in SOD enzyme activity and PxSODC gene expression differ between different temperature-evolved strains, suggesting the involvement of additional genes in the regulation of SOD enzyme activity. This hypothesis is supported by our transcriptomic analysis identified additional SOD genes (Supplemental Table S5).

In the SODC gene family, two SOD genes (Px04C00505 and Px13C00423) underwent non-synonymous mutations and showed reduced expression in the hot and cold strains at different developmental stages, while the expression of three SOD genes without mutations remained relatively stable or increased at most developmental stages. The evolution of protein functions is driven by mutations, which in some cases can switch directly from one function to another through single amino acid changes (Nobeli et al., 2009; Patsch et al., 2024). Under extreme temperatures, the expression of these genes in the 3^rd-instar P. xylostella larvae of the hot and cold strains trends to be similar to those of the ancestral strain. This indicates that adverse environmental conditions increase intracellular oxidative stress, which requires regulation of SOD expression and enzyme activity to scavenge superoxide anions (Islam et al., 2022). Maintaining high levels of SOD enzyme activity requires additional energy, placing a strain on cellular energy metabolism and resource allocation (Emre et al., 2013). SOD genes with non-synonymous mutations, similar to PxSODC, can lead to the change in protein structure or function, affecting enzyme activity and allowing for faster O₂^- scavenging at lower transcript levels, reducing resource requirements. The three unmutated SOD genes, if mutated, might adversely affect the moth or have other functions, such as involvement in cellular signaling pathways (Mondola et al., 2016).

Metabolomic analysis of different developmental stages in the ancestral and mutant strains revealed that after the loss of PxSODC gene, the metabolism of P. xylostella underwent temperature-adaptive adjustments involving lipids, nucleotides, carbohydrates, coenzymes and vitamins, and amino acids. This study also showed that DNA methylation plays a key role in the temperature adaptation of P. xylostella. While DNA methylation may be associated with gene activation, its main function remains the inhibition of gene expression (Stroud et al., 2015; Wang et al., 2018). Transcriptomic analysis also showed that more genes were up-regulated in the hot and cold strains compared to the ancestral strain, highlighting the role of DNA methylation in regulating gene expression to help P. xylostella maintain physiological functions and survive. In addition to directly regulating the expression of temperature responsive genes, epigenetic effects can also indirectly affect the response of insects to temperature challenges (Reynolds, 2017). By adding the molecular data of epigenetic markers, underlying mechanisms on the adaptive adaptation can be more easily elucidated. Therefore, further work is required to better understand the impact of non-genetic effects on adaptation to future climates including how they interact with genetic adaptive capacity.

This study elucidates the molecular mechanisms underlying adaptation of Plutella xylostella to both high- and low-temperature environments and functionally validates differentially expressed genes identified in ancestral, hot-evolved, and cold-evolved strains. Nevertheless, thermal adaptation in arthropods may engage distinct, temperature-specific biological pathways; accordingly, future work will prioritize the characterization of strain-specific differentially expressed genes. Beyond functional validation of the canonical stress-associated gene PxSODC, additional genes harboring nonsynonymous mutations warrant detailed investigation to clarify their roles within the broader regulatory network. Importantly, our findings also underscore a critical role for DNA methylation in thermal adaptation in P. xylostella. Building on this insight, we plan to generate PxDnmt1 knockout lines and apply whole-genome bisulfite sequencing (WGBS) across thermal strains to identify key methylation sites, thereby enabling a more comprehensive assessment of their contributions to temperature adaptation in P. xylostella and the generality of these mechanisms across other species.

Materials and methods

Insects

The founding population of P. xylostella was established from field-collected larvae on cabbage (Brassica oleracea var. capitata) in organic farms in Fuzhou, Fujian Province, China (26°05’N, 119°18’E) in July 2012. The collection site was confirmed through farmer interviews and local agricultural records to have no history of insecticide application for at least five years. Approximately 300-500 larvae were collected from multiple host plants distributed across a 2-hectare area to maximize genetic diversity and minimize founder effects. The field-collected population was reared in the laboratory under controlled conditions of 26°C and 60% relative humidity with a 12-h light/12-h dark cycle, without exposure to insecticides. This setup was referred to as the favorable temperature environment based on a previous study on the relationship between temperature and the P. xylostella development (Liu et al., 2002). Population size was maintained at >500 individuals per generation to minimize inbreeding and genetic drift. This laboratory-acclimated population was designated as the ancestral strain and served as the baseline for all subsequent experiments. Eggs and larvae were reared in sterile plastic disposable Petri dishes (90 mm) with an artificial diet containing 68 g yeast powder, 20.4 g agar, 127.5 g raw wheat germ, 3.4 g potassium sorbate, 3.4 g methyl paraben, 34 g sucrose, 10.2 g powder of radish seed, 1.7 g vitamin premix, 3.4 g ascorbic acid, 3.4 mL cola oil, 0.34 mL linoleic acid and 850 mL water The adults were allowed to mate and lay eggs in disposable paper cups, and were fed with a 10% honey solution.

The previous study on the relationship between temperature and developmental rate shows that P. xylostella can survive and develop at the temperatures 32°C at the maximum, 26°C as the optimum, and 10°C for the minimum (Liu et al., 2002). To generate thermally adapted populations, we established 18 replicate populations from the ancestral strain and randomly assigned them to three thermal regimes: (1) a hot-evolved treatment with cycling temperatures of 32°C/27°C (12 h light/12 h dark), (2) a cold-evolved treatment with cycling temperatures of 15°C/10°C (12 h light/12 h dark), and (3) a control treatment maintained at constant 26°C. All other environmental conditions (humidity, photoperiod, diet) remained identical across treatments. The populations were maintained with non-overlapping generations and a census population size of approximately 500 individuals per replicate.

The thermal selection experiment was conducted for approximately three years. To minimize maternal effects and ensure that observed differences were due to genetic adaptation rather than developmental plasticity, all populations were reared for two additional generations under common garden conditions (26°C) prior to phenotypic and molecular assays. For each treatment, we randomly selected six replicate populations for downstream experiments, resulting in 18 experimental cohorts in total.

Development of age-stage-specific two-sex life tables

To study the life history and population fitness of the three different P. xylostella strains, age-stage-specific life tables were developed. A total of 90 newly laid eggs from each strain were transferred to three 90 mm diameter Petri dishes (30 eggs per dish) (Chi et al., 2020).

The stage-specific number of individuals was recorded daily (food replaced every two days). During the pupal stage, each pupa was placed in a perforated 1.5 mL centrifuge tube. After eclosion, adults were transferred to 50 mL plastic cups for mating and oviposition and fed with a 10% honey solution. We continuously monitored and recorded daily oviposition and the number of surviving adults within the population until all individuals died.

For the individual life tables, 120 randomly selected newly laid eggs from each strain (i.e., 120 eggs/strain/temperature) were placed under different temperature conditions (favorable, hot, and cold environments). The eggs were individually transferred to 30 mm diameter petri dishes (artificial diet replaced every two days). The number of surviving larvae and pupae was recorded daily, and each pupa was placed individually in a perforated 1.5 ml centrifuge tube. One newly emerged adult female and one newly emerged adult male were placed in a 25 mL plastic cup and fed with a 10% honey solution for mating and oviposition. The number of eggs laid by each adult female was recorded daily until death. The longevity of both males and females was recorded.

Life history and population fitness parameters were calculated to assess the life history traits of P. xylostella evolved in different temperature environments based on the recorded data of age-stage-specific two-sex life tables of different strains. The numerical computation was done using the TWOSEX-MSChart software (http://140.120.197.173/ecology/prod2.htm) and the bootstrap method with 100,000 replications to obtain standard errors of the fitness parameters. The paired BT method was also used to assess the difference between the different treatments. A P value of less than 0.05 indicates a statistically significant difference (Chi, 1988).

Metabolomic and transcriptomic profiling

Samples collected for metabolomic profiling included: (1) the 3^rd-instar larvae from different strains (AS, HS and CS) maintained at the favorable temperature (26°C); and (2) the 2-day-old eggs, 1-day-old 3^rd-instar larvae, 2-day-old pupae and newly emerged adults from the ancestral and SODC-MU (MU1 and MU2) strains (see section Deletion of the targeted gene using the CRISPR/Cas9 genome editing) at the favorable temperature. Six biological replicates were used for metabolomic profiling and three replicates for transcriptomic profiling for each of the strains.

Metabolomic profiling

Stored samples were thawed on ice and weighed (50 ± 2 mg) into 1.5 mL centrifuge tubes, to which three pre-cooled steel balls (3 mm) and 500 µL of pre-cooled 70% methanol (Merck, Germany) were added. Each of the samples was homogenized in a pre-cooled tissue homogenizer (25 HZ, 5 min) (Tissuelyser, Qiagen), and then the homogenized sample was left on ice for 15 min, then centrifuged at 12,000 rpm at 4°C for 10 min. The supernatant was transferred to a new 1.5 mL centrifuge tube and stored overnight at −20°C. The following day, the samples were centrifuged at 12,000 rpm for 3 minutes at 4°C. The supernatant was then collected with a sterile syringe and filtered through a 0.22 µm filter (Waters, USA) into an HPLC sample vial. The instrumental system for data acquisitions mainly used ultra-high performance liquid chromatography and tandem mass spectrometry (multiple reaction monitoring mode). The chromatography and tandem mass spectrometry conditions were as described by Li et al. (2021) Significant metabolites were identified using the criteria: |log₂Fold Change| ≥ 1; VIP ≥ 1; P < 0.05. All differential metabolites were assigned to different pathways by KEGG analysis.

Transcriptomic profiling

mRNA libraries were constructed for each of the samples and sequenced on the Illumina HiSeq4000 platform (Illumina, San Diego). Raw reads obtained from sequencing were filtered, low quality reads were removed using adapters, and clean reads were obtained for subsequent information analysis. Clean reads were aligned to the P. xylostella genome using HISAT2 (V2.1.0) (http://121.37.197.72:5010), with sequence alignment performed using the software’s default parameters. Gene expression levels were measured using FPKM (fragments per kilobase of transcript per million fragments mapped). The P-value was corrected for multiple hypothesis testing following Benjamini and Hochberg (1995). Differential expression analysis between samples was performed using DESeq2. Differentially expressed genes were identified using the criteria of |log₂Fold Change| ≥ 1 and FDR (false discovery rate) < 0.05. All differentially expressed genes were classified into different pathways by KEGG analysis.

Weighted gene co-expression network analysis (WGCNA)

In the transcriptome data set, genes with an FPKM value below 0.1 were filtered out. The WGCNA package in R (v2.0) was used to calculate weight values and a soft threshold was determined based on the scale-free network principle. A gene clustering tree was constructed based on gene expression correlations, and gene modules were identified based on these clustering relationships. The minimum number of genes per module was set to 30, and the cut height threshold was set to 0.25 to merge potentially similar modules, with other parameters set to default values. Trait data, comprising differential metabolites between the hot and cold strains relative to the ancestral strain, were used for correlation analysis. Correlation was considered to be significant if the absolute correlation coefficient was greater than 0.8 and P < 0.05. Genes within specific modules were compared with differentially expressed genes and those common to both were considered as candidate genes.

RNA extraction and cDNA synthesis

Total RNA was extracted using the Eastep® Super Total RNA Extraction Kit (Promega, Shanghai) according to the manufacturer’s instructions. RNA integrity and quality were assessed using a NanoDrop 2000 spectrophotometer (GE Healthcare, USA) and 2% agarose gel electrophoresis. Total RNA (2000 ng) was reverse transcribed into cDNA using the Reverse-Transcription System Kit (Promega, Shanghai).

Gene cloning

The following gene sequences were amplified: (1) differentially expressed genes identified from specific modules; and (2) genes from the PxSODC gene family. Reference sequences for these genes were obtained from the P. xylostella genome database and full-length primers were designed using Primer Premier 6.0 (Supplemental Table S7). PCR amplifications were performed using cDNA from the ancestral, hot and cold strains of P. xylostella (at least six samples) as templates, using Phanta® Max Super-Fidelity DNA Polymerase (Vazyme, China).

The 50 µL PCR reaction mixture contained 2× reaction buffer (25.0 µL), dNTP mix (1.0 µL), 20 µM upstream primer F (2.0 µL), 20 µM downstream primer R (2.0 µL), DNA polymerase (1.0 µL), nuclease-free water (15.0 µL) and cDNA (4.0 µL). Amplification conditions were as follows 95°C for 3 minutes, followed by 35 cycles of 95°C for 15 seconds, 58°C for 15 seconds and 72°C for a gene-specific extension time, with a final extension at 72°C for 5 minutes. The integrity of the PCR products was verified by 2% agarose gel electrophoresis and purified using the Omega gel extraction kit (USA). The purified fragments were cloned into the pESI-Blunt vector using the Hieff CloneTM Zero TOPO-Blunt Simple Cloning Kit (TIANGEN, China) and transformed into competent DH5α cells (YEASEN, Shanghai). The ligated product was transferred to DH5α competent cells and plated on LB+100 μg/ml ampicillin plates and incubated overnight at 37°C. A single colony was picked and placed in 500 µL of liquid LB with 100 µg/mL ampicillin. A positive clone was sent to Sangon Biotech (Shanghai, China) for sequencing. Sequence alignments of candidate genes from different strains were performed using Snap Gene software.

Sequence analysis and phylogenetic tree construction

To determine the exons and introns of the target genes, their sequences were aligned with gDNA from the P. xylostella genome database. Protein sequences were analyzed using a protein sequence analysis and classification tool (InterPro, http://www.ebi.ac.uk/interpro/). The relative molecular masses and isoelectric points of the proteins were predicted using Expasy (http://expasy.org/tools/dna.html). The secondary structures of the proteins were predicted using PSIPRED (http://bioinf.cs.ucl.ac.uk/psipred). A phylogenetic tree was constructed using the neighbor-joining method in MEGA X software, and the robustness of the tree was verified by bootstrap analysis (bootstrap = 1000 replicates).

qRT-PCR analysis

For real-time quantitative PCR analysis, we collected samples as follows: (1) we randomly selected 22 genes from the transcriptomes of the 3^rd-instar larvae of ancestral, hot, and cold strains to validate the transcriptome data; (2) we collected samples from the eggs, 3^rd-instar larvae, pupae and adult males and females of the ancestral, hot, cold and SODC-MU (MU1 and MU2) strains to assess the transcription levels of genes including PxSODC, Px04C00505, Px13C00423, Px20C00248, Px15C00224 and Px15C00223; (3) we collected the 3^rd-instar larvae of the ancestral and hot strains exposed to different high temperature treatments (32°C, 34°C, 36°C) to analyze the transcription levels of the above-mentioned genes; and (4) we collected the 3^rd-instar larvae of the ancestral and cold strains exposed to different low temperature treatments (12°C, 10°C, 8°C) for similar assessments.

qRT-PCR primers were designed using Primer Premier 6.0, with PxRPL³²as the reference gene (Supplemental Table S7). The reaction mixture of 20 μL contained 10 μL 2× Taq Pro Universal qPCR Master Mix (Vazyme, Nanjing, China), 0.4 μL of each primer, 7.15 μL nuclease-free water and 2.0 μL cDNA. The QuantStudio Real-Time PCR System (Thermo, USA) protocol was as follows: 95°C for 30 s; 40 cycles of 95°C for 10 s, 60°C for 30 s; followed by melting curve analysis of 95°C for 15 s, 60°C for 1 min, 95°C for 15 s. Each sample contained 3 biological replicates and 3 technical replicates, and the relative expression of genes was calculated using the 2^-ΔΔCt method.

Deletion of the targeted gene using the CRISPR/Cas9 genome editing

We generated stable homozygous mutant strains of P. xylostella with the PxSODC gene deleted using the CRISPR/Cas9 system. The target site was designed based on the 5’-N20NGG-3’ motif (underscore indicates PAM sequence), and the potential off-target effect of sgRNAs was predicted using Cas-OFFinder (http://www.rgenome.net/cas-offinder). The in vitro transcription template for sgRNA was generated from a single nucleotide strand under the following conditions: 95°C for 3 min, followed by 35 cycles of 95°C for 15 s, 68°C for 15 s and 72°C for 30 s, with a final extension at 72°C for 5 min. The amplified product was purified by gel extraction. The sgRNA was obtained by in vitro transcription of the gel-purified product using the HiScribe™ T7 Quick High Yield RNA Synthesis Kit (New England Biolabs, USA). The reaction mixture contained 2.5 μL NTP buffer mix, 0.5 μL T7 RNA polymerase mix, 65 ng gel-purified product, made up to 5 μL with nuclease-free water. After overnight incubation at 37°C, 0.5 μL of DNase was added to remove DNA and the product was incubated at 37°C for 20 minutes to yield sgRNA. The sgRNA was purified by phenol-chloroform extraction and stored at −80°C.

We prepared a 10 μL reaction mixture containing 300 ng/μL sgRNA and 200 ng/μL Cas9 protein (GenCrispr, Nanjing), 1 μL 10× reaction buffer and nuclease-free water to make up to 10 μL and incubate at 37°C for 25 minutes. The mixture was injected into freshly laid eggs using the Olympus SZX16 microinjection system (Olympus, Japan) and the entire microinjection was completed within 30 minutes of eggs being laid. After injection, the eggs were placed in a Petri dish and the number of eggs hatched was recorded. Adult gDNA was extracted using the TiANamp Genomic DNA Kit (TIANGEN, China). Specific primers were designed for PCR amplification (Supplemental Table S7), with conditions as follows: 95°C for 3 min, followed by 34 cycles of 95°C for 15 s, 58°C for 15 s and 72°C for 15 s, with a final extension at 72°C for 5 min. The sequence of the PCR products was checked by Sangon Biotech (Shanghai) Co., Ltd.

The injected eggs were referred to as the G0 generation. These were reared to adulthood, crossed with the ancestral (non-injected) adults and used to extract genomic DNA from G0 adults after oviposition (the resulting progeny representing the G1 generation). PCR products flanking the two sgRNA target sites were amplified as mentioned above to determine genotypes and identify heterozygotes (individuals with double peaks in the sequence chromatogram starting from the sgRNA target site). The G1 generation was self-crossed to produce the G2 generation, and all G1 adults were genotyped based on PCR amplification for individual identification. G2 progeny derived from G1 heterozygotes with the same allelic mutation were selected. The G2 generation was then self-crossed to produce the G3 generation, retaining those with the same type of homozygous mutations to establish homozygous lines. If the G3 generation remained heterozygous, self-crossing continued until homozygous mutations were obtained (Wang et al., 2020). By the end, three mutants were obtained and called SODC-MU (MU1, MU2 and MU3) strains.

RNA interference mediated silencing of target genes

Gene-specific primers containing T7 promoter sequences were designed (Supplemental Table S7), and PCR was performed using total P. xylostella cDNA as a template. The PCR products were purified using a gel extraction kit. Double-stranded RNA (dsRNA) was synthesized by in vitro transcription using the TranscriptAid T7 High Yield Transcription Kit (Vazyme, Nanjing, China). The dsRNA was diluted to 2 μg/μL using DEPC-treated water (Beyotime, Shanghai, China). A volume of 500 nL diluted dsRNA (dsEGFP or dsDnmt1) was injected into pupae using a Nanoliter 2000 microinjector (World Precision Instruments LLC, USA). Total RNA was extracted 24 hours after injection and reverse transcribed to cDNA. Gene knockdown efficiency was analyzed by qPCR using beetles injected with dsEGFP as controls. The experiment was performed in three independent biological replicates (Zhou et al., 2024).

Assessing the response to high temperature

To assess the response of different P. xylostella strains (ancestral strain, hot strain and mutants) to extremely high temperatures, 2-day-old eggs, 1-day-old 3^rd-instar larvae and 2-day-old pupae were individually placed in 90 mm diameter Petri dishes. Adult females and males were placed individually in perforated 1.5 mL centrifuge tubes. Based on a preliminary trial on the stage-specific temperature tolerance limit of P. xylostella (eggs, larvae, pupae, and both male and female adults of the ancestral and hot strains were placed in different temperature environments ranging from 40 to 45°C), pupae from the ancestral and mutant strains were exposed to 43°C while eggs, larvae, and adults were exposed to 42°C for periods ranging from 30 to 180 minutes. After treatment, all replicate samples were transferred to an environment maintained at 26°C, where survival was observed and recorded. Survival was defined as the successful development of eggs, larvae and pupae to the next stage, while adults had to show movement of an appendage or mouthparts. Experiments were performed with six biological replicates, with each replicate contained 20 individuals.

We randomly selected 20 female adults injected with dsRNA to determine their critical thermal maximum (CTMax). A thermistor probe (Omega, USA) was inserted into a 1.5 mL centrifuge tube, which was suspended inside a 50 mL centrifuge tube with the opening sealed with cotton. This assembly was then placed in a 2 L glass beaker containing 1000 mL water, with the beaker top sealed with insulating foam board. The entire setup was positioned on a thermostatically controlled magnetic stirrer, where the temperature inside the 1.5 mL tube was increased at a constant rate of 0.5°C/min. When the temperature reached 26°C, female adults were quickly transferred into the 1.5 mL centrifuge tube containing the temperature probe, and their behavioral responses were continuously monitored as temperature increased. The CTMax was recorded when moths exhibited spasms, lost their crawling or flying ability, and remained motionless at the bottom of the tube, typically lying ventral side up (in most cases) or dorsal side up (in fewer instances). Although antennae and limbs might still exhibit slight tremors at this point, the insects typically died within seconds (Li et al., 2024b).

Measurement of the supercooling and freezing points

To investigate the cold hardiness of different P. xylostella strains (including the ancestral, cold and mutant trains), we randomly selected 40 pupae from each strain to examine their supercooling and freezing points. A thermistor probe from a subcooling point tester (Omega, USA) was attached to a pupa, secured with conductive tape and placed in a centrifuge tube, with the tube mouth sealed by cotton. The centrifuge tubes were then placed in a 50 mL plastic cup filled with ethylene glycol (antifreeze), and the cup was stored in an ultra-low temperature freezer set at −70°C, with the temperature first dropping rapidly and then decreasing at a rate of 0.10°C per second until the supercooling point was reached. By recording temperature changes at intervals of every second, the supercooling and freezing points of pupae were determined based on the inflection point of body temperature. The same experimental approach was also applied to P. xylostella injected with dsRNA.

Detection of oxidative stress indicator

Samples were collected for the detection of superoxide dismutase (SOD) and superoxide anion (O₂^-) in the following conditions: (1) the eggs, 3^rd-instar larvae, pupae, and adult males and females of the ancestral, hot, cold, and SODC-MU (MU1 and MU2) strains at the favorable temperature (26°C); (2) the 3^rd-instar larvae of the ancestral, hot, and SODC-MU (MU1 and MU2) strains after 2 hours of heat stress at 32°C, 34°C, and 36°C; and (3) the 3^rd-instar larvae of the ancestral, cold, and SODC-MU (MU1 and MU2) strains after 2 hours of cold stress at 12°C, 10°C, and 8°C. The experiment was performed with four independent biological replicates. The levels of SOD and O₂^- were measured using commercial assay kits (Comin, China) according to the manufacturer’s instructions.

Detection of 5-methylcytosine concentration

Pupae and female adults were collected for detection of 5-methylcytosine (5-mC) concentration after injection with dsRNA. The levels of 5-mC were measured using a commercial insect 5-methylcytosine (5-mC) ELISA detection kit (Shanghai Enzyme-linked Biotechnology Co., Ltd., China) according to the manufacturer’s instructions. The experiment was performed with four independent biological replicates.

Data analysis

Statistical analysis was performed using SPSS version 23.0. First, the Shapiro-Wilk test was used to assess the normality of data distribution. For the normally distributed data, independent samples t-tests were used to analyze the difference between the control and treated groups. For the data involving three or more groups with homogeneity of variances, one-way ANOVA followed by Tukey’s multiple comparison test was used, and if variances were uneven, Tamhane’s T2 test was used. For the non-normally distributed data, the Mann-Whitney test was used to analyze the difference between two groups, or the Kruskal-Wallis test was used to compare the differences among three or more groups (Lei et al., 2024).

Data availability

The raw sequence data generated in this study have been deposited in the Genome Sequence Archive at the National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences, under accession number CRA024611 (https://ngdc.cncb.ac.cn/gsa). The associated metadata are available in the OMIX database under accession numbers OMIX009807 (https://ngdc.cncb.ac.cn/omix/preview/kkl093nl) and OMIX009846 (https://ngdc.cncb.ac.cn/omix/preview/gXwZCXMI).

Acknowledgements

This work was financially supported by the central government-guided local science and technology development project (2022L3087), the Fujian Natural Science Fund for Distinguished Young Scholars (2022J06013), and the research project grants of Fujian Agriculture and Forestry University (KFXH23021, KFB24114A), the State Key Laboratory of Agriculture and Forestry Biosecurity, the International Joint Research Laboratory of Ecological Pest Control, Ministerial and Provincial Joint Innovation Centre for Safety Production of Cross-Strait Crops, and the “111” Program in China.

Additional information

Author contributions

Conceptualization, G.L., H.Z., Y.C., F.Y., M.Y., and S.Y.; Data curation, G.L., H.Z., Z.M., Y.D., M.Y., and S.Y.; Formal analysis, Z.M. and Y.D.; Investigation, Z.M. and Y.D.; Methodology, G.L., H.Z., Z.M., Y.D., Y.C., F.Y., M.Y., and S.Y.; Project administration, S.Y.; Supervision, L.V., G.M.G., and S.Y.; Writing – original draft, G.L., H.Z., M.Y., and S.Y.; Writing – review & editing, G.L., H.Z., M.Y., L.V., G.M.G., and S.Y.

Funding

the central government-guided local science and technology development project (2022L3087)

• Shijun You

the Fujian Natural Science Fund for Distinguished Young Scholars (2022J06013)

• Shijun You

the research project grants of Fujian Agriculture and Forestry University (KFXH23021)

• Shijun You

the research project grants of Fujian Agriculture and Forestry University (KFB24114A)

• Shijun You

Additional files

Supplementary file 1. Tables S1-S7.

Supplementary file 2. Figures S1-S9.