# Comprehensive annotation of olfactory and gustatory receptor genes and transposable elements revealed their evolutionary dynamics in aphids

Sergio Gabriel Olvera-Vazquez, Xilong Chen, Aurélie Mesnil, Camille Meslin, Fabricio Almeida-Silva, Johann Confais, Yann Bourgeois, Gianluca Lombardi, Célia Lougmani, Karine Alix, et al.

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

## To cite this version:

Sergio Gabriel Olvera-Vazquez, Xilong Chen, Aurélie Mesnil, Camille Meslin, Fabricio Almeida-Silva, et al.. Comprehensive annotation of olfactory and gustatory receptor genes and transposable elements revealed their evolutionary dynamics in aphids. Molecular Biology and Evolution, 2025, 42 (12), pp.msaf238. 10.1093/molbev/msaf238. anses-05299488

HAL Id: anses-05299488

https://anses.hal.science/anses-05299488v1

Submitted on 8 Oct 2025

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

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

<figure><img src="image_2.png" /></figure>
Comprehensive annotation of olfactory and gustatory receptor genes and transposable elements revealed their evolutionary dynamics in aphids

Sergio Gabriel Olvera-Vazquez¹ a), Xilong Chen¹,² b), Aurélie Mesnil¹,² c), Camille Meslin³,
d), Fabricio Almeida-Silva⁴, e), Johann Confais⁵, f), Yann Bourgeois⁶, g), Gianluca Lombardi⁷, h),
Célia Lougmani¹, i), Karine Alix¹, j), Nicolas Francillonne⁵, k), Nathalie Choisne⁵, l), Stephane
Cauet⁸, m), Jean-Christophe Simon⁹, n), Christelle Buchard¹⁰, o), Nathalie Rodde⁹, p), David
Ogereau¹⁰, q), Claire Mottet¹¹, r), Alexandre Degrave¹², s), Elorri Segura¹¹, t), Alessandra
Carbone⁷, 13u), Barrès Benoît¹¹, v), Emmanuelle Jacquin-Joly³, w), William Marande⁸ x), Dominique
Lavenier¹⁴, y), Fabrice Legeai⁹, z), Amandine Cornille¹, 2, zz)

¹) Université Paris Saclay, INRAE, CNRS, AgroParisTech, GQE -Le Moulon, 91190, Gif-sur-Yvette, France

2) New York University Abu Dhabi, Abu Dhabi, United Arab Emirates

³) INRAE, Sorbonne Université, CNRS, IRD, UPEC, Université Paris Cité, Institut d'Ecologie et des Sciences de l'Environnement de Paris (iEES-Paris), 78026 Versailles, cedex, France

⁴) Department of Plant Biotechnology and Bioinformatics, Ghent University, 9052 Ghent, Belgium

5) Plant Bioinformatics Facility, BioinfOmics, INRAE, Université Paris-Saclay, 78026 Versailles, France

6) DIADE, University of Montpellier, CIRAD, IRD, 34394, Montpellier, France

7) Laboratory of Computational, Quantitative and Synthetic Biology, UMR 7238, Sorbonne University, Paris, France

8) Centre National de Ressources Génomiques Végétales-CNRGV, INRAE, Castanet, Tolosan, France

9) Institut de Génétique, Environnement et Protection des Plantes (IGEPP), Domaine de la Motte, 35653 Le Rheu Cedex - France

¹⁰) Laboratoire Évolution, Génomes, Comportement et Écologie, 1198 Gif-Sur-Yvette Cedex, France
1 $^{11)}$Université de Lyon, Anses, INRAE, USC CASPER, Lyon, France

2 $^{12)}$Université Angers, Institut Agro, INRAE, IRHS, SFR QUASAV, F-49000 Angers, France

3 $^{13)}$Institut Universitaire de France

4 $^{14)}$IRISA / INRIA. Campus de Beaulieu. 35042 Rennes cedex, France

5

6

7 a) solvera@koppert.com.mx

8 b) xc3350@nyu.edu

9 c) am15503@nyu.edu

10 d) camille.meslin@inrae.fr

11 e) fabricio.almeidasilva@psb.vib-ugent.be

12 f) johann.confais@inrae.fr

13 g) yann.bourgeois@ird.fr

14 h) gianluca.lombardi@sorbonne-universite.fr

15 i) c.lougmani@gmail.com

16 j) karine.alix@inrae.fr

17 k) nathalie.choisne@inrae.fr

18 l) nicolas.francillonne@inrae.fr

19 m) stephane.cauet@inrae.fr

20 n) jean-christophe.simon@inrae.fr

21 o) christelle.buchard@inrae.fr

22 p) nathalie.rodde@inrae.fr

23 q) david.ogereau@egce.cnrs-gif.fr

24 r) claire.mottet@anses.fr

25 s) alexandre.degrave@agrocampus-ouest.fr

26 t) elorri.segura@anses.fr

27 u) alessandra.carbone@lip6.fr

28 v) benoit.barres@anses.fr

29 w) emmanuelle.joly@inrae.fr

30 x) william.marande@inrae.fr

31 y) lavenier@irisa.fr
1 z) fabrice.legeai@inrae.fr

2 zz) amandine.cornille@nyu.edu

3

4 Corresponding author: amandine.cornille@nyu.edu

5

6

ACCEPTED MANUSCRIPT
# 1 Abstract

2 Gene duplication and transposable elements (TEs) are major drivers of genomic innovation that 3 can fuel adaptation. While the roles of duplication and TE-driven diversification are documented 4 in plant pathogens, they remain insufficiently explored in insect pests such as aphids, where 5 olfactory (OR) and gustatory receptor (GR) genes are key to host recognition. We analyzed 521 6 OR and 399 GR genes, alongside TEs, across 12 aphid genomes with varying host ranges. Aphid 7 lineages with broader host ranges exhibited higher evolutionary rates, driven by gene family 8 expansions linked to host interaction, including lipid metabolism, immune function, and 9 transposase activity. OR and GR genes evolved through proximal and tandem duplications and 10 were shaped by diversifying selection, with bursts of positive selection followed by prolonged 11 purifying selection, consistent with adaptation to novel hosts. Younger TEs were significantly 12 enriched near OR genes compared to GRs and other genomic regions, suggesting a catalytic role 13 of TEs in their diversification. However, OR proteins encoded by TE-associated ORs exhibited 14 reduced functional potential. In contrast, GR proteins encoded by TE-associated GRs retained 15 signatures of adaptation, as inferred from deep learning models predicting functionally important 16 protein regions. These findings suggest that TE activity may facilitate functional innovation in 17 GRs while alleviating constraints or pseudogenization in ORs. This study reveals how 18 duplication, selection, and TE dynamics shape gene evolution in insect pests. It also provides the 19 first chromosome-scale genome assembly of *Dysaphis plantaginea*, with comprehensive 20 annotations and functional predictions of OR/GR genes, bridging adaptive evolution with 21 mechanistic insights. 22

23

24 **Keywords:** insect, aphid, adaptation, chemoreception, speciation, crop, transposable elements, in 25 silico functional prediction.
# 1 Introduction

2 Studying the genomic bases of parasite adaptation to their host can contribute to understanding
3 the evolutionary processes underlying population adaptation (Gladieux et al. 2014; Simon et al.
4 2015; Moller and Stukenbrock 2017; Ebert and Fields 2020; Märkle et al. 2021; Cornille et al.
5 2022; Shih et al. 2023) and help tackle pressing issues such as parasite outbreaks (Kirk et al.
6 2013; Cacciò et al. 2018; Desvignes et al. 2022; Gowler et al. 2023). The questions related to
7 parasite adaptation to their host involve identifying the number and nature of genes involved, as
8 well as their mode of evolution, which can occur through drift (Barton and Mallet 1996; White et
9 al. 2021; Lynch 2023) or natural selection (Fay and Wittkopp 2008; Karlsson et al. 2014; Arnab
10 et al. 2023). Recently, the role of repetitive regions of the genomes - called transposable
11 elements (TEs hereafter) - in driving genomic variation at loci, fueling the arms race between
12 parasites and their hosts, has also gained attention (Duplessis et al. 2011; Hartmann et al. 2017;
13 Plissonneau et al. 2018; Rocafort et al. 2022; O'Donnell et al. 2025).

14 Gene duplication is a source of genomic variation that can fuel parasite adaptation to its
15 host (Wagner 1998; Magadum et al. 2013). When a newly duplicated gene is divergent, it can be
16 immediately favored by selection if the host cannot recognize it, improving the parasite's fitness.
17 The freshly duplicated gene is then fixed by positive selection (Kondrashov et al. 2002;
18 Kondrashov 2012). Periods during which genes evolve under positive selection that may occur
19 recurrently or episodically (also called diversifying selection) and be followed by prolonged
20 episodes of purifying selection (leading to an average rate of sequence evolution, dN/dS<1
21 (Kimura 1977; Kimura 1980)) on a long-term evolutionary time scale. Alternatively, a newly
22 duplicated gene can have redundant functions immediately, and one of the copies could
23 accumulate deleterious mutations and become a pseudogene (Assis and Bachtrog 2013). The role
24 of positive selection and purifying selection in driving the molecular radiation of genes involved
25 in the parasite adaptation to their host remains to be explored beyond model systems such as
26 viruses, bacteria, fungal pathogens, or insect models such as drosophila (McBride 2007; Assis
27 and Bachtrog 2013; McKenzie and Kronauer 2018; Remick et al. 2023; Lehman et al. 2024).
28
29 Studying the TE content and localization in parasite genomes is a way to understand the
30 evolutionary dynamics of genes involved in the parasite's adaptation to its hosts. The selfish
31 proliferation of TEs that drives their insertion dynamics can also drive the evolutionary dynamics
31 of duplicated genes involved in the parasite's adaptation to its host (Gilbert et al. 2021;
Oggenfuss et al. 2023). TEs are supposed to be catalyzers of proximal duplicated gene copies, i.e., genes near each other but separated by several genes (10 or fewer genes) (Liu and Wessler 2017; Dupeyron et al. 2019). In the TE-rich genomes of some plant pathogens, genome structure alternating gene-rich and TE-rich regions are observed and suggested to facilitate the generation of variation that fuels the arms race between parasites and their hosts (Seidl and Thomma 2017; Frantzeskakis et al. 2019; Schrader and Schmitz 2019; Frantzeskakis et al. 2020; Torres et al. 2020). While the role of TEs on parasite gene diversification is well documented in plant pathogens, their role in shaping the diversity of genes involved in insect-plant interactions is much less characterized (Stergiopoulos and De Wit 2009; Simon et al. 2015; Torres et al. 2020; Baril et al. 2023; Seong and Krasileva 2023). Yet, insects have a massive variation in TE and chemosensory gene content (Elliott and Gregory 2015; Gilbert et al. 2021).

Several chemosensory gene families play a key role in host recognition in herbivorous insects. Olfactory receptors (ORs) and gustatory receptors (GRs) are crucial signal-transducing receptors in insect olfaction and gustation (Jacquin-Joly and Merlin 2004; Vosshall and Stocker 2007; Vosshall and Stocker 2007; Engsontia et al. 2014). Other chemosensory genes, such as insect-binding proteins (IBPs) and odorant-binding proteins (OBPs), assist in transporting odorants (Pelosi et al. 2006). Since ORs and GRs trigger neural responses, they are essential targets for evolutionary studies. Additionally, functional data on ORs and GRs have shown their connection to speciation, making them central to understanding insect adaptation (Jacquin-Joly and Merlin 2004; Almeida et al. 2014; Engsontia et al. 2014). OR and GR duplicated genes form highly diversified families whose rapid evolution is driven by gene duplication (Niimura and Nei 2007; Kulmuni et al. 2013). Phylogenetic analyses reveal the origin of ORs and GRs. GRs emerged early in metazoan evolution and are the oldest chemosensory gene family (Robertson 2019; Robertson et al. 2019). In arthropods, OR genes likely originated from the diversification of the ancestral gustatory receptor (GR-like) family, as revealed by phylogenetic analyses (Robertson 2015; Eyun et al. 2017; Robertson 2019). However, the fate of OR and GR genes post-duplication is less documented and shows mixed patterns across insect groups (Hallem et al. 2006; McBride 2007; Joseph and Carlson 2015). In lepidopterans, gene gains and occasional losses appear linked to host adaptation, with polyphagous species such as *Helicoverpa armigera* Hübner and *Spodoptera frugiperda* showing GR expansions associated with broader host ranges (Xu et al. 2016; Gouin et al. 2017; Pearce et al. 2017). Evidence of positive selection in these
genes suggests host-driven adaptation, as seen in the *Spodoptera* species (Meslin et al. 2022).

Hemipterans, such as the pea aphid, show rapid OR and GR gene expansions alongside positive selection, supporting a model of adaptation by host specialization (Smadja et al. 2009; Smadja et al. 2012). Conversely, in dipterans like *Drosophila* sp., OR/GR diversity seems to be shaped more by genetic drift than selection, influenced by species endemism and smaller effective populations (Gardiner et al. 2008; Hansson and Stensmyr 2011). Therefore, research on the evolution of OR and GR genes still lags behind a few insect model species and is limited primarily due to technical challenges. High sequence divergence across and within species complicates the automated annotation of OR and GR genes, necessitating manual annotation and extensive curation to ensure accuracy. In addition, the functional validation of chemosensory genes, including OR and GR genes, remains experimentally challenging (Li et al. 2023), but in silico functional prediction methods using deep learning approaches offer a new avenue (Haley et al. 2025). Recent advances in deep learning approaches, such as MuLAN (Lombardi and Carbone 2024), enable in silico inference of residue-level functional potential directly from protein sequences. MuLAN leverages a local attention mechanism trained on protein-protein interaction datasets to predict functionally important residues, such as ligand-binding sites or structural interfaces. This opens the possibility of connecting molecular signatures of selection with their mechanistic and functional relevance, and then guiding experimental assays without requiring exhaustive testing of tens or hundreds of candidates.

The relationship between TEs and chemosensory gene evolution also remains an open question. TEs are prevalent across insect genomes and vary significantly in content and type among species (McCullers and Steiniger 2017; Gilbert et al. 2021; Ye et al. 2022; Baril et al. 2023). In some insects, high densities of TEs are linked to duplicated genes tied to adaptive traits, suggesting that TEs may promote gene duplications that support evolutionary responses to environmental pressures (Gilbert et al. 2021). For example, in the ant *Cardiocondyla obscurior*, OR genes are located in TE-dense regions, or “TE islands”, implying that TEs might facilitate OR gene duplication (McKenzie and Kronauer 2018). Similarly, TEs are found near insecticide-resistance genes in the pest *H. armigera*, suggesting that TEs may be instrumental in driving adaptive gene duplications (Klai et al. 2020). However, a detailed analysis of TE density and diversity surrounding GR and OR genes beyond a few model insects is still needed to fully understand the relationship between TEs and OR/GR evolution. TE-mediated tandem
duplications may significantly shape the evolution of OR and GR gene families, contributing to their adaptation to diverse host plants.

Due to their diverse host range, aphids are ideal for studying the genomic basis of host adaptation (Peccoud et al. 2010; Simon et al. 2015; Shih et al. 2023). OR and GR genes are key receptors for aphids, which use chemicals to screen their environmental landscape and detect their host (Robertson 2019; Robertson et al. 2019). The sequencing of several genomes of aphid species and the annotations of OR and GR genes in two aphid models, *Acyrthosiphon pisum* and *Myzus persicae* Sulzer, provided the first hints on the diversity and evolution of chemosensory genes and TE content in aphids (Smadja et al. 2009; The International Aphid Genomics Consortium 2010; He et al. 2015; Thorpe et al. 2018; Julca et al. 2020; Mathers et al. 2020; Biello et al. 2021; Byrne et al. 2022; Zhang et al. 2022). The role of TEs in the evolution of aphids has been investigated in only a few studies (Bouallegue et al. 2017; Ahmad et al. 2021; Baril et al. 2023). Comparative studies on the evolutionary dynamics of OR and GR genes across diverse aphid species with different host plants remain scarce. Additionally, the evolutionary dynamics of TEs in aphid genomes, particularly near OR and GR genes, is poorly understood.

This study examines the evolutionary forces that shape gene repertoire evolution in 12 aphid species with varying host plants. We focus on the OR and GR multigene families, particularly in the Aphidinae subfamily. Notably, the Aphidinae subfamily, which includes aphids with complex life cycles that colonize diverse host plants, provides a valuable comparative framework for understanding the evolutionary trajectories of OR and GR gene families during divergence. Additionally, we analyze TE diversity and dynamics to understand their role in OR and GR evolution. Using a newly assembled genome of the rosy apple aphid (Dysaphis plantaginea Passerini) and 11 publicly available aphid genomes from both the Aphidinae and Eriosomatinae subfamilies—along with one aphid-like outgroup species—we systematically annotated OR and GR genes and classified TEs using standardized pipelines. We address four key questions: (1) Do aphid species in the Aphidinae subfamily show expansion of the OR and GR gene family, or any other gene family that could be related to host adaptation? (2) What duplication mechanisms drive the evolution of OR and GR genes? (3) What types of selective pressures act on OR and GR genes post-duplication? (4) Do TE content and dynamics vary across aphid species, and more specifically, is there TE enrichment near OR and GR genes?
Beyond analyzing duplication mechanisms, selection regimes, and TE associations, we further incorporated MuLAN-based deep learning predictions to evaluate the functional potential of OR and GR residues under selection. This integrative approach offers a novel mechanistic OR and GR by linking evolutionary pressures with the predicted functional relevance of chemosensory receptors. Together, these analyses provide insight into the genomic bases driving aphid adaptation and speciation.

## Results

### Annotation of GR and OR genes in the aphid and aphid-like genomes

We retrieved 11 high-quality genome assemblies from the Aphidinae and Eriosomatinae subfamilies, along with one aphid-like species, *Daktulosphaira vitifoliae* Fitch (Phyloxerinae subfamily) (BUSCO completed gene values > 95.7%, Tables S1 and S2). All assemblies were retrieved on the BIPAA platform (https://bipaa.genouest.org/sp/). Additionally, we included a newly assembled genome of *D. plantaginea* (version 3; see details in Text S1, Tables S2 and S3). We then identified the OR and GR gene repertoire in the 13 selected genome assemblies.

Through a systematic homology search and manual curation across the 13 genomes (see details of the pipeline in Figure S1 and Text S2), we identified 668 OR and 580 GR genes. Among these, 521 OR and 399 GR genes were newly annotated in 10 aphid species (*A. gossypii*, *E. lanigerum*, *D. noxia*, *D. plantaginea*, *M. cerasi*, *M. persicae*, *P. nigronervosa*, *R. maidis*, *R. padi*, and *S. miscanthi*; Table S3). In the remaining three species (*A. glycines*, *A. pisum*, and *D. vitifoliae*), we confirmed 147 OR and 181 GR genes that were previously annotated, without identifying any additional OR or GR genes. On average, we detected 44 ± 16.85 GR and 51 ± 12.78 OR genes per species (Table S1, Figure S2), and the number of OR and GR genes was significantly lower only in *D. noxia* (*P*<0.001, Model 1, Table S4, Figure S3). Although not statistically significant (*P* > 0.05, Model 1, Table S5, Figure S3), the number of OR genes was higher than the GR genes. Additionally, there was no significant correlation between genome size and OR or GR gene repertoires (Spearman's correlation: $ρ_{GR}$ = 0.23, non-significant; $ρ_{OR}$ = -0.02, non-significant).

### Gene gains and losses, including OR and GR genes, in the Aphidinae subfamily
As a preliminary step before analyzing the evolutionary history of gene repertoires in aphids, we confirmed the genetic relationships among the studied species based on single-copy genes detected with OrthoFinder (Emms and Kelly 2019). From the 13 assemblies, OrthoFinder identified 5,868 orthogroups shared among all species, including 891 single-copy orthogroups (Figure S4 and Table S5). The phylogenetic tree constructed using these 891 single-copy orthogroups, representing 240,958 conserved amino acid residues and scaled with TimeTree, was highly statistically supported (Figure 1). The ultrametric tree provided a more comprehensive and well-supported phylogenetic representation, consistent with the known taxonomy of the species (Choi et al. 2018; Mathers et al. 2021; Byrne et al. 2022; Huang et al. 2023).

Subsequently, we employed CAFE (Mendes et al. 2021) to analyze the evolutionary dynamics of gene repertoires in the 12 aphid species (excluding *Daktulosphaira vitifoliae* because it is an aphid-like species), including the previously detected OR and GR genes. For CAFE, we used the ultrametric tree inferred earlier and the multi-copy orthogroups identified by OrthoFinder as inputs. OrthoFinder identified 28,672 multi-copy orthogroups (Figure S2 and Table S5). Following the CAFE pipeline's recommendations, we removed gene copies identified in more than one species or those with fewer than 100 gene copies to avoid bias in our inferences (using a customized Python script from the CAFE pipeline). As a result, we retained 22,948 orthogroups, which included 42 GR orthogroups (with 578 GR genes out of the 580 previously annotated) and 47 OR orthogroups (comprising 617 OR genes out of the 668 previously annotated) for further analyses.

We statistically analyzed the evolutionary rates (gene gain and loss) among gene families (orthogroups) across the phylogeny. We compared models assuming a single evolutionary rate along the phylogeny with models allowing evolutionary rates to vary based on taxonomic rank (genus, subfamily, family, species) or host preferences (details in Table S6). Among these, CAFE provided the most substantial support for the model with subfamily-specific evolutionary rates (log-likelihood -lnL = 243,475, ε = 0.02878, P-value < 0.05; Table S7). In this model, the Aphidinae subfamily exhibited the highest estimated gene gain and loss rate (birth-and-death rate, λ = 0.0092, P-value < 0.05, Figure 1, Table S7), which was 3.8 times greater than that of the Eriosomatinae subfamily (λ = 0.0024, P-value < 0.05, Figure 1). This result indicates that gene family turnover was higher in Aphidinae than in Eriosomatinae and Phylloxerinae (Figure 1).
We focused on 147 significantly expanded and 51 significantly contracted gene families (*P*-value < 0.01) in the Aphidinae subfamily, containing at least 50 genes, for gene ontology (GO) analysis. GO analysis of expanded gene families revealed pathways related to metabolism, transcription, lipid molecules and interactions, immune responses, and transposase activity (Figure S5, Table S8). The GO analysis of contracted gene families included general functions (Figure S6, Table S9). Notably, a large gene gain was observed in the clade of *M. persicae*, *M. cerasi*, and *D. plantaginea* clade (+1,148 genes). This clade comprises species that infest fruit trees, but also exhibit a wide range of other life-history traits, including variation in life cycle, reproductive strategies, and host specialization (see Table S6 for details).

Of the 42 GR orthogroups, 11 exhibited significant expansions, and five showed significant contractions across the 12 aphid species. Six of the 47 OR orthogroups had significant expansions and nine significant contractions. We identified four expanded chemosensory receptors (two GRs and two ORs) in the Aphidinae (Figure 1). The highest gain of GR and OR genes was observed in *A. pisum* and *M. cerasii* Fabricius, which are known for their broad host-plant range and ecological diversity (Blackman and Eastop 2008).

These findings highlight significant changes in OR and GR genes within a species subfamily known for its broad host range. However, their expansion and contraction occurred less than in other gene families associated with metabolism, transcription, lipid signaling, molecular interactions, immune responses, and transposase activity. To further understand the evolutionary dynamics of OR and GR genes, we investigated their origins and the selective pressures acting upon them to assess whether they evolved faster than the rest of the genome.

## Weak micro-synteny and evolution by tandem and proximal duplication of GR and OR genes

First, we examined the preserved order of genes on the same chromosome (Duran et al. 2009) and the arrangement between related species (also called micro-synteny). Second, we explored the mechanisms driving the evolution of OR and GR genes, focusing on gene duplication events, including proximal and tandem duplications, that contribute to the evolution of these families.

We investigated the micro-synteny of OR and GR protein sequences in the 12 aphid genomes with the syntenet R package (Almeida-Silva et al. 2023). A total of 216 out of 668 ORs
(33.8%) and 228 out of 580 GRs (39.5%) belonged to syntenic clusters (37 and 60 clusters, respectively) across the 12 aphid species. Syntenic clusters were, therefore, not deeply conserved across taxa (i.e., not specific species; Table S10). The GR genes were more syntenic than the OR genes (Table S10). Syntenet identified 70 and 31 genes with putative chemosensory functions as they were clustered with manually annotated GR and OR genes, some of which were clustered into orthogroups by OrthoFinder (Table S10). However, these genes were excluded from further analyses as they were not detected during the systematic manual annotation of OR and GR genes. Consequently, only those genes confirmed through manual annotation as OR/GR genes were considered for detailed analysis.

We investigated the mode of evolution of OR and GR genes with the R package doubletrouble (Almeida-Silva and Van de Peer 2025). We identified and classified them as dispersed duplicates (gene copies neither neighboring nor colinear), proximal duplicates (gene copies separated by 10 or fewer genes on the same chromosome), tandem duplicates (gene copies closely adjacent to each other in the same chromosome), and segmental duplicates (gene copy block with high sequence identity; (Qiao et al. 2019; Abdullaev et al. 2021)), depending on their copy number and genomic distribution. Most OR and GR genes originated from tandem and proximal duplications (Figure 2). For some species, most genes originated from “dispersed duplication,” meaning they could not be classified into segmental, tandem, or proximal duplications. This is often because dispersed duplicates were duplicated too long ago, and subsequent genomic rearrangements have obscured their original duplication mechanisms, such as tandem duplications that were later broken apart.

Therefore, the OR and GR gene families evolved mostly by tandem duplications and are not deeply conserved across taxa, but the GR genes are more conserved than the ORs. Our results suggested different evolutionary histories between OR and GR genes. We therefore further investigated selective pressures acting on OR and GR genes following tandem duplication.

## A faster evolution of OR and GR genes compared to genome-wide levels

We calculated the substitution rate of protein sequences (ω) and the ratio of nonsynonymous (dN) to synonymous (dS) substitutions for OR and GR genes, comparing them to the rest of the genome (hereafter referred to as OG for single-copy orthologous genes). This analysis aimed to
1 determine whether OR and GR genes evolved under neutral evolution, positive selection, or
2 purifying selection.

3 First, we defined the subfamilies of ORs and GRs to investigate the selective pressure.
4 We classified the subfamilies of OR and GR genes, building a Maximum Likelihood tree on 365
5 OR amino-acid sequences and 381 GR amino-acid sequences from our set of 12 aphid and one
6 aphid-like genomes (Figures S7 and S8, Table S11). The GR genes clustered into 12 subfamilies
7 (monophyletic groups with high bootstrap values > 80%, Figure S8), and the OR genes clustered
8 into 24 subfamilies (high bootstrap values > 80%, Figure S9). We detected two known GR
9 subfamilies, the sugar and fructose receptors (GR01 and GR02, respectively, Figure S8), based
10 on the similarities with previously reported ones in *A. pisum* and *A. glycines*. We also detected
11 the conserved odorant co-receptor (Orco or OR01, Figure S9) subfamily based on the clustering
12 with the already reported Orco genes of *A. pisum*, *A. glycines*, and *D. vitifoliae*. The high number
13 of OR subfamilies relative to GR subfamilies and the star-like topology of the phylogenetic tree
14 (Figure S9) suggested that OR genes are evolving faster than GR genes or may have diverged
15 more recently, as previously suggested (Robertson 2015; Eyun et al. 2017; Robertson 2019). We
16 explore these hypotheses in further detail below.

17 We used the codeml program in the PAML package v4 (Yang 2007) to estimate the
18 substitution rate of the protein sequences (ω) of the 12 and 24 GR and OR gene subfamilies,
19 respectively. Codeml computes the ratio (ω) for various models of sequence evolution (tree
20 branch or/and codon-site). The most likely model explaining the data is chosen based on a
21 likelihood ratio test computed among models. Here, we used the branch and site models. We did
22 not use the branch-site model as we had no previous information about the phylogeny of these
23 genes or an apparent chemosensory gene family clustering to test for selection. For most GR and
24 OR subfamilies, the free ratio (different ω values among branches, Yang, 2007) was the most
25 statistically supported model (Table S12), indicating that within each OR and GR subfamily,
26 sequences evolved with their ω (Table S12). One of the exceptions was the fructose receptor
27 genes (GR02), which displayed a single evolutionary rate among the sequences (Table S12). Of
28 the 381 GR sequences analyzed, 57.7% evolved under purifying selection, 21.62% under
29 positive selection, and 20.61% neutrally. Of the 365 OR sequences, 46.11% evolved under
30 purifying selection, 25.47% under positive selection, and 27.34% under neutral evolution. The
signatures of selection were not OR/GR subfamily-specific: each subfamily presented each of the three signatures of selection (Tables S13). However, lineage-specific models (branch models) have a limited ability to detect short episodes of positive selection that affect only a few amino acids or amino acids under recurrent diversifying selection. We, therefore, used site-specific models implemented in codeml on sequences of each GR and OR subfamily. The site-specific models can account for site rate variation and detect recent episodes of positive selection (diversifying selection). We used three pairs of site-specific models to test for recurrent, diversifying selection: M0 (one ratio) and M3 (Discrete), M1 (Neutral) and M2 (Selection), and M7 (Beta) and M8 (Beta & ω) and used the parameter cleadata=1 to remove sites containing gaps, to avoid bias in analysis. For the OR17 family, no sites remained for analysis when using this parameter setting in codeml, due to the high proportion of gaps in the alignment. The likelihood ratio test showed higher support for models assuming neutral or evolution by purifying selection (M0, M1a, and M7) models than the alternative ones supporting diversifying selection (M3, M2a, and M8) (Tables S14 and S15), except the gene groups GR03, GR04, GR06, and OR20 (Table S15). The posterior probability (PPs > 95%) of the Bayes empirical bases (BEB) distribution of the inferred positively selected sites under the M8 model further indicated 12 sites in the GR03, two sites in the GR04, one site in the GR06 and one site in the OR20 (these sites overlapped in M2a and M3; Table S16) under positive selection. Altogether, branch and site models indicate that OR and GR genes mainly evolved under purifying selection.

We extracted the nonsynonymous-to-synonymous substitution ratio ($\omega$ or $dN/dS$) from the best-fit models. We compared these values between GR and OR genes and orthologous genes across the entire genome. We used 50 randomly selected single-copy gene families identified by OrthoFinder to represent the rest of the genome, as single-copy genes provide a more accurate representation of background evolutionary rates, avoiding the potential confounding effects of gene duplication and copy number variation seen in multi-copy gene families. We ran the same pipeline explained above for GR and OR genes to estimate $\omega$ for each OG family. Fourteen percent of the single-copy orthologous gene groups evolved under the same $\omega$, while 86% fitted the branch models with a specific $\omega$ (Table S17). For the site models, all OG families (except one out of the 49) evolved under the same $\omega$ (M0), including 72% that evolved under purifying selection ($\omega$ < 1; Tables S18-S21). We then statistically compared the OG's $dN/dS$ ratio ($\omega$) to
those of the OR and GR gene subfamilies. The OR and GR genes exhibited significantly weaker purifying selection compared to single-copy genes (Model 1, Figure S10, Table S22). Our results indicate that GR and OR genes mainly evolve under purifying selection (Figure S10), with nonsignificant differences in selection pressure on OR genes compared to GRs. We also observed a negative correlation between the *dN/dS* ratio and divergence among species, indicating that genes with the lowest divergence were under the highest evolutionary rates (*p*-value<0.0004, r² = -0.09, Figure S11).

Purifying selection was the dominant signature observed in OR and GR genes, affecting over 50% of these genes. However, these signatures are notably weaker than those in other genomic regions (OG), indicating that OR and GR genes in aphids evolve at higher rates than the rest of the genomes. Additionally, respectively 4 and 37% of OR and GR genes show evidence of recent positive selection, together with the observed negative correlation between the evolutionary rates of GR and OR genes and species divergence, suggesting that these genes undergo rapid episodes of positive selection, followed by sustained periods of purifying selection, which aligns with a pattern of diversifying selection.

## Functional relevance of OR and GR genes

To strengthen the functional relevance of our evolutionary findings, we examined whether OR and GR genes, which show signatures of positive selection in the three codeml models—M1a vs. M2a, M7 vs. M8, and M0 vs. M3—also exhibit hallmarks of functional significance at the protein structural level. Three GR (GR03, 04, 06) and one OR (OR20) families consistently showed signatures of positive selection in the three codeml models (Table S15). Using site-specific attention scores estimated by the MuLAN functional prediction deep learning framework (Lombardi and Carbone 2024), we observed that positively selected sites identified via codeml in each of the four families consistently displayed significantly higher attention scores, except in OR20 (Figure 3a), compared to neutrally evolving sites. Notably, regions under positive selection in GR families overlapped with high-attention surface patches (Figure 3b) on sites consistent with structural markers such as S1–S4 helices (del Mármol et al. 2021), which are ligand-binding pockets, suggesting adaptive diversification of functionally important residues in GR. In contrast, OR20 exhibited fewer positively selected sites and a more diffuse distribution of functional patches, potentially reflecting weaker adaptive constraints or insufficient statistical
power due to the limited number of positively selected sites. This non-random association
suggests that diversifying selection preferentially targets functionally essential residues in GR
genes, supporting a mechanistic link between adaptive evolution and GR functional
diversification.

## Variable TE contents and evolutionary dynamics among aphid species

We predicted and annotated TEs as a preliminary step to examine their enrichment near OR and GR genes (see the next result part). The analysis revealed that TE prediction and annotation were not significantly correlated with genome assembly quality (see Tables S23 and S24 and Materials and Methods for details). We observed significant variation in TE content and evolutionary dynamics across aphid species. While no correlation between TE content and genome size was found (r²=0.02, NS, Table S3), differences in the proportion of TE types were noted, particularly Terminal Inverted Repeat (TIR) class II elements. Kimura distance-based divergence revealed varying evolutionary histories among species, with some species, such as *D. plantaginea* and *M. persicae*, showing recent TE activity, particularly in TIR elements. In contrast, others exhibited older and more stable TE dynamics. These differences suggest distinct evolutionary pressures on TE content across aphid species (Table S3, Figures 3 and S12).

## Younger transposable elements are more frequently enriched near OR genes than GR genes and other genomic regions

To investigate the role of TEs in shaping the diversity of GR and OR genes in aphids, we assessed TE enrichment in their vicinity. TEs were detected using the LOLA tool (Sheffield and Bock 2016), and their relationships to gene locations were evaluated using TEgrip (Meguerditchian et al., 2021). We observed a stronger enrichment of TEs near OR genes compared to GR genes when considering 2 kb and 10 kb flanking regions. However, this enrichment was not uniformly detected across all aphid species (Fisher's test, Bonferroni correction, q-values < 0.05; Figure 4, Tables S25). In *D. plantaginea*, 62.5% of GR genes and 88.6% of OR genes were associated with enriched TEs. Notably, TEs tended to encompass GR genes or lie proximally (<2 kb), while OR genes were more often flanked by nearby TEs without direct overlap. In *D. vitifoliae*, six distinct TEs were enriched around OR loci, with two elements
found near 31.7% and 55.5% of OR genes. Similarly, *R. padi* and *M. cerasi* exhibited enrichment patterns around OR genes (Table S26).

We further examined the evolutionary age of TEs associated with OR and GR genes compared to those in other genomic regions, using single-copy orthologous groups (OGs) as a reference. TE age was estimated via Kimura distance as a proxy for sequence divergence. TEs associated with OR genes were significantly younger than those linked to GR genes or other regions of the genome (Kruskal-Wallis test, *P* < 0.01; Figure 4), indicating more recent TE activity around OR loci. Despite this recent TE activity near OR genes, we mainly detected evidence of purifying selection across TE-associated OR and GR genes, except for GR03 (two TEs associated with four genes), GR04 (one TE linked to one gene), and GR06 (one TE linked to two genes) (Table S26). This suggests that recent TE insertions around OR genes have not driven widespread adaptive selection.

To assess whether TE proximity affects the predicted functional potential of chemosensory genes, we compared site-level attention scores derived from the MuLAN deep learning framework between TE-associated and non-TE-associated genes. Statistically significant differences were observed in both GRs (*P*-value = 3.7×10⁻⁶) and ORs (*P*-value < 2.2×10⁻¹⁶; Figure S15). Interestingly, GR proteins encoded by TE-associated GR genes tended to exhibit slightly higher functional scores, suggesting that TE activity may occasionally promote functional innovation. In contrast, TE-associated OR genes displayed lower attention scores, indicating that recent TE insertions might be more disruptive in OR loci, favoring pseudogenization or relaxed structural constraints rather than adaptive innovation.

Taken together, these findings reveal a mechanistic difference in how TE activity influences the evolution of OR and GR genes. OR loci appear evolutionarily younger—marked by more recent TE insertions, lower synteny, and star-like phylogenetic topologies—yet they show fewer adaptive signatures and reduced predicted functional potential. Conversely, older GR loci retain adaptive hotspots linked to functionally important encoded residues, where occasional TE insertions may have facilitated innovation rather than degradation.
# 1 Discussion

Our goal was to identify the repertoire and mode of evolution of chemosensory receptor genes that enable phytophagous insects to adapt to their host plants, as well as to explore the potential role of TEs in driving the evolutionary dynamics of these genes. We focused on aphids, a key family of phytophagous insects, for which a comprehensive analysis of gene family evolution has not yet been conducted beyond studies focused on crop pest species. We investigated gene family evolution, focusing specifically on OR and GR genes, which are essential for host recognition in phytophagous insects, as well as the diversity and evolutionary dynamics of TEs in aphids. This comprehensive study required extensive manual annotation, covering hundreds of OR and GR genes and TEs across the genomes of 12 aphid species, including several significant crop pests. Among these, we analyzed a newly assembled, high-quality genome of *D. plantaginea*, a key pest affecting fruit trees. Our findings indicate that gene families related to aphid-plant interactions, particularly OR and GR gene families, have undergone expansion through proximal duplications in an aphid subfamily characterized by its wide host range. Additionally, the evolution of these OR and GR genes is faster than the rest of the genome and influenced by diversifying selection, characterized by rapid episodes of positive selection followed by more extended periods of purifying selection, reflecting the adaptive dynamics involved in host recognition and specialization (Smadja et al. 2009; Engsontia et al. 2014; Gouin et al. 2017). To complement these evolutionary insights, we performed an *in silico* functional validation of positively selected sites using a deep learning-based predictor of structurally and functionally relevant residues, linking adaptive changes to potential mechanistic effects. The presence of younger TEs associated with OR genes, compared to GR genes and other genomic regions, along with the star-like structure of the OR phylogenetic tree and low synteny, suggests that OR genes have experienced more recent evolutionary changes. This aligns with observations in other model systems (Robertson 2019). The higher TE enrichment near OR and GR genes indicates their potential role in catalyzing tandem duplications, contributing to the diversification of these gene families. Additionally, we observed species-specific TE content and evolutionary dynamics, particularly in aphids infesting fruit trees, such as *D. plantaginea* and *M. persicae*. This study enhances our understanding of the evolutionary mechanisms that drive the adaptation
1 of phytophagous to hosts, as well as the role of TEs in shaping genomic architecture (Smadja et al. 2009; Engsontia et al. 2014; Gouin et al. 2017; Ye et al. 2022; Baril et al. 2023).

3

4 The impact of host range on chemosensory gene evolution in aphids

5 We found that the Aphidinae subfamily exhibited 3.8 times higher evolutionary rates compared to Eriosomatinae, driven by changes in gene families (including OR and GR gene families) associated with host interactions, including those related to lipid metabolism, immune responses, and transposase activity. The Aphidinae is the most species-rich subfamily of aphids and includes over 200 genera, characterized by various life cycles, including shifts between woody perennial primary hosts and multiple herbaceous secondary hosts (Blackman and Eastop 2000; Blackman and Eastop 2008). The higher evolutionary rates of duplicated gene families associated with plant interactions in this aphid subfamily, which exhibits a high host plant diversity, may suggest that host diversity can drive speciation in this subfamily (Peccoud et al. 2010). Previous studies have indicated a positive correlation between host range and speciation rates (Weingartner et al. 2006; Hardy and Otto 2014; Nicholson et al. 2015; Robertson et al. 2019). While OR and GR genes, as well as other gene families, exhibited duplications in some aphid species that could be linked to host diversity, further investigation is needed to determine how these expansions correlate with ecological adaptation. Complex factors, including seasonal and geographic variability, influence host breadth in aphids (Blackman and Eastop 2008). Some species are specialists, depending on a single host throughout their life cycle, whereas others are generalists, feeding on multiple hosts across different seasons (Dixon 1977). Therefore, while gene family expansions and contractions provide insights into aphid genomic evolution, additional ecological analyses are required to clarify their adaptive role.

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

Purifying selection is the dominant signature observed in OR and GR genes, affecting over 50% of these genes. However, these signatures are notably weaker than those in other genomic regions (OG), indicating that OR and GR genes in aphids evolve at higher rates. Additionally, approximately 20-30% of OR and GR genes show evidence of recent positive selection. These results, with the observed negative correlation between the evolutionary rates of GR and OR genes and species divergence, suggest that these genes undergo rapid episodes of positive selection, followed by sustained periods of purifying selection, consistent with a pattern
1 of diversifying selection. The observed co-occurrence of positive selection signatures and
2 elevated attention scores at specific sites within OR and GR genes suggests that these residues
3 are not only evolving adaptively but may also be functionally important. Because attention
4 scores from deep learning frameworks such as MuLAN are predictive of structurally or
5 functionally relevant positions in proteins, their enrichment in positively selected sites
6 strengthens the hypothesis that adaptive molecular changes with potential mechanistic relevance
7 drive chemosensory receptor diversification in aphids. Interestingly, not all families exhibiting
8 signatures of positive selection showed a corresponding enrichment in functionally relevant sites.
9 For example, OR20 presented limited or no significant differences in attention scores between
10 positively selected and neutral sites. This lack of signal may reflect either a low number of
11 selected sites in this family or structural regions under relaxed constraints, suggesting that not all
12 adaptive episodes necessarily lead to stable functional innovation. Such discrepancies underscore
13 the importance of combining evolutionary and structural predictions to distinguish between
14 transient diversification and functionally significant adaptation. This integrative evidence
15 supports the idea that duplicated OR and GR genes may undergo neo-functionalization or sub-
16 functionalization (Hallem et al. 2006; Assis and Bachtrog 2013; Meslin et al. 2022) via selection
17 acting on critical structural sites.

18 The interaction between host plants and herbivorous insects can generate the evolutionary
19 dynamics we inferred, with recurrent episodes of positive selection followed by long periods of
20 purifying selection that drive the diversification of chemosensory genes. In host–parasite
21 systems, such as those involving aphids and their host plants, variation in these genes is critical
22 for adapting to plant-produced defensive volatiles or tastant compounds (Simon et al. 2015; Züst
23 and Agrawal 2016; Beran and Petschenka 2022). Aphids with enhanced olfactory capabilities are
24 better equipped to detect and navigate toward suitable hosts from a distance, thereby increasing
25 their fitness. Once on the plant, improved gustatory perception enables aphids to identify optimal
26 feeding sites (Joseph and Carlson 2015). The faster evolutionary rate observed in OR genes
27 compared to GR genes likely reflects their primary role in long-distance host recognition,
28 whereas GRs function at a later stage, facilitating feeding after host location. However,
29 population genomic studies of selection footprints in candidate genes or genome-wide scans are
30 still needed to fully characterize these adaptive processes over shorter evolutionary timescales
1 (Simon et al. 2015; Shih et al. 2023) than our study. While such studies exist for *A. pisum*
2 (Smadja et al. 2009; Smadja et al. 2012; Eyres et al. 2016; Nouhaud et al. 2018), they must now
3 be extended to a broader range of aphid species and populations to reveal the generality of these
4 patterns.

5 Although this was not the primary objective of our study, our comprehensive annotation
6 of OR and GR genes across several aphid species provided an opportunity to compare their
7 repertoires among species. We did not find any variation in the size of the OR/GR gene
8 repertoires among aphid species. As we used a standardized pipeline on high-quality reference
9 genomes followed by manual curation, this lack of differences cannot be due to unidentified OR
10 and GR genes. Some studies concluded that the numbers of OR and GR genes are associated
11 with the host breadth of each aphid species (Quan et al. 2019; Robertson et al. 2019)(Quan et al.
12 2019; Robertson 2019). In other insect species, various conclusions are drawn from different
13 host specializations, such as the specialist fig wasp *Ceratosolen some* Mayr (Xiao et al. 2013),
14 which displays 46 OR and 6 GR genes, or the specialized human body louse *Pediculus humanus*
15 Haeckel, which has 10 ORs and 6 GRs (Kirkness et al. 2010). Despite the variable host
16 preferences of the studied aphid species, a higher sampling of more aphid species across various
17 subfamilies would help to better answer whether the increase in gene repertoires of OR and GR
18 genes occurred in association with host breadth. Even other chemosensory genes (e.g., CSPs and
19 OBPs, (Vogt 2003; Pelosi et al. 2006)) could be investigated, but additional manual annotation
20 would be required. Finally, specific genes may undergo pseudogenization following
21 specialization in the host plants in phytophagous conditions. This adaptation may result from a
22 reduced need for specific sensory capabilities when they specialize in particular plant types.
23 Here, we decided to focus on non-degraded sequences; these genomic fossils or evolutionary
24 relics may help better understand each genome's evolution rate, as they generally evolve without
25 selective constraints (Podlaha and Zhang 2010).

26
27 **Transposable elements and the evolution of GR and OR gene duplication in aphids**
28
29 TE content and organization in aphid genomes are crucial data for understanding the
30 evolutionary dynamics of duplicated genes that contribute to insect adaptation to their hosts. Our
31 research revealed that OR and GR genes in aphids likely originated from tandem duplications.
Similar tandem duplications have been observed in *Drosophila* species related to genes involved in defense against pathogens, insecticide resistance, chorion development, cuticular peptides, and lipases (Rogers et al. 2015).

High densities of transposons are frequently associated with duplicated genes in various insect species, indicating that TEs may facilitate the formation of proximal duplicated gene copies (Liu and Wessler 2017; Dupeyron et al. 2019; Gilbert et al. 2021). We found only a few enrichments of TEs near OR and GR genes. However, we found that TEs near OR and GR genes are younger than elsewhere in the genome, indicating heightened TE activity in these regions, particularly in species like *D. plantaginea*. The rapid turnover of TEs, characterized by high birth and death rates and strong purifying selection, may account for the limited detection of TE enrichment (Nei et al. 2000). TEs only persist as long as active elements can move to new positions within a genome. Over time, individual transposon lineages are often lost within specific clades. Transposons may have played a significant role in catalyzing the tandem duplication of OR and GR genes, thereby contributing to the diversity observed within these gene families. For example, in the ant *Cardiocondyla obscurior*, OR genes are located in TE-rich regions, suggesting a significant role for TEs in their duplication (McKenzie and Kronauer 2018). In the moth *Spodoptera spp.*, some TEs are significantly enriched near GR gene clusters (Meslin et al. 2022). TEs are also close to insecticide resistance genes in *H. armigera*. However, evidence of this effect appears to be limited to genes that have undergone recent diversification.

Interestingly, OR genes located near TEs showed significantly lower functional attention scores compared to ORs situated further away, while GR genes did not display this trend. At first glance, this seems counterintuitive because TE-proximal ORs also show signs of being evolutionarily younger. A plausible explanation is that these ORs are recently duplicated copies generated by TE-mediated events that have not yet stabilized into fully functional receptor structures. Such young duplicates may be evolving under relaxed or neutral selective pressures, allowing rapid sequence divergence without necessarily yielding structurally optimized or functionally mature proteins. Only a few GR families with TE associations (GR03, GR04, GR06) retained clear evidence of adaptive diversification. Moreover, TE insertions near promoters or within introns could modulate transcriptional activity, alternative splicing, or epigenetic silencing of OR/GR genes, ultimately influencing gene expression plasticity. Such effects have been documented in other insect models and may represent an underexplored layer
1 of functional diversification in aphids (Rebollo et al. 2012). The reduced attention scores
2 observed in TE-associated ORs might also reflect altered or disordered protein structures,
3 possibly resulting from sequence disruption or relaxed structural constraints following
4 duplication. Taken together, these findings suggest that TE activity has a dual role. In GRs, TE
5 insertions may facilitate adaptive diversification by introducing regulatory or structural novelty,
6 leading to functional innovation. In ORs, however, the same process may be more disruptive,
7 leading to relaxed constraints, reduced structural stability, and an increased likelihood of
8 pseudogenization. This explains why ORs appear evolutionarily younger (with more recent TE
9 insertions) but show fewer adaptive signatures and lower functional attention scores. In contrast,
10 GRs—despite being older—retain adaptive hotspots linked to functionally important residues.
11 Future work integrating protein structure prediction (e.g., AlphaFold, ProfileView),
12 transcriptomics, and ligand-binding assays will be needed to determine which OR proteins
13 encoded by TE-associated OR genes eventually gain new functions (neo-functionalization) and
14 which are destined to become nonfunctional.

15 More broadly, we also detected variations in the evolutionary dynamics of TEs among
16 different aphid species. However, it is premature to establish a direct connection between these
17 TE dynamics and the ecological or evolutionary histories of the species involved. Further
18 investigation of the transposon insertion landscape along the genome of aphid populations should
19 provide a deeper understanding of the selection and the demographic forces shaping the
20 adaptation of aphids to their host plant. Such a population genomics lens would provide insights
21 into the role of TEs in aphid evolution and, more generally, into the balance between
22 demography and selection in the evolution of insects (Lockton et al. 2008; Oggenfuss et al.
23 2021).

24

## Conclusion

26 This study reveals how duplication, selection, and TE dynamics shape gene evolution in insect
27 pests. Our findings indicate that OR and GR genes evolve faster than the rest of the genome
28 following gene duplication, shaped by diversifying selection with episodic bursts of positive
29 selection interspersed with periods of purifying selection. This pattern supports the hypothesis
30 that shifts to novel host plants in aphids drive adaptation to new environments, eventually
leading to species diversification. Additionally, our results suggest that TEs may play a dual role in shaping key genes associated with aphid adaptation—facilitating innovation in GRs while exerting more disruptive effects on ORs, which may lead to relaxed structural constraints or pseudogenization. Future research incorporating a broader range of aphid genomes with varying host breadths will help clarify the contributions of gene duplication and TE dynamics to host plant adaptation. Comparative analyses between highly polyphagous species and host-specialized aphids could provide further insights into the relationship between host range and the evolution of OR and GR genes. Additionally, investigating other gene families involved in aphid-host interactions, such as digestion genes, detoxification genes (Baril et al. 2023), other chemosensory genes (e.g., CSPs and OBPs (Vogt 2003; Pelosi et al. 2006)), and chemosensory ionotropic receptors (Rytz et al. 2013), alongside population genomic analyses and TE dynamics, will deepen our understanding of insect adaptation to their host. This study also delivers valuable genomic resources, including a high-quality genome assembly for *D. plantaginea* (the rosy apple aphid) and comprehensive annotations of ORs, GRs, and TEs. OR genes under recent positive selection present promising targets for environmentally sustainable pest control strategies based on GR and OR-functioning disruption. We anticipate that these genomic datasets will be essential for future investigations into the molecular mechanisms of host adaptation in aphids and other phytophagous insects. Importantly, by integrating evolutionary models with deep learning-based functional predictions, this study bridges adaptive molecular evolution with mechanistic relevance, providing a framework to prioritize candidate receptors for future functional assays.

22

## Material and methods

24

### Aphid genomic resources

We assembled a high-quality reference genome of *D. plantaginea*, the rosy apple aphid (see details in Text S1, Tables S2 and S3). We obtained the genomes of twelve aphid species and one aphid-like species from previous publications (Table S1): *Aphis glycines* (bt1; PRJNA594673), *Acyrthosiphon pisum* LSR1 (v2; PRJNA13657), *Aphis gossypii* (PRJNA934319), *Diuraphis noxia* (PRJNA226156), *Eriosoma lanigerum* (PRJNA362407), *M. persicae* (clone O;

30
1 PRJNA1106641), *Myzus cerasi* (PRJEB24287), *Pentalonia nigronervosa* (PRJNA573804),
2 *Rhopalosiphum maidis* (PRJNA1000643), *Rhopalosiphum padi* (PRJNA1037176), and *Sitobion
3 miscanthi* (PRJNA908645), and one aphid-like species of the Phylloxeridae family,
4 *Daktulosphaira vitifoliae* Fitch (PRJNA929591). We selected aphid species identified as
5 polyphagous or generalists with a broad host range to study gene evolution in a subfamily
6 exhibiting the greatest diversity of life history traits, primarily in terms of host breadth, when
7 applicable. The quality of the genomes was assessed using N50 (the shortest contig length
8 needed to cover 50% of the genome) and BUSCO (Benchmarking Universal Single-Copy
9 Orthologs) scores, as detailed in the publications (Table S1). We downloaded the genome
10 assemblies and annotations of the aphid species from the AphidBase (Bioinformatics platform
11 for Agroecosystem Arthropods called BIPAA (Legeai et al. 2010):
12 https://bipaa.genouest.org/is/aphidbase/).

13

14 Annotation of chemosensory receptor genes

15 Detailed annotation pipeline is presented in Text S2 and Figure S1. Briefly, we utilized
16 previously reported annotated OR and GR genes for A. *glycines*, *A. pisum*, and *D. vitifoliae*
17 (Robertson et al. 2019; Rispe et al. 2020) (Table S1) as input for Exonerate (Slater and Birney
18 2005), the InsectOR pipeline (Karpe et al. 2021), and Scipio (Keller et al. 2008) to predict OR
19 and GR genes across the 13 genomes mentioned above (Figure S1). We manually curated and
20 filtered the predicted gene models on the Apollo server (Lee et al. 2013) hosted on the BIPAA
21 platform using BLAST on NCBI to identify similarity percentages with previously reported
22 OR/GR gene sequences (Robertson et al. 2019) and confirmed the presence of conserved protein
23 domains (Sayers et al. 2019). After annotating the candidate OR/GR genes for each aphid
24 species, we initiated a second round of gene detection using the InsectOR pipeline (Karpe et al.
25 2021) to refine and improve the final set of OR and GR gene annotations. As a final validation,
26 we aligned the protein sequences for each gene family using MAFFT v7.4 (Kuraku et al. 2013;
27 Katoh et al. 2019) with the standard parameters and a gap-opening penalty of two. We removed
28 or split sequences with extremely short lengths or highly divergent sequences compared to
29 known protein sequences using available OR and GR sequences from *A. pisum* and *A. glycines*
30 reported by Robertson et al. as references (2019).

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30
We performed Spearman rank correlation tests between OR or GR repertoires and genome size estimated from the assemblies. We statistically tested for differences in numbers of OR and GR genes among aphid species using a generalized linear model (glm.nb function) as follows:

$$ \text{Gene count}_{ijk} = \mu + \text{species}_i + \text{families}_j + \varepsilon_{ijk} \quad (\text{Model 1}) $$

Where Gene countijk represented the number of genes of species i from the OR or GR family j, while μ denoted the overall mean. Species and families were treated as fixed effects, and εijk was the residual term. The residual was fitted to a Negative Binomial distribution using the R package MASS, version 7.3-58.3 (Venables and Ripley 2013).

We incorporated the annotated GR and OR genes into the previously retrieved genome annotation files (GFF) for the 12 aphid species and one aphid-like species. To eliminate overlapping regions, we utilized BCFTools version 1.5 (Danecek et al. 2021). These updated files served as input for detecting orthogroups and conducting the analyses described below.

## Ortholog identification and species tree inferences

We used the combined set of predicted proteins from the 13 species to infer orthogroups using OrthoFinder v2.4 (Emms and Kelly 2019). This analysis allowed us to identify conserved gene families and assess the evolutionary relationships among the species, facilitating a comprehensive understanding of gene repertoire evolution in the aphid lineages studied (Emms and Kelly 2019). We then constructed a phylogenetic tree of the 13 species using single-copy orthologous groups identified by OrthoFinder v2.4 (Emms and Kelly 2019). This approach enabled us to accurately depict the evolutionary relationships among the aphid species accurately and provided a foundation for further analyses of gene evolution. We aligned and concatenated the single-copy orthologs using MAFFT (Kuraku et al. 2013; Katoh et al. 2019), ensuring a comprehensive alignment for subsequent phylogenetic analyses. We used TrimAl (Capella-Gutiérrez et al. 2009) to remove the poorly aligned regions with a conservation threshold of 60% and a gap threshold of 90%. We constructed a maximum likelihood phylogenetic tree using IQ-TREE v2.0.6 (Minh et al. 2020) based on the aligned single-copy orthologous sequences and
conducted 10,000 bootstrap replicates. The phylogenetic tree was edited with FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/), with *D. vitifoliae* serving as the outgroup.

## Expansion and contraction of gene families

We utilized the Computational Analysis of Gene Family Evolution (CAFE v5; Mendes et al. 2021)) to identify rapidly evolving gene families, which indicate a high rate of genomic turnover (gains and losses per gene per million years). AFE estimates the birth-death parameter $\lambda$ (probability that any gene will be gained or lost at each node and terminal branch, accounting for gene family expansions and contractions) along a provided phylogenetic tree scaled to time with the same length of root tips and gene family count estimated with OrthoFinder (Emms and Kelly 2019).

First, we constructed an ultrametric species tree from the IQ-TREE analysis using the r8s software (http://sourceforge.net/projects/r8s/). We used divergence times and a reference species tree created with the TimeTree website (www.timetree.org) as inputs, along with five available time references (Kim et al. 2011; Ren et al. 2013; Hardy et al. 2015).

We ran the CAFE software using the ultrametric tree and gene family counts inferred from OrthoFinder. To ensure accuracy, we employed a Python script to eliminate gene families with gene copies in only one species or that contained more than 100 gene copies in one or more species. This correction in gene counting helps to avoid potential errors in the estimation of evolutionary rates ((Mendes et al. 2021); https://github.com/hahnlab/cafe_tutorial/). We conducted the CAFE analysis using several modes to estimate the net gain and loss rates. These included 1) a single parameter ($\lambda$) across the entire phylogeny and 2) multiple parameters for each gene family. We established a single evolutionary rate (gene birth-death rate) as our null hypothesis and proposed different evolutionary rates along each branch of the phylogeny as our alternative hypothesis. The evolutionary hypotheses tested included factors such as host preference (tree versus herbaceous host), specialization (generalist vs. specialist), phylogeny, and taxonomic rank (family, subfamily, tribe, and genus). Finally, we evaluated the fit of the models by comparing the log-likelihoods of the alternative and null models, utilizing a likelihood-ratio test to identify significant differences (LRT; Magis et al. 2010; Thissen et al. 2013)) script on RStudio v1.3.1.
We performed a Gene Ontology (GO) enrichment analysis on the 147 expanded and 51 significantly contracted gene groups identified by CAFE. We employed the TopGO R package (Alexa and Rahnenfuhrer 2010) with the “elim” method to identify genes associated with biological processes or functions based on Gene Ontology (GO) terms.

## Micro-synteny analyses and mode of evolution of GR and OR genes

Expanded gene families often occur as tandem arrays, a genomic architecture that can contribute to increased gene birth and death rates, increasing copy number variation among species (Ohno 2013). Therefore, we examined how genomic organization varies between GR and OR subfamilies in aphid species to generate insights into the molecular evolution shaping OR and GR subfamily diversity. We built a synteny network using the R/Bioconductor package Syntenet (Almeida-Silva et al. 2023). We used chemosensory protein sequences and the genomic coordinates for each species (gff files). We explored the synteny of OR and GR genes, focusing only on gene clusters containing the OR and GR genes previously annotated (Table S1). We visualized the phylogenomic profiles of the syntenic genes as a heat map. We performed the protein alignment using DIAMOND (Buchfink et al. 2015) in sensitive mode (e-value = 1e-10, top hits = 5). To identify and classify duplicated gene pairs, we used the doubletrouble R package (Almeida-Silva and Van de Peer 2025), applying its default parameters for detecting whole genome duplication (WGD), small-scale duplications (SSD), and paralogous genes based on the syntenic blocks.

## Selective pressures acting on OR and GR genes

We inferred selection pressures within OR and GR gene subfamilies defined from the phylogenetic trees we built. OR and GR subfamilies were monophyletic groups supported by high bootstrap values (>80%). We also visually inspected the tree topology. We then defined the foreground and background branches for the selection analysis. We detailed the methods below.

First, we built two phylogenetic trees for GR and OR genes. Fructose and Orco genes were used as outgroups (Robertson 2019) to root the trees. From the protein files of each GR and OR genes, we extracted their coding sequences (CDS) from the original FASTA files of the respective genome assemblies. We aligned the sequences using MUSCLE (Edgar 2004) with default parameters and a maximum of 16 iterations. We utilized Gblocks (Castresana 2000) to
curate the sequences, restricting contiguous non-conserved positions and removing any lines that
introduced significant gaps in the alignment. We also removed sequences <150 pb or presenting
stop codons (Figure S7 and Table S11). The sequences were then converted to codon alignments
using PAL2NAL (Suyama et al. 2006). We constructed separate phylogenetic trees for the OR
and GR genes using IQ-TREE v2.0 (Minh et al. 2020), employing the model finder option,
ultrafast bootstrap with 1,000 replicates, and a branch support test (SH-aLRT) with 1,000
replicates. The OR and GR subfamilies were visually identified on the phylogenetic trees as
monophyletic groups with bootstrap support exceeding 80%.

We used the maximum likelihood approach implemented in codon-based analysis (codeml) with PAML version 4.2 (Yang 2007) to assess positive and purifying selection for each GR and OR subfamily. Codeml computes the likelihood of a large set of molecular evolution models that estimate $\omega$, the ratio of the non-synonymous mutation rate ($dN$) over the synonymous mutation rate ($dS$) (Kimura 1980). The most likely model of molecular evolution explaining the data is chosen based on a likelihood ratio test (LTR) statistic $2\hat{\lambda}L$ (with a $\chi^2$ distribution with degrees of freedom equal to the differences in the number of parameters) computed among models. Once the most likely model was selected, evidence of positive selection was inferred from the estimates of $\omega$: values of $\omega > 1$ indicated positive selection, while values of $\omega < 1$ indicated purifying selection.

We utilized the branch (free-ratio) and the site models. We did not use the branch-site model as we neither had any phylogenetic background nor an apparent OR or GR gene family clustering to test for selection. The branch model assumes different $\omega$ ratios for different branches on the phylogeny (Nielsen and Yang 1998; Yang 1998) and was used to detect positive selection acting on particular sequences of the subfamilies without averaging $\omega$ throughout the phylogenetic tree. Using a Likelihood Ratio Test, the branch model was compared to a null hypothesis ($\omega=1$, one-ratio M0; (Nielsen and Yang 1998)).

We also tested the site model, which can identify positively selected sites in a multiple-sequence alignment (Yang and Nielsen 2002). The models can be site-class-specific, all of which assume that the $\omega$ ratio is the same across branches of the phylogeny but different among sites in the alignment. These codon substitution models are M0 (one-ratio), M1a (nearly neutral), M2a (positive selection), M3 (discrete), M7 (beta), M8 (beta and $\omega > 1$), and M8a (beta and $\omega = 1$). We applied the LTR statistics to evaluate the fit of the nested models to the data for positive
selection (Nielsen and Yang 1998) confronting M0 vs. M3, M1a vs. M2a, and M7 vs. M8 models. Support for positive selection can be identified if M2a provides a better fit than M1a or M8 provides a better fit than M7 or M8a (Yang et al. 2000; Anisimova et al. 2001). We confirmed signals of positive selection when a better fit of the selected model was supporting M2a and M8.

The distribution of $dN/dS$ ratio estimated from branch models for each gene among OR, GR, and Orthologous Groups (OG) families and species was analyzed with the following the general linear model (function *glm.nb* in R package *lme4*) as follows:

$$dN/dS_{ijkl} = \mu + \text{species}_i + \text{subfamily}_j + \text{family}_k + \varepsilon_{ijkl} \quad (\text{Model 2}).$$

Where $dN/dS_{ijkl}$ is the $dN/dS$ ratio of the sequence from subfamily j in the OR/GR/OG family k of species i, $\mu$ represents the overall mean, species and families are treated as fixed effects, and $\varepsilon_{ijkl}$ is the residual. The residual was fitted to a Negative Binomial using the R package Mass, version 7.3-58.3 (Venables and Ripley 2013). We compared the fit of three models to the data using the *anova* function (aov R package, (Chambers et al. 2017)). The statistical significance of an effect was assessed with an ANOVA-type test and Tukey's Honest Significance Test (R function Tukey-HSD) to compare $dN/dS$ among gene families (OG, GR, OR) pairwise. Predicted $dN/dS$ values from the best were then plotted with the *plot_model* function (sjPlot R package, (Lüdecke 2013)).

## Repetitive DNA element annotation

Examining transposable element (TE) content and evolutionary patterns around OR and GR genes in various aphid species can provide insights into the role of TEs driving OR and GR gene diversification. However, variability in TE annotation methods across aphid genomes poses a challenge for direct comparisons. Standardizing annotation tools is crucial for understanding TE-mediated evolution in these genes and clarifying the genetic mechanisms underlying host adaptation in aphids.

We used the REPET package version 3.0 (https://urgi.versailles.inra.fr/Tools/REPET) (Quesneville et al. 2005; Flutre et al. 2011) for the de novo identification and annotation of TEs in the 13 genomes under survey. REPET includes two main pipelines: Tedenovo (Flutre et al.
2011), which detects repeated elements in a genome to build a TE reference library, and TEannot (Quesneville et al. 2005), which annotates TEs based on a TE consensus library. Consensus sequences were built from at least three repeated sequences identified in the genome using the TEdenovo pipeline (Flutre et al. 2011). The GROUPER program (Quesneville et al. 2005) was used for the TEdenovo High Scoring Pairs (HSP) clustering step. These consensus sequences were then classified using the Wicker code (Wicker et al. 2007) and PASTEC (Hoede et al. 2014). We used the TEannot pipeline (Quesneville et al. 2005) to annotate all repetitive elements in the genome, utilizing the TE consensus library from TEdenovo. TE annotation was done using similarity searches with BLASTER (Quesneville et al. 2003), RepeatMasker (Smit et al. 2015), and CENSOR (Jurka et al. 1996). Short simple repeats (SSRs) were annotated using TRF (Benson 1999) and RepeatMasker (Smit et al. 2015) and removed from the TE annotations. Fragments of the same copy were joined using MATCHER (Quesneville et al. 2003; Quesneville et al. 2005) and the “long join procedure” (Ahmed et al. 2011). To improve the quality of the TE reference library from TEdenovo, we applied a “second TEannot process” (Jamilloux et al. 2017). The first round of TEannot identifies consensuses that annotate at least one full-length copy (FLC), with both fragmented and unfragmented annotations aligned over more than 95% of the consensus TE sequence. The second round of TEannot is then performed using only those consensuses that annotated at least one FLC.

We test whether the TE content order differed among aphid species using a general linear model (glm) with a quasipoisson distribution and ANOVA test with R/Rstudio (http://cran.r-project.org/mirrors.html) as follows:

$$TE \ count_{ijk} = Species_i + TE \ order_j + \varepsilon_{ijk} \ (Model \ 3).$$

The TE count represented the number of TE copies, species indicated the species analyzed, and TEorder referred to the TE orders corresponding to each TE copy. We also retrieved the sequences of the TE copies identified by REPET to prepare the consensus. We compared the identity of the TE copies among species using an upset plot (UpSetR R package).
## 1 TE family divergence

We compared TE content and age along the 12 aphid and the one aphid-like genomes. We calculated the Kimura distances (sequence divergence) using scripts of the RepeatMasker package (Flynn et al. 2020). We used the RepeatMasker annotation as input to generate the landscape of aphid species TEs. Kimura's two-parameter (K2P) model is a nucleotide substitution model for estimating genetic distances and phylogenetic relationships, considering transitional and transversional substitution rates (Kimura 1980). We plot the proportion of base pair repeats (Y-axis) at different levels of divergence (X-axis) calculated as a 2-p Kimura distance using the ggplot2 R package (Wickham 2016). We also utilized the ggplot2 R package to create a plot describing the aphid TE content distribution using the Kimura substitution level as a function of the percentage of the genome. We defined the TE content according to the Wicker TE classification (Wicker et al. 2007) utilized in the PASTEC classifier pipeline from the REPET package (Quesneville et al. 2005).

## 14 TE enrichment near OR and GR genes

We used the locus overlap analysis within the LOLA R package Bioconductor (Sheffield and Bock 2016) to test for TE enrichment near OR/GR genes in the 12 aphid and the one aphid-like genomes. We built three different datasets to launch the analysis. The query set contained genomic regions of 10 Kb around each OR/GR gene. The region universe (or the control genomic region) comprised gene coordinates of all genic areas of the genome and expanded 10 Kb around each gene. We tested for TE enrichment within OR/GR regions compared to the control regions using the Fisher's Exact Test with false discovery rate correction (q-value<0.1).

We used the TEGRIP algorithm (Meguerditchian et al. 2021) to characterize relationships between enriched TEs and the chemosensory genes to which they are associated: "subset" means that the TE is included inside the gene, "superset" that the gene is included inside the TE, "upstream" indicate a TE located upstream of the gene, otherwise "downstream". Finally, "overlap" denotes at least one base of the gene overlapping the downstream or upstream part of the TE. As TEgrip only returns the closest TE to the gene, we used a modified script to conserve all TEs (https://forgemia.inra.fr/sergio.olvera-vazquez/aphid_chemosensory genes/-/blob/main/Genomic_architecture_TE/relationships_analysis.R). We also calculated Kimura distances between each copy of enriched TE and the consensus sequences using ape R package (Paradis and Schliep 2019) (https://forgemia.inra.fr/sergio.olvera-
vazquez/aphid_chemosensory genes/-/blob/main/Kimura_distance/Kimura_distances.R). For TE
detected enriched and without classification following the REPET pipeline, we used DeepTE to
improve its classification (Yan et al. 2020). Statistical differences between the distribution of
Kimura distances of TE copies associated with GR genes and those associated with OR genes
were computed using Wilcoxon signed-rank test (Wilcoxon 1945).

Identification of functionally relevant residues using MuLAN

To assess the potential functional significance of sites within OR and GR receptor genes, which show signatures of positive selection in all three codeml comparisons — M1a vs. M2a, M7 vs. M8, and M0 vs. M3, with all being significant —, the deep learning framework MuLAN (Mutation-driven Light Attention Networks) (Lombardi and Carbone 2024) was used. Although initially developed to predict the impact of amino acid mutations on protein-protein interactions, MuLAN can also be applied to individual protein sequences to estimate residue-level attention scores. These scores, normalized between 0 and 1, are inferred from local sequence context and have been shown to correlate with binding or functionally important regions. Full-length coding sequences of each OR and GR gene family were first translated into amino acid sequences. Sequences within each gene family were aligned using MAFFT v7.505 with the E-INS-i strategy, which is suitable for proteins containing multiple conserved domains. Per-residue attention scores were then computed by running MuLAN on each sequence independently. These scores were subsequently mapped onto the multiple sequence alignment and averaged across aligned positions within each gene family. The resulting attention profiles were plotted alongside codon sites previously inferred to be under positive selection, as identified by codeml analyses (see “Codon-based estimation of selection rates” section). Comparative analyses were performed to evaluate whether attention scores differed significantly between positively selected and non-selected sites, using Wilcoxon-Mann-Whitney U tests (Wilcoxon 1945). These plots included (i) average attention score profiles per gene family, (ii) overlays of codon sites inferred under positive selection (posterior probability > 0.95 under M8 or M2a models), and (iii) boxplots of attention scores grouped by selection regime (positive, purifying, neutral).
To evaluate whether TE proximity influences the functional potential of OR and GR genes, attention scores derived from the MuLAN deep learning model were compared between genes located near TEs and those located further away. Distributions of attention scores for TE-associated and non-TE-associated genes were compared separately for OR and GR gene families using the Wilcoxon rank sum test with continuity correction. Statistical significance thresholds were set at p < 0.05, and results were visualized using boxplots stratified by gene family and TE association status. This analysis was designed to identify whether local TE context correlates with reduced or increased residue-level functional signal, potentially reflecting structural disruption.

All MuLAN analyses were conducted using the published implementation, and the alignment and plotting pipeline is available in the project's GitHub repository (https://gitlab.lcqb.upmc.fr/lombardi/mulan).

## Data availability

The newly assembled genome of *D. plantaginea* is publicly available on the BIPAA platform (https://bipaa.genouest.org/sp/dysaphis_plantaginea/). The TE consensus sequence and annotations are available on RepetDB (https://urgi.versailles.inrae.fr/repetdb); see detailed URLs in Table S27. The OR/GR CDS sequence and gene structure annotations files are available on the ZENODO repository: DOI 10.5281/zenodo.12750441. The attention scores computed for OR and GR protein sequences are available on the ZENODO repository: https://doi.org/10.5281/zenodo.16257464. All analysis files and codes are available there: https://github.com/CornilleEclecticLab/aphid-comparative-genomics.

## Acknowledgments

We thank the Bioinformatics Platform for Agroecosystem Arthropods (BIPAA) platform, particularly Fabrice Legeai (Institut de Génétique, Environnement et Protection des Plantes, INRAE, Institut Agro, Université Rennes) and the Centre National de Ressources Génomiques Végétales-CNRGV, especially William Marande, Nathalie Rodde, and Stéphane Cauet (INRAE), for the *de novo* assembly of *D. plantaginea*. We thank the Plant Bioinformatics
1 Facility, BioinfOmics, INRAE, Université Paris-Saclay. We thank Hugh David Loxdale for the
discussion on specialization in aphids.

3

4 Author contributions.

5 AC obtained funding, and AC, EJJ, and CM designed the experiment. CB, DO, AD, ES, BB, FL,
6 and AC prepared the material for genome assemblies; SGOV, AC, AM, GL, XC, FAS, JC, CL,
7 IE, CM, EJJ, and YB analyzed the data; All co-authors discussed the results. SGOV, AM, and
8 AC wrote the manuscript with critical input from other co-authors.

9

10 Funding

11 This research was funded by the ATIP-CNRS Inserm, IDEEV, LabEx BASC, the PARSADA
12 program (ARDECO platform, FranceAgriMer: PAEDP0924005531), and Tamkeen, under the
13 grant AD454 from the New York University Abu Dhabi Research Institute.

14

15 Conflict of interest

16 The authors of this preprint declare that they have no financial conflict of interest with the
17 content of this article.

18
**Figure 1.** Gene family expansion and contraction events in 12 aphid species and the outgroup *Daktulosphaira vitifoliae* detected using CAFE. The internal branch with the most significant number of gene gains (expansions) and losses (contractions) corresponds to the most recent common ancestor of aphid species from the *Aphidinae* subfamily. Species abbreviations are as follows: AGLY = *Aphis glycines* clone BT1; AGOS = *Aphis gossypii*; APIS = *Acyrthosiphon pisum* clone LSR1; DNOX = *Diuraphis noxia*; DPLA = *Dysaphis plantaginea*; DVIT = *Daktulosphaira vitifoliae* (outgroup); ELAN = *Eriosoma lanigerum*; MCER = *Myzus cerasi*; MPER = *Myzus persicae* clone O; PNIG = *Pentalonia nigronervosa*; RMAI = *Rhopalosiphum maidis*; RPAD = *Rhopalosiphum padi*; SMIS = *Sitobion miscanthi*. Gene gains (+) are shown in green, and gene losses (-) in purple. Significant gains and losses of olfactory receptor (OR) and gustatory receptor (GR) genes are highlighted in orange and blue, respectively. All branches are supported by 100% bootstrap values. ● indicates the estimated evolutionary rate (birth-and-death rate) for each branch.

**Figure 2.** Mode of evolution of olfactory and gustatory receptor genes (OR and GR, respectively) distribution in the 12 aphid species and the grape phylloxera *D. vitifoliae*. *a.* Percentage of duplication modes of each aphid species inferred with R package doubletrouble. *b.* Dispersed duplication: genes cannot be classified as segmental, tandem, or proximal.

**Figure 3.** Functional relevance of positively selected sites in aphid gustatory and olfactory receptors. (A) Site-level attention scores computed with the MuLAN deep learning framework for gustatory receptor families GR03, GR04, GR06, and olfactory receptor family OR20. Sites under positive selection (detected by codeml models M1a vs. M2a, M7 vs. M8, and M0 vs. M3) show significantly higher attention scores than neutrally evolving sites in GR03, GR04, and GR06 (Wilcoxon-Mann-Whitney test, *p* < 2.2 × 10¹⁶). In contrast, no significant difference was detected for OR20. The dashed horizontal line indicates the median attention score across all sites. (B) Predicted 3D structures with the Boltz-1 model (Wohlwend et al. 2025) of GR03, GR04, GR06, and OR20 mapped with MuLAN attention scores. Residues are color-coded from blue (low attention score) to red (high attention score), highlighting predicted functional and binding sites on ligand-binding pockets, consistent with structural markers such as S1-S4 helices. Positive selection sites are represented as spheres in the predicted structures. For each gene family (GR03, GR04, GR06, OR20), a single representative protein sequence was chosen
for 3D structure visualization because members within the same family showed highly similar predicted structures.

**Figure 4. Transposable element (TE) content in the twelve aphid species and the grape phylloxera in a phylogenetic context (IQtree, with all bootstraps > 0.99).** A) Proportion of TEs in the genome for each analyzed species. The purple portion represents the fraction of the genome occupied by TEs. This proportion is based on the respective TE library for each species. B) Proportion of Class I (yellow to red) and II (green to blue) consensus TEs in each analyzed species. C) Kimura distance distribution for the most recently active TE across the analyzed species (0% divergence). The x-axis represents the Kimura distance for each TE category. "Order" is used instead of "category" for precision. The TE categories (e.g., DNA_MITE, LINE, SINE, etc.) are indicated on the x-axis. Species abbreviations: AGLY = *Aphis glycines* clone BT1; AGOS = *Aphis gossypii*; APIS = *Acyrthosiphon pisum* clone LSR1; DNOX = *Diuraphis noxia*; DPLA = *Dysaphis plantaginea*; DVIT = *Daktulosphaira vitifoliae*; ELAN = *Eriosoma lanigerum*; MCER = *Myzus cerasi*; MPER = *Myzus persicae* clone O; PNIG = *Pentalonia nigronervosa*; RMAI = *Rhopalosiphum maidis*; RPAD = *Rhopalosiphum padi*; SMIS = *Sitobion miscanthi*.

**Figure 5. Distribution and age of transposable element (TE) enrichment near olfactory and gustatory receptor genes in 12 aphid species.** (A) TE Enrichment analysis. This panel illustrates the enrichment of individual TE copies near olfactory receptor (OR) and gustatory receptor (GR) genes based on 15 annotated consensus TE sequences. The analysis identifies significant enrichment within 10 kb of OR and GR genes, determined by mapping TE copies to their genomic locations. Since consensus TE sequences represent generalized forms, they do not directly indicate enrichment but serve as a reference for identifying individual TE insertions. TE-gene relationships are categorized as follows: Subset (light orange): TE located within the gene. Superset (light yellow): Gene located within the TE. Downstream (red): TE located downstream of the gene. Upstream (light green): TE located upstream of the gene. Overlap (orange): Gene overlaps with TE upstream/downstream regions by at least one base. Among the 12 aphid species analyzed, only five species—*Aphis glycines* (AGLY), *Dysaphis plantaginea* (DPLA), *Daktulosphaira vitifoliae* (DVIT), *Myzus cerasi* (MCER), and *Rhopalosiphum padi* (RPAD)—exhibited significant TE enrichment near OR and GR genes. (B) TE Age Distribution (Kimura
Distance Analysis) The distribution of Kimura distances is shown for TE copies associated with OR and GR genes compared to all TE copies in the genome. The Kimura distance measures sequence divergence, estimating TE insertion age. TE copies are categorized as follows: All TE copies (black), TE copies associated with GR genes (blue), TE copies associated with OR genes (yellow). Statistical differences between these distributions were assessed using Kruskal-Wallis tests, with significance levels indicated as ns: Not significant (p > 0.05), p < 0.05, p < 0.01, p < 0.005, p < 0.001. This analysis provides insights into whether TE insertions near OR and GR genes are evolutionarily recent or ancient across aphid species.

9
10
11
1 References

2 Abdullaev ET, Umarova IR, Arndt PF. 2021. Modelling segmental duplications in the human

3 genome. BMC Genomics 22:496.

4 Ahmad A, Wallau GL, Ren Z. 2021. Characterization of Mariner transposons in seven species of

5 Rhus gall aphids. Sci Rep 11:16349.

6 Ahmed I, Sarazin A, Bowler C, Colot V, Quesneville H. 2011. Genome-wide evidence for local

7 DNA methylation spreading from small RNA-targeted sequences in Arabidopsis. Nucleic

8 Acids Research 39:6919-6931.

9 Alexa A, Rahnenfuhrer J. 2010. topGO: Enrichment analysis for Gene Ontology. R package

10 version 2.18.0.

11 Almeida FC, Sánchez-Gracia A, Campos JL, Rozas J. 2014. Family Size Evolution in

12 Drosophila Chemosensory Gene Families: A Comparative Analysis with a Critical

13 Appraisal of Methods. Genome Biology and Evolution 6:1669-1682.

14 Almeida-Silva F, Van de Peer Y. 2025. doubletrouble: an R/Bioconductor package for the

15 identification, classification, and analysis of gene and genome duplications.

16 Bioinformatics 41:btaf043.

17 Almeida-Silva F, Zhao T, Ullrich KK, Schranz ME, Van De Peer Y. 2023. syntenet: an

18 R/Bioconductor package for the inference and analysis of synteny networks. Martelli PL,

19 editor. Bioinformatics 39:btac806.

20 Anisimova M, Bielawski JP, Yang Z. 2001. Accuracy and Power of the Likelihood Ratio Test in

21 Detecting Adaptive Molecular Evolution. Molecular Biology and Evolution 18:1585-

22 1592.

23 Arnab SP, Amin MR, DeGiorgio M. 2023. Uncovering Footprints of Natural Selection Through

24 Spectral Analysis of Genomic Summary Statistics. Kim Y, editor. Molecular Biology and

25 Evolution 40:msad157.

26 Assis R, Bachtrog D. 2013. Neofunctionalization of young duplicate genes in Drosophila.

27 Proceedings of the National Academy of Sciences 110:17409-17414.

28 Baril T, Pym A, Bass C, Hayward A. 2023. Transposon accumulation at xenobiotic gene family

29 loci in aphids. Genome Research 33:1718-1733.

30 Barton NH, Mallet J. 1996. Natural Selection and Random Genetic Drift as Causes of Evolution

31 on Islands [and Discussion]. Philos Trans R Soc Lond B Biol Sci 351:785.

32 Benson G. 1999. Tandem repeats finder: a program to analyze DNA sequences. Nucleic Acids

33 Research 27:573-580.

34 Beran F, Petschenka G. 2022. Sequestration of Plant Defense Compounds by Insects: From

35 Mechanisms to Insect-Plant Coevolution. Annual Review of Entomology 67:163-180.
1

2

3

4

Biello R, Singh A, Godfrey CJ, Fernández FF, Mugford ST, Powell G, Hogenhout SA, Mathers TC. 2021. A chromosome-level genome assembly of the woolly apple aphid, *Eriosoma lanigerum* Hausmann (Hemiptera: Aphididae). *Molecular Ecology Resources* 21:316–326.

5

6

Blackman RL, Eastop VF. 2000. Aphids on the world's crops: an identification and information guide (No. Ed. 2). John Wiley & Sons Ltd.

7

8

Blackman RL, Eastop VF. 2008. Aphids on the world's herbaceous plants and shrubs, 2 volume. set. John Wiley & Sons

9

10

11

Bouallègue M, Filée J, Kharrat I, Mezghani-Khemakhem M, Rouault J-D, Makni M, Capy P. 2017. Diversity and evolution of mariner-like elements in aphid genomes. *BMC Genomics* 18:494.

12

13

Buchfink B, Xie C, Huson DH. 2015. Fast and sensitive protein alignment using DIAMOND. *Nat Methods* 12:59–60.

14

15

16

17

Byrne S, Schughart M, Carolan JC, Gaffney M, Thorpe P, Malloch G, Wilkinson T, McNamara L. 2022. Genome sequence of the English grain aphid, *Sitobion avenae* and its endosymbiont *Buchnera aphidicola*. Scofield S, editor. *G3 Genes|Genomes|Genetics* 12:jkab418.

18

19

20

Cacciò SM, Chalmers RM, Dorny P, Robertson LJ. 2018. Foodborne parasites: Outbreaks and outbreak investigations. A meeting report from the European network for foodborne parasites (Euro-FBP). *Food and Waterborne Parasitology* 10:1–5.

21

22

Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. 2009. *trimAl*: a tool for automated alignment trimming in large-scale phylogenetic analyses. *Bioinformatics* 25:1972–1973.

23

24

Castresana J. 2000. Selection of Conserved Blocks from Multiple Alignments for Their Use in Phylogenetic Analysis. *Molecular Biology and Evolution* 17:540–552.

25

26

27

Chambers JM, Freeny AE, Heiberger RM. 2017. Analysis of variance; designed experiments. In: *In Statistical models in S*. Routledge. p. 145–193. Available from: https://doi.org/10.1201/9780203738535

28

29

30

Choi H, Shin S, Jung S, Clarke DJ, Lee S. 2018. Molecular phylogeny of Macrosiphini (Hemiptera: Aphididae): An evolutionary hypothesis for the Pterocomma-group habitat adaptation. *Molecular Phylogenetics and Evolution* 121:12–22.

31

32

33

Cornille A, Ebert D, Stukenbrock E, Rodríguez de la Vega RC, Tiffin P, Croll D, Tellier A. 2022. Unraveling coevolutionary dynamics using ecological genomics. *Trends in Genetics* Jun 14.

34

35

36

Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM, et al. 2021. Twelve years of SAMtools and BCFtools. *GigaScience* 10:giab008.
1 Desvignes T, Lauridsen H, Valdivieso A, Fontenele RS, Kraberger S, Murray KN, Le François
2 NR, Detrich HW, Kent ML, Varsani A, et al. 2022. A parasite outbreak in notothenioid
3 fish in an Antarctic fjord. *iScience* 25:104588.

4 Dixon AFG. 1977. Aphid Ecology: Life Cycles, Polymorphism, and Population Regulation.
5 *Annual Review of Ecology and Systematics*:329-353.

6 Dupeyron M, Singh KS, Bass C, Hayward A. 2019. Evolution of Mutator transposable elements
7 across eukaryotic diversity. *Mobile DNA* 10:12.

8 Duplessis S, Cuomo CA, Lin Y-C, Aerts A, Tisserant E, Veneault-Fourrey C, Joly DL, Hacquard
9 S, Amselem J, Cantarel BL, et al. 2011. Obligate biotrophy features unraveled by the
10 genomic analysis of rust fungi. *Proc. Natl. Acad. Sci. U.S.A.* 108:9166-9171.

11 Duran C, Edwards D, Batley J. 2009. Genetic Maps and the Use of Synteny. In: Gustafson JP, Langridge P, Somers DJ, editors. Plant Genomics. Vol. 513. Methods in Molecular Biology™. Totowa, NJ: Humana Press. p. 41–55. Available from:
14 http://link.springer.com/10.1007/978-1-59745-427-8_3

15 Ebert D, Fields PD. 2020. Host–parasite co-evolution and its genomic signature. *Nature Reviews Genetics* 21:754–768.

17 Edgar RC. 2004. MUSCLE: multiple sequence alignment with high accuracy and high
18 throughput. *Nucleic Acids Research* 32:1792–1797.

19 Elliott TA, Gregory TR. 2015. Do larger genomes contain more diverse transposable elements?
20 *BMC Evol Biol* 15:69.

21 Emmns DM, Kelly S. 2019. OrthoFinder: phylogenetic orthology inference for comparative
22 genomics. *Genome Biol* 20:238.

23 Engsontia P, Sangket U, Chotigeat W, Satasook C. 2014. Molecular Evolution of the Odorant
24 and Gustatory Receptor Genes in Lepidopteran Insects: Implications for Their Adaptation
25 and Speciation. *Journal of Molecular Evolution* 79:21–39.

26 Eyres I, Jaquiéry J, Sugio A, Duvaux L, Gharbi K, Zhou J, Legeai F, Nelson M, Simon J, Smadja
27 CM, et al. 2016. Differential gene expression according to race and host plant in the pea
28 aphid. *Molecular Ecology* 25:4197–4215.

29 Eyun S, Soh HY, Posavi M, Munro JB, Hughes DST, Murali SC, Qu J, Dugan S, Lee SL, Chao
30 H, et al. 2017. Evolutionary History of Chemosensory-Related Gene Families across the
31 Arthropoda. *Molecular Biology and Evolution* 34:1838–1862.

32 Fay JC, Wittkopp PJ. 2008. Evaluating the role of natural selection in the evolution of gene
33 regulation. *Heredity* 100:191–199.

34 Flutre T, Duprat E, Feuillet C, Quesneville H. 2011. Considering Transposable Element
35 Diversification in De Novo Annotation Approaches. Xu Y, editor. *PLoS ONE* 6:e16526.
1 Flynn JM, Hubley R, Goubert C, Rosen J, Clark AG, Feschotte C, Smit AF. 2020.

2 RepeatModeler2 for automated genomic discovery of transposable element families.

3 *Proc. Natl. Acad. Sci. U.S.A.* 117:9451–9457.

4 Frantzeskakis L, Di Pietro A, Rep M, Schirawski J, Wu C, Panstruga R. 2020. Rapid evolution in plant–microbe interactions – a molecular genomics perspective. *New Phytologist* 225:1134–1142.

7 Frantzeskakis L, Kusch S, Panstruga R. 2019. The need for speed: compartmentalized genome evolution in filamentous phytopathogens. *Molecular Plant Pathology* 20:3–7.

9 Gardiner A, Barker D, Butlin RK, Jordan WC, Ritchie MG. 2008. *Drosophila chemoreceptor gene evolution: selection, specialization and genome size*. *Molecular Ecology* 17:1648–1657.

12 Gilbert C, Peccoud J, Cordaux R. 2021. Transposable Elements and the Evolution of Insects.

13 *Annu. Rev. Entomol.* 66:355–372.

14 Gladieux P, Ropars J, Badouin H, Branca A, Aguileta G, De Vienne DM, Rodríguez De La Vega RC, Branco S, Giraud T. 2014. Fungal evolutionary genomics provides insight into the mechanisms of adaptive divergence in eukaryotes. *Molecular Ecology* 23:753–773.

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

21 Gawler CD, Essington H, O’Brien B, Shaw CL, Bilich RW, Clay PA, Duffy MA. 2023. Virulence evolution during a naturally occurring parasite outbreak. *Evol Ecol* 37:113–129.

24 Haley OC, Tibbs-Cortes L, Hayford RK, Harding S, Woodhouse M, Cannon E, Gardiner J, Portwood J, Sen TZ, Kim H-S, et al. 2025. Why do some predicted protein structures fold poorly? Benchmarking AlphaFold, ESMFold, and Boltz in maize. *bioRxiv* 2025.07.05.663230.

28 Hallem EA, Dahanukar A, Carlson JR. 2006. Insect odor and taste receptors. *Annu. Rev. Entomol.* 51:113–135.

30 Hansson BS, Stensmyr MC. 2011. Evolution of Insect Olfaction. *Neuron* 72:698–711.

31 Hardy NB, Otto SP. 2014. Specialization and generalization in the diversification of phytophagous insects: tests of the musical chairs and oscillation hypotheses. *Proc. R. Soc. B.* 281:20132960.

34 Hardy NB, Peterson DA, Von Dohlen CD. 2015. The evolution of life cycle complexity in aphids: Ecological optimization or historical constraint?: ADAPTIVENESS OF COMPLEX LIFE CYCLE. *Evolution* 69:1423–1432.
1 Hartmann FE, Sánchez-Vallet A, McDonald BA, Croll D. 2017. A fungal wheat pathogen evolved host specialization by extensive chromosomal rearrangements. *The ISME Journal* 11:1189–1204.

4 He Z, Zhang Z, Guo W, Zhang Y, Zhou R, Shi S. 2015. De Novo Assembly of Coding Sequences of the Mangrove Palm (Nypa fruticans) Using RNA-Seq and Discovery of Whole-Genome Duplications in the Ancestor of Palms. *PLoS ONE* 10:e0145385.

7 Hoede C, Arnoux S, Moisset M, Chaumier T, Inizan O, Jamilloux V, Quesneville H. 2014. PASTEC: An Automatic Transposable Element Classification Tool. *Cordaux R*, editor. *PLoS ONE* 9:e91929.

10 Huang T, Liu Y, He K, Francis F, Wang B, Wang G. 2023. Chromosome-level genome assembly of the spotted alfalfa aphid Therioaphis trifolii. *Sci Data* 10:274.

12 Jacquin-Joly E, Merlin C. 2004. Insect olfactory receptors: contributions of molecular biology to chemical ecology. *Journal of Chemical Ecology* 30:2359–2397.

14 Jamilloux, Daron, Choulet, H. Quesneville. 2017. De Novo Annotation of Transposable Elements: Tackling the Fat Genome Issue. *Proceedings of the IEEE* 105:474–481.

16 Joseph RM, Carlson JR. 2015. Drosophila Chemoreceptors: A Molecular Interface Between the Chemical World and the Brain. *Trends in Genetics* 31:683–695.

18 Julca I, Marcet-Houben M, Cruz F, Vargas-Chavez C, Johnston JS, Gómez-Garrido J, Frias L, Corvelo A, Loska D, Cámara F, et al. 2020. Phylogenomics Identifies an Ancestral Burst of Gene Duplications Predating the Diversification of Aphidomorpha. *Battistuzzi FU*, editor. *Molecular Biology and Evolution* 37:730–756.

22 Jurka J, Klonowski P, Dagman V, Pelton P. 1996. Censor—a program for identification and elimination of repetitive elements from DNA sequences. *Computers & Chemistry* 20:119–121.

25 Karlsson EK, Kwiatkowski DP, Sabeti PC. 2014. Natural selection and infectious disease in human populations. *Nat Rev Genet* 15:379–393.

27 Karpe SD, Tiwari V, Ramanathan S. 2021. InsectOR—Webserver for sensitive identification of insect olfactory receptor genes from non-model genomes. *PLOS ONE* 16:e0245324.

29 Katoh K, Rozewicki J, Yamada KD. 2019. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. *Briefings in Bioinformatics* 20:1160–1166.

31 Keller O, Odronitz F, Stanke M, Kollmar M, Waack S. 2008. Scipio: Using protein sequences to determine the precise exon/intron structures of genes and their orthologs in closely related species. *BMC Bioinformatics* 9:278.
1 Kim H, Lee S, Jang Y. 2011. Macroevolutionary Patterns in the Aphidini Aphids (Hemiptera: Aphididae): Diversification, Host Association, and Biogeographic Origins. Moreau CS, editor. *PLoS ONE* 6:e24749.

4 Kimura M. 1977. Preponderance of synonymous changes as evidence for the neutral theory of molecular evolution. *Nature* 267:275-276.

6 Kimura M. 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. *J Mol Evol* 16:111-120.

8 Kirk H, Dorn S, Mazzi D. 2013. Molecular genetics and genomics generate new insights into invertebrate pest invasions. *Evolutionary Applications* 6:842-856.

10 Kirkness EF, Haas BJ, Sun W, Braig HR, Perotti MA, Clark JM, Lee SH, Robertson HM, Kennedy RC, Elhaik E, et al. 2010. Genome sequences of the human body louse and its primary endosymbiont provide insights into the permanent parasitic lifestyle. *Proc. Natl. Acad. Sci. U.S.A.* 107:12168-12173.

14 Klai K, Chénais B, Zidi M, Djebbi S, Caruso A, DENIS F, CONFAIS J, BADAWI M, CASSE N, MEZGHANI KHEMAKHEM M. 2020. Screening of Helicoverpa armigera Mobilome Revealed Transposable Element Insertions in Insecticide Resistance Genes. *Insects* 11.

17 Kondrashov FA. 2012. Gene duplication as a mechanism of genomic adaptation to a changing environment. *Proceedings of the Royal Society B: Biological Sciences* 279:5048-5057.

19 Kondrashov FA, Rogozin IB, Wolf YI, Koonin EV. 2002. Selection in the evolution of gene duplications. *Genome Biology* 3:research0008.1.

21 Kulmuni J, Wurm Y, Pamilo P. 2013. Comparative genomics of chemosensory protein genes reveals rapid evolution and positive selection in ant-specific duplicates. *Heredity* 110:538-547.

24 Kuraku S, Zmasek CM, Nishimura O, Katoh K. 2013. aLeaves facilitates on-demand exploration of metazoan gene family trees on MAFFT sequence alignment server with enhanced interactivity. *Nucleic Acids Research* 41:W22-W28.

27 Lee E, Helt GA, Reese JT, Munoz-Torres MC, Childers CP, Buels RM, Stein L, Holmes IH, Elsik CG, Lewis SE. 2013. Web Apollo: a web-based genomic annotation editing platform. *Genome Biol* 14:R93.

30 Legeai F, Shigenobu S, Gauthier J-P., Colbourne J, Rispe C, Collin O, Richards S, Wilson ACC, Murphy T, Tagu D. 2010. AphidBase: a centralized bioinformatic resource for annotation of the pea aphid genome. *Insect Molecular Biology* 19:5-12.

33 Lehman SS, Verhoeve Victoria I., Driscoll Timothy P., Beckmann John F., Gillespie Joseph J. 2024. Metagenome diversity illuminates the origins of pathogen effectors. *mBio* 15:e00759-23.
1 Li Z, Capoduro R, Bastin-Héline L, Zhang S, Sun D, Lucas P, Dabir-Moghaddam D, François
2 M-C, Liu Y, Wang G, et al. 2023. A tale of two copies: Evolutionary trajectories of moth
3 pheromone receptors. *Proceedings of the National Academy of Sciences*
4 120:e2221166120.

5 Liu K, Wessler SR. 2017. Functional characterization of the active Mutator-like transposable
6 element, Mutat1 from the mosquito Aedes aegypti. *Mobile DNA* 8:1.

7 Lockton S, Ross-Ibarra J, Gaut BS. 2008. Demography and weak selection drive patterns of
8 transposable element diversity in natural populations of *Arabidopsis lyrata*. *Proc. Natl.
9 Acad. Sci. U.S.A.* 105:13965-13970.

10 Lombardi G, Carbone A. 2024. MuLAN: Mutation-driven Light Attention Networks for
11 investigating protein-protein interactions from sequences. *bioRxiv* :2024.08.24.609515.

12 Lüdecke D. 2013. sjPlot: Data Visualization for Statistics in Social Science. :2.8.16. Available
13 from: https://CRAN.R-project.org/package=sjPlot

14 Lynch M. 2023. Mutation pressure, drift, and the pace of molecular coevolution. *Proc. Natl.
15 Acad. Sci. U.S.A.* 120:e2306741120.

16 Magadum S, Banerjee U, Murugan P, Gangapur D, Ravikesavan R. 2013. Gene duplication as a
17 major force in evolution. *J Genet* 92:155-161.

18 Magis D, Béland S, Tuerlinckx F, De Boeck P. 2010. A general framework and an R package for
19 the detection of dichotomous differential item functioning. *Behavior Research Methods*
20 42:847-862.

21 Märkle H, John S, Cornille A, Fields PD, Tellier A. 2021. Novel genomic approaches to study
22 antagonistic coevolution between hosts and parasites. *Mol Ecol* 30:3660-3676.

23 del Mármol J, Yedlin MA, Ruta V. 2021. The structural basis of odorant recognition in insect
24 olfactory receptors. *Nature* 597:126-131.

25 Mathers TC, Mugford ST, Hogenhout SA, Tripathi L. 2020. Genome Sequence of the Banana
26 Aphid, Pentalonia nigronervosa Coquerel (Hemiptera: Aphididae) and Its Symbionts. *G3
27 Genes* | *Genomes* | *Genetics* 10:4315-4321.

28 Mathers TC, Wouters RHM, Mugford ST, Swarbreck D, Van Oosterhout C, Hogenhout SA.
29 2021. Chromosome-Scale Genome Assemblies of Aphids Reveal Extensively Rearranged
30 Autosomes and Long-Term Conservation of the X Chromosome. *Ouangraoua A*, editor.
31 *Molecular Biology and Evolution* 38:856-875.

32 McBride CS. 2007. Rapid evolution of smell and taste receptor genes during host specialization
33 in *Drosophila sechellia*. *Proc. Natl. Acad. Sci. U.S.A.* 104:4996-5001.

34 McCullers TJ, Steiniger M. 2017. Transposable elements in *Drosophila*. *Mobile Genetic Elements* 7:1-18.
1 McKenzie SK, Kronauer DJC. 2018. The genomic architecture and molecular evolution of ant odorant receptors. *Genome Res.* 28:1757-1765.

3 Meguerditchian C, Ergun A, Decroocq V, Lefebvre M, Bui Q-T. 2021. Pipeline to detect the positional and directional relationship between transposable elements and adjacent genes in host genome. *Peer Community Journal* [Internet] 1. Available from: https://peercommunityjournal.org/articles/10.24072/pcjournal.57/

7 Mendes FK, Vanderpool D, Fulton B, Hahn MW. 2021. CAFE 5 models variation in evolutionary rates among gene families. Robinson P, editor. *Bioinformatics* 36:5516-5518.

10 Meslin C, Mainet P, Montagné N, Robin S, Legeai F, Bretaudieu A, Johnston JS, Koutroumpa F, Persyn E, Monsempès C, et al. 2022. Spodoptera littoralis genome mining brings insights on the dynamic of expansion of gustatory receptors in polyphagous noctuidae. *Sethuraman A*, editor. *G3 Genes* | *Genomes* | *Genetics* 12:jkac131.

14 Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, Von Haeseler A, Lanfear R. 2020. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era. *Teeling E*, editor. *Molecular Biology and Evolution* 37:1530-1534.

17 Moller M, Stukenbrock EH. 2017. Evolution and genome architecture in fungal plant pathogens. *Nature Reviews Microbiology* advance online publication:756-771.

19 Nei M, Rogozin IB, Piontkivska H. 2000. Purifying selection and birth-and-death evolution in the ubiquitin gene family. *Proceedings of the National Academy of Sciences* 97:10866-10871.

22 Nicholson SJ, Nickerson ML, Dean M, Song Y, Hoyt PR, Rhee H, Kim C, Puterka GJ. 2015. The genome of Diuraphis noxia, a global aphid pest of small grains. *BMC Genomics* 16:429.

25 Nielsen R, Yang Z. 1998. Likelihood Models for Detecting Positively Selected Amino Acid Sites and Applications to the HIV-1 Envelope Gene. *Genetics* 148:929-936.

27 Niimura Y, Nei M. 2007. Extensive Gains and Losses of Olfactory Receptor Genes in Mammalian Evolution. *Hahn M*, editor. *PLoS ONE* 2:e708.

29 Nouhaud P, Gautier M, Gouin A, Jaquiéry J, Peccoud J, Legeai F, Mieuzet L, Smadja C, Lemaitre C, Vitalis R, et al. 2018. Identifying Genomic Hotspots Of Differentiation And Candidate Genes Involved In The Adaptive Divergence Of Pea Aphid Host Races. *Molecular Ecology* 27.

33 O'Donnell S, Rezende G, Vernadet J-P, Snirc A, Ropars J. 2025. Harboring Starships: The Accumulation of Large Horizontal Gene Transfers in Domesticated and Pathogenic Fungi. *Genome Biology and Evolution* 17:evaf125.
1 Oggenfuss U, Badet T, Wicker T, Hartmann FE, Singh NK, Abraham L, Karisto P, Vonlanthen
2 T, Mundt C, McDonald BA, et al. 2021. A population-level invasion by transposable
3 elements triggers genome expansion in a fungal pathogen. *eLife* 10:e69249.

4 Oggenfuss U, Feurtey A, Reyes-Avila CS, Gluck-Thaler E, Puccetti G, Glad HM, Abraham LN,
5 Stalder L, Tralamazza SM, González-Sáyer SM, et al. 2023. Genome Evolution in Fungal
6 Plant Pathogens: From Populations to Kingdom-Wide Dynamics. In: Pöggeler S, James
7 T, editors. Evolution of Fungi and Fungal-Like Organisms. Cham: Springer International
8 Publishing. p. 103–121. Available from: https://doi.org/10.1007/978-3-031-29199-9_5

9 Ohno S. 2013. Evolution by gene duplication. Springer Science & Business Media Available
10 from: https://doi.org/10.1007/978-3-642-86659-3

11 Paradis E, Schliep K. 2019. ape 5.0: an environment for modern phylogenetics and evolutionary
12 analyses in R. *Bioinformatics* 35:526–528.

13 Pearce SL, Clarke DF, East PD, Elfekih S, Gordon KHJ, Jermiin LS, McGaughran A, Oakeshott
14 JG, Papanikolaou A, Perera OP, et al. 2017. Genomic innovations, transcriptional
15 plasticity and gene loss underlying the evolution and divergence of two highly
16 polyphagous and invasive Helicoverpa pest species. *BMC Biology* 15:63.

17 Peccoud J, Simon J-C, von Dohlen C, Coeur d'acier A, Plantegenest M, Vanlerberghe-Masutti F,
18 Jousselin E. 2010. Evolutionary history of aphid-plant associations and their role in aphid
19 diversification. *Comptes Rendus Biologies* 333:474–487.

20 Pelosi P, Zhou J-J, Ban LP, Calvello M. 2006. Soluble proteins in insect chemical
21 communication. *Cellular and Molecular Life Sciences CMLS* 63:1658–1676.

22 Plissonneau C, Hartmann FE, Croll D. 2018. Pangenome analyses of the wheat pathogen
23 Zymoseptoria tritici reveal the structural basis of a highly plastic eukaryotic genome.
24 *BMC Biol* 16:5.

25 Podlaha O, Zhang J. 2010. Pseudogenes and Their Evolution. In: Encyclopedia of Life Sciences.
26 1st ed. Wiley. Available from:
27 https://online library.wiley.com/doi/10.1002/9780470015902.a0005118.pub2

28 Qiao X, Li Q, Yin H, Qi K, Li L, Wang R, Zhang S, Paterson AH. 2019. Gene duplication and
29 evolution in recurring polyploidization–diploidization cycles in plants. *Genome Biol*
30 20:38.

31 Quan Q, Hu X, Pan B, Zeng B, Wu N, Fang G, Cao Y, Chen X, Li X, Huang Y, et al. 2019.
32 Draft genome of the cotton aphid Aphis gossypii. *Insect Biochemistry and Molecular
33 Biology* 105:25–32.

34 Quesneville H, Bergman CM, Andrieu O, Autard D, Nouaud D, Ashburner M, Anxolabehere D.
35 2005. Combined Evidence Annotation of Transposable Elements in Genome Sequences.
36 *PLoS Comp Biol* 1:e22.
1 Quesneville H, Nouaud D, Anxolabéhere D. 2003. Detection of New Transposable Element

2 Families in Drosophila melanogaster and Anopheles gambiae Genomes. *Journal of*

3 *Molecular Evolution* 57:S50–S59.

4 Rebollo R, Romanish MT, Mager DL. 2012. Transposable Elements: An Abundant and Natural

5 Source of Regulatory Sequences for Host Genes. *Annual Review of Genetics* 46:21–42.

6 Remick BC, Gaidt MM, Vance RE. 2023. Effector-Triggered Immunity. *Annual Review of*

7 *Immunology* 41:453–481.

8 Ren Z, Zhong Y, Kurosu U, Aoki S, Ma E, Von Dohlen CD, Wen J. 2013. Historical

9 biogeography of Eastern Asian–Eastern North American disjunct Melaphidina aphids

10 (Hemiptera: Aphididae: Eriosomatinae) on Rhus hosts (Anacardiaceae). *Molecular*

11 *Phylogenetics and Evolution* 69:1146–1158.

12 Rispe C, Legeai F, Nabity PD, Fernández R, Arora AK, Baa-Puyoulet P, Banfill CR, Bao L,

13 Barberà M, Bouallègue M, et al. 2020. The genome sequence of the grape phylloxera

14 provides insights into the evolution, adaptation, and invasion routes of an iconic pest.

15 *BMC Biol* 18:90.

16 Robertson HM. 2015. The Insect Chemoreceptor Superfamily Is Ancient in Animals. *Chemical*

17 *Senses* 40:609–614.

18 Robertson HM. 2019. Molecular Evolution of the Major Arthropod Chemoreceptor Gene

19 Families. *Annual Review of Entomology* 64:227–242.

20 Robertson HM, Robertson ECN, Walden KKO, Enders LS, Miller NJ. 2019. The

21 chemoreceptors and odorant binding proteins of the soybean and pea aphids. *Insect*

22 *Biochemistry and Molecular Biology* 105:69–78.

23 Rocafort M, Bowen JK, Hassing B, Cox MP, McGreal B, De La Rosa S, Plummer KM,

24 Bradshaw RE, Mesarich CH. 2022. The Venturia inaequalis effector repertoire is

25 dominated by expanded families with predicted structural similarity, but unrelated

26 sequence, to avirulence proteins from other plant-pathogenic fungi. *BMC Biol* 20:246.

27 Rogers RL, Cridland JM, Shao L, Hu TT, Andolfatto P, Thornton KR. 2015. Tandem

28 Duplications and the Limits of Natural Selection in Drosophila yakuba and Drosophila

29 simulans.Palsson A, editor. *PLoS ONE* 10:e0132184.

30 Rytz R, Croset V, Benton R. 2013. Ionotropic Receptors (IRs): Chemosensory ionotropic

31 glutamate receptors in Drosophila and beyond. *Insect Biochemistry and Molecular*

32 *Biology* 43:888–897.

33 Sayers EW, Agarwala R, Bolton EE, Brister JR, Canese K, Clark K, Connor R, Fiorini N, Funk

34 K, Hefferon T, et al. 2019. Database resources of the National Center for Biotechnology

35 Information. *Nucleic Acids Research* 47:D23–D28.
1 Schrader L, Schmitz J. 2019. The impact of transposable elements in adaptive evolution.

2 *Molecular Ecology* 28:1537–1549.

3 Seidl MF, Thomma BPHJ. 2017. Transposable Elements Direct The Coevolution between Plants

4 and Microbes. *Trends in Genetics* 33:842–851.

5 Seong K, Krasileva KV. 2023. Prediction of effector protein structures from fungal

6 phytopathogens enables evolutionary analyses. *Nat Microbiol* 8:174–187.

7 Sheffield NC, Bock C. 2016. LOLA: enrichment analysis for genomic region sets and regulatory

8 elements in R and Bioconductor. *Bioinformatics* 32:587–589.

9 Shih P-Y, Sugio A, Simon J-C. 2023. Molecular Mechanisms Underlying Host Plant Specificity

10 in Aphids. *Annu. Rev. Entomol.* 68:431–450.

11 Simon J-C, d'Alencon E, Guy E, Jacquin-Joly E, Jaquiyer J, Nouhaud P, Peccoud J, Sugio A,

12 Streiff R. 2015. Genomics of adaptation to host-plants in herbivorous insects. *Briefings in*

13 *Functional Genomics* 14:413–423.

14 Slater GSC, Birney E. 2005. Automated generation of heuristics for biological sequence

15 comparison. *BMC Bioinformatics* 6:31.

16 Smadja C, Shi P, Butlin RK, Robertson HM. 2009. Large gene family expansions and adaptive

17 evolution for odorant and gustatory receptors in the pea aphid, *Acyrthosiphon pisum*. *Mol*

18 *Biol Evol* [Internet] 26. Available from: http://dx.doi.org/10.1093/molbev/msp116

19 Smadja CM, Canbäck B, Vitalis R, Gautier M, Ferrari J, Zhou JJ. 2012. Large-scale candidate

20 gene scan reveals the role of chemoreceptor genes in host plant specialization and

21 speciation in the pea aphid. *Evolution* [Internet] 66. Available from:

22 http://dx.doi.org/10.1111/j.1558-5646.2012.01612.x

23 Smit A, Hubley R., Green P. 2015. RepeatMasker Open-4.0. 2013–2015. Available from:

24 http://www.repeatmasker.org

25 Stergiopoulos I, De Wit PJGM. 2009. Fungal Effector Proteins. *Annu. Rev. Phytopathol.* 47:233–

26 263.

27 Suyama M, Torrents D, Bork P. 2006. PAL2NAL: robust conversion of protein sequence

28 alignments into the corresponding codon alignments. *Nucleic Acids Research* 34:W609–

29 W612.

30 The International Aphid Genomics Consortium. 2010. Genome Sequence of the Pea Aphid

31 *Acyrthosiphon pisum*. Eisen JA, editor. *PLoS Biol* 8:e1000313.

32 Thissen D, Steinberg L, Wainer H. 2013. Use of item response theory in the study of group

33 differences in trace lines. In: In Test validity. Routledge. p. 147–169.
1 Thorpe P, Escudero-Martinez CM, Cock PJA, Eves-van Den Akker S, Bos JIB. 2018. Shared
2 Transcriptional Control and Disparate Gain and Loss of Aphid Parasitism Genes. Moran
3 N, editor. *Genome Biology and Evolution* 10:2716–2733.

4 Torres DE, Oggenfuss U, Croll D, Seidl MF. 2020. Genome evolution in fungal plant pathogens: looking beyond the two-speed genome model. *Fungal Biology Reviews* 34:136–143.

6 Venables WN, Ripley BD. 2013. Modern applied statistics with S-PLUS. Springer Science & Business Media

8 Vogt RG. 2003. 14 - Biochemical diversity of odor detection: OBPs, ODEs and SNMPs. In: Blomquist G, Vogt R, editors. Insect Pheromone Biochemistry and Molecular Biology. San Diego: Academic Press. p. 391–445. Available from: https://www.sciencedirect.com/science/article/pii/B9780121071516500165

12 Vosshall LB, Stocker RF. 2007. Molecular Architecture of Smell and Taste in Drosophila.
13 *Annual Review of Neuroscience* 30:505–533.

14 Wagner A. 1998. The fate of duplicated genes: loss or new function? *Bioessays* 20:785-788.

15 Weingartner E, Wahlberg N, Nylin S. 2006. Dynamics of host plant use and species diversity in *Polygonia* butterflies (Nymphalidae). *Journal of Evolutionary Biology* 19:483–491.

17 White PS, Arslan D, Kim D, Penley M, Morran L. 2021. Host genetic drift and adaptation in the evolution and maintenance of parasite resistance. *J of Evolutionary Biology* 34:845–851.

19 Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, Flavell A, Leroy P, Morgante M, Panaud O, et al. 2007. A unified classification system for eukaryotic transposable elements. *Nat Rev Genet* 8:973–982.

22 Wickham H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York
23 Available from: https://ggplot2.tidyverse.org

24 Wilcoxon F. 1945. Individual Comparisons by Ranking Methods. *Biometrics Bulletin* 1:80-83.

25 Wohlwend J, Corso G, Passaro S, Getz N, Reveiz M, Leidal K, Swiderski W, Atkinson L, Portnoi T, Chinn I, et al. 2025. Boltz-1 Democratizing Biomolecular Interaction Modeling. *bioRxiv*:2024.11.19.624167.

28 Xiao J-H, Yue Z, Jia L-Y, Yang X-H, Niu L-H, Wang Z, Zhang P, Sun B-F, He S-M, Li Z, et al. 2013. Obligate mutualism within a host drives the extreme specialization of a fig wasp genome. *Genome Biol* 14:R141–R141.

31 Xu W, Papanicolaou A, Zhang H-J, Anderson A. 2016. Expansion of a bitter taste receptor family in a polyphagous insect herbivore. *Scientific Reports* 6:23666.

33 Yan H, Bombarely A, Li S. 2020. DeepTE: a computational method for de novo classification of transposons with convolutional neural network. *Bioinformatics* 36:4269–4275.
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

Yang Z. 1998. On the Best Evolutionary Rate for Phylogenetic Analysis. *Cannatella* D, editor.

*Systematic Biology* 47:125–133.

Yang Z. 2007. *PAML 4: Phylogenetic Analysis by Maximum Likelihood*. *Molecular Biology*

*and Evolution* 24:1586–1591.

Yang Z, Nielsen R. 2002. Codon-Substitution Models for Detecting Molecular Adaptation at

Individual Sites Along Specific Lineages. *Molecular Biology and Evolution* 19:908–917.

Yang Z, Nielsen R, Goldman N, Pedersen A-MK. 2000. Codon-Substitution Models for

Heterogeneous Selection Pressure at Amino Acid Sites. *Genetics* 155:431–449.

Ye X, Yang Y, Zhao C, Xiao S, Sun YH, He C, Xiong S, Zhao X, Zhang B, Lin H, et al. 2022.

Genomic signatures associated with maintenance of genome stability and venom turnover

in two parasitoid wasps. *Nature Communications* 13:6417.

Zhang S, Gao X, Wang L, Jiang, Su H, Jing T, Yang y. 2022. Chromosome-level genome

assemblies of two cotton-melon aphid *Aphis gossypii* biotypes unveil mechanisms of

host adaption. *Molecular Ecology Resources* 22:1120–1134.

Züst T, Agrawal AA. 2016. Mechanisms and evolution of plant resistance to aphids. *Nature*

*Plants* 2:15206.
<figure><img src="image_3.png" /><figcaption>Figure 1
159x120 mm (x DPI)</figcaption></figure>

1

2

3

4
<figure><img src="image_4.png" /><figcaption>Figure 2
159x89 mm (x DPI)</figcaption></figure>

1

2

3

4
<figure><img src="image_5.png" /><figcaption>Figure 3
150x160 mm (x DPI)</figcaption></figure>

1
2
3
4
<figure><img src="image_6.png" /><figcaption>Figure 4
159x90 mm (x DPI)</figcaption></figure>

1
2
3
4
A)

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

B)

<figure><img src="image_8.png" /><figcaption>Figure 5
159x95 mm (x DPI)</figcaption></figure>

1
2
3