Genomic basis of adaptation to constant and fluctuating environments in a global pest of cereals

genome and its annotation are available via the ORCAE platform 28 (https://bioinformatics.psb.ugent.be/orcae/overview/Aceto). Scripts used to analy se the 29 evolve-and-resequence experiment are available at 30 https://github.com/konczal/WCM_EvolveResequence 31 32 Funding statement: This study was supported by the National Science Centre (NSC), Poland, 33 research grants no. 2017/27/N/NZ8/00305 and 2019/32/T/NZ8/00151, which were awarded 34 to Alicja Laska-Modzelewska. The experiment measuring frequency changes of candidate SNPs 35 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 2 in smooth brome was supported by the National Science Centre (NSC), Poland, research grant 36 no. 2021/41/B/NZ8/01703 awarded to Anna Skoracka. 37 38 Conflict of interest disclosure: The authors declare no conflict of interest. 39 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 3

40 In agricultural systems, spatially and temporally heterogeneous environments are expected to 41 favour the evolution of generalist pests and pathogens, yet the genomic mechanisms 42 underlying niche-breadth expansion remain poorly understood. Here, we address this gap by 43 assembling one of the smallest known eukaryotic genomes—from the global cereal pest Aceria 44 tosichella—and by resequencing generalist and specialist populations which evolved in the lab 45 on constant versus fluctuating hosts. A previous study demonstrated that generalist 46 populations retained a wider niche, including the ability to exploit a refuge host, compared to 47 specialist populations that showed a narrowing of the ecological niche. Genomic scans 48 identified 640 SNPs that significantly differed in frequency between treatments, including a 49 single highly differentiated 120-kbp region containing 13 genes. Allele-frequency shifts across 50 many SNPs in this region paralleled phenotypic divergence in the ability to utilise refuge hosts, 51 and gene annotations suggested their roles in starvation resistance and nutrient signalling. 52 However, experimental validation of the top candidate gene did not support its strong direct 53 effect on survival on the refuge host, implying a more complex, likely polygenic basis for 54 adaptation. Consistent with this interpretation, numerous significant SNPs occurred outside 55 the focal region, with enrichment for genes involved in xenobiotic transport, detoxification 56 pathways, and ABC transporters. Together, these results indicate that generalism in A. 57 tosichella relies on complex molecular mechanisms that enable survival on suboptimal or 58 refuge hosts, relying to a considerable extent on the ability to deal with stress arising, for 59 example, from toxic compounds present in suboptimal hosts. 60 61 Key words: cereal pest, evolve-and-resequence, environmental heterogeneity, host specificity, 62 niche breadth, wheat curl mite 63 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 4 1. Introduction 64 Generalists are frequently expected to cope better with environmental change, yet this 65 adaptive advantage may come at the cost of reduced efficiency in exploiting specific resources 66 compared to specialists (Egas et al. 2004; Sexton et al. 2017). This trade -off leads to the 67 prediction that spatially and temporally heterogeneous environments favour the evolution of 68 generalists, whereas stable environments promote speciali sation (Egas et al. 2004; Levins 69 1968). However, such trade-offs appear to be far from universal (Bono et al. 2017, 2020; 70 Gompert et al. 2015; Remold 2012). While some studies support this prediction (Condon et al. 71 2014; Kassen 2002; Sant et al. 2021), others do not (Ketola et al. 2013; Saarinen et al. 2018). 72 Therefore, to understand why some adaptations involve trade-offs while others do not, it is 73 necessary to investigate the genetic mechanisms underlying the evolution of ecological niche 74 breadth (VanWallendael et al. 2019). 75 The genetic basis of such differences is only just beginning to be revealed. Recent studies 76 suggest that generalists may rely on generic response s, such as the upregulation of major 77 stress response pathways—for example, heat shock proteins (Leonard and Lancaster 2022; 78 Olazcuaga et al. 2023) —or increased transcriptional plasticity (Birnbaum and Abbot 2020). 79 Elucidating these mechanisms is particularly important for understanding how pests adapt to 80 croplands, which constitute the foundation of human food security. Modern agriculture relies 81 heavily on monocultures , which provide immense resource abundance but support low 82 biodiversity. Monocultures typically favour specialist pests, whereas crop rotation is thought 83 to limit their spread. Indeed, specialist herbivores tend to be more prevalent under 84 monotonous conditions and concentrated resources than in polycultures (Altieri 1999; Andow 85 1991). Conversely, crop rotation or pulsed resources (e.g. due to harvesting) are expected to 86 favour generalists capable of exploiting diverse resources (Kennedy and Storer 2000). On the 87 one hand, this may affect their population growth and subsequent damage by imposing the 88 costs of being a generalist; on the other, it allows them to invade new hosts of agricultural or 89 ecological value (Litovska et al. 2026; Paredes et al. 2021; Poveda et al. 2025) . Therefore, 90 understanding how agricultural practices shape the evolution of pest niche breadth —and 91 thereby their capacity to exploit crops and persist in agroecosystems —is of major scientific 92 and practical importance. 93 We recently demonstrated that the one of the most devastating global cereal pests, the wheat 94 curl mite (Aceria tosichella , Keifer ), can rapidly adapt to different levels of environmental 95 variability. Bread wheat (Triticum aestivum), which is typically grown in monocultures, is the 96 primary host plant of the wheat curl mite. After harvest, however, the mite must locate 97 alternative or refuge resources. This means that wild populations experience temporally 98 heterogeneous conditions , where spatial homogeneity during t he growing season is 99 intermingled with more heterogenous conditions outside of it. We thus employed replicated 100 experimental evolution to test whether temporally homogeneous environments (a single plant 101 species: either wheat or barley , Hordeum vulgare ) versus heterogeneous environments 102 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 5 (alternating these two plant species) drive the evolution of speciali sation or generalisation 103 (Skoracka et al. 2022). 104 In our experiment, a fter forty-five generations of evolution i n a stable host environment, 105 specialised phenotypes evolve d, showing improved performance on either the original host 106 (wheat) or an alternative host (barley). By contrast, a fluctuating host environment, favoured 107 the mites ’ ability to exploit multiple plant species not encountered during experimental 108 evolution, including smooth brome (Bromus inermis)—a wild grass species that A. tosichella 109 uses as a temporal refuge during early spring, autumn, and winter when cereal crops are 110 unavailable. However, mite populations cannot persist on brome in a long term (Laska et al. 111 2021; Skoracka et al. 2022; Figures 1A-C). These findings confirm that cereal pests can adjust 112 their niche breadth in response to agricultural practices. 113 Interestingly, while both barley specialists and generalists retained the ability to persist on 114 brome, wheat specialists lost this capacity (Figure 1C). This pattern suggests two interrelated 115 phenomena: (1) a trade -off between adaptation to optimal (wheat) and suboptimal 116 environments (barley, brome), and (2) a genetic overlap between the mechanisms underlying 117 adaptation to variable environments and those enabling persistence on suboptimal hosts. 118 Elucidating the genomic changes that drive speciali sation and generali sation is therefore 119 crucial to determine whether such genetic trade -offs exist and whether shared molecular 120 pathways can be identified among populations adapting to suboptimal hosts and those 121 evolving in variable environments. 122 In this study, we present a de novo genome assembly of the wheat curl mite and investigate 123 the genetic mechanisms underlying the its adaptation to different environments using an 124 evolve-and-resequence approach. By examining if there is an overlap between SNPs favoured 125 by selection on an alternative host (barley) and those favoured in the alternating h ost 126 treatment, we aim to identify candidate variants responsible for the persistence of mites on 127 alternative or refuge hosts. 128 129 2. Materials and Methods 130 2.1. Genome 131 2.1.1. Sample origin 132 The obligate phytophagous arthropod, eriophyoid mite Aceria tosichella (Keifer), commonly 133 known as the wheat curl mite (hereafter WCM) was used as a study system. WCM is a cryptic 134 species complex, consisting of several different genotypes that can be distinguished via DNA 135 barcoding (Skoracka et al. 2018a). For this study, we used the MT-1 genotype, which is a major 136 global pest of wheat (Skoracka et al. 2018b). 137 An inbred strain (isoline) of the WCM genotype MT-1 was established through ten generations 138 of successive inbreeding. The founding individuals were obtained from population collected 139 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 6 from wheat (Triticum aestivum ) in Poland (52°02'36"N, 16°46'02"E) in 2016. Mites were 140 identified taxonomically using stereomicroscopy and molecularly by barcoding the 141 cytochrome oxidase subunit I (COI) and the D2 fragment of 28S rDNA (NCBI accession 142 numbers: JF920077 and JF920097, respectively). A laboratory colony was established under 143 controlled room conditions to provide the female nymphs required for the inbreeding process. 144 WCM reproduces via arrhenotokous parthenogenesis, in which fertilised eggs develop into 145 females and unfertilised eggs develop into males (Miller et al. 2012). Consequently, a single 146 virgin female can establish a new population by producing sons with which she subsequently 147 mates. The life cycle of WCM includes the egg, larva, nymph and adult stages. Both larvae and 148 nymphs enter an immobile quiescent stage for 1-2 days prior to moulting. 149 To initiate the isoline (Generation 0; G0), a single quiescent female nymph was transferred to 150 a 4–5-day-old wheat plant housed within an isolator. After 10 –12 days, a quiescent female 151 nymph from the next generation (G1)—the offspring of a mother-son mating—was identified 152 and transferred to a fresh wheat plant. This procedure was repeated for ten generations. To 153 account for potential line collapse during inbreeding, 20 independent lines were initiated 154 simultaneously. 155 After ten generations of inbreeding, one successful colony was allowed to expand for three 156 weeks. Approximately 3,000 individuals were then harvested by cutting infested wheat leaves 157 into 15 mm fragments and placing them into 1.5 ml Eppendorf tubes filled with 70% ethanol. 158 The tubes were shaken for one minute, after which the leaf material was immediately removed 159 to prevent chlorophyll dissolution. The samples were centrifuged to form a pellet; the 160 supernatant was then discarded, and the pellet was left to air -dry. The resulting mite pellet 161 was used for genomic DNA extraction and subsequent sequencing. All colonies and 162 experimental lines were maintained at 22°C, 40% RH, and a 12/12 h L/D photoperiod. 163 2.1.2. High molecular weight genomic DNA (HMW gDNA) extraction 164 Following the evaporation of ethanol from the tube containing the mite pellet (~3000 mite 165 individuals), 180 µl of ATL buffer was added. The pellet was then homogenised using a micro-166 pestle, which was subsequently rinsed with an additional 30 µl of ATL buffer to ensure 167 maximum sample recovery . The solution was supplemented with 20 µl of P roteinase K , 168 vortexed thoroughly, and incubated at 56°C. The following day, the tube was mixed by gentle 169 inversion and centrifuged at maximum speed for 5 seconds . 200 µl of the supernatant was 170 then transferred to a new 1 .5 ml tube. To remove contaminating RNA, 4 µl of RNase A was 171 added to the supernatant , followed by a 10-minute incubation at room temperature with 172 gentle mixing. 173 After RNA digestion, 60 µl of 5M NaCl ( final concentration ~30%) was added to the sample. 174 The tube was inverted five times to precipitate the proteins and then centrifuged at 1,100 x g 175 for 15 minutes at 4°C. The supernatant containing the HMW gDNA was carefully transferred 176 to a new 1.5 ml tube using a wide-bore pipette tip. The sample was centrifuged again at 1,100 177 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 7 x g for 10 minutes at 4°C, and 180 µl of the supernatant was transferred to a 1 .5 ml tube 178 containing 360 µl of 100% ethanol using a wide-bore tip. The probe was inverted 20 times, 179 rotated for 20 minutes on a laboratory rotator (setting 2.5), and centrifuged at 6,250 x g for 5 180 minutes at 4°C for. 181 The resulting HMW gDNA p ellet aggregated on the tube wall. The s upernatant was gently 182 transferred to a separate 1.5 ml tube as a backup in case of not sufficient DNA precipitation. 183 The pellet was resuspended in 200 µl of water by gently rocking the tube 10 times. To facilitate 184 a second precipitation, 5 M NaCl was added to reach a final concentration of 0.25 M (1/20 of 185 the sample volume) and rocked ten times. Subsequently, two volumes of 100% ethanol were 186 added; the sample was gently inverted 20 times, rotated for 20 minutes (setting 2.5), and 187 centrifuged at 6 ,250 x g for 5 minutes at 4°C. The supernatant was again collected as a 188 secondary backup. The HMW gDNA pellet was air-dried at room temperature for 10 minutes 189 and then dissolved in 35 µl of TE buffer for one hour at room temperature. DNA concentration 190 was quantified using Qubit flourometer (Invitrogen). 191 2.1.3. Library construction and sequencing` 192 Libraries were prepared using the SQK-LSK109 kit for FLO-MIN106 (R9.4.1) flow cells, following 193 the manufacturer’s instructions (Oxford Nanopore Technologies) with the following 194 modifications: 1) 15 µl of gDNA dissolved in TE buffer, and 32 µl of nuclease-free water was 195 used in the library preparation step 4 (instead of 47 µl DNA suggested in a protocol) ; 2) SFB 196 was used twice to obtain all fragments of DNA at the step 13 of library preparation and 197 repeated two more times with LFB eluted to the 10 µl EB (points 15-16); 3) during step 9 of 198 priming and loading the SpotON flow cell 2,5 µl of load beads, 20 µl more of DNA library (from 199 2x10 previous steps) and 3 µl of nuclease-free water was taken. 200 Sequencing was scheduled for 72 h; however, the run was terminated after ~44 h due to a 201 significant drop in sequencing yield. The run generated 5,263,860 long reads that passed initial 202 quality control (QC) during basecalling. Following adapter removal using Porechop, 5,248,467 203 reads remained. The resulting read lengths reached up to 276 kb, with an N50 of 3,532 bp and 204 an L50 of 687,157. 205 2.1.4. Assembling and annotating the reference genome 206 Basecalling was performed using Guppy (v5.0.7) on GPU. The basecalled reads were trimmed 207 of adapter sequences using Porechop. Contigs were assembled using Flye (v 2.8.1), resulting 208 in 120 contigs. PurgeHaplotigs was run on the Flye assembly aiming at removing haplotigs. 209 Minimap2 was used to map reads to the draft primary genome assembly. Sequence polishing 210 was performed five times with racon and nanopolish software using the ONT reads . Final 211 quality control was performed with quast software. Assembly completeness was evaluated 212 using BUSCO (v5.0.8) with the arachnida_odb12 lineage dataset. 213 We sequenced transcriptomes of around 2000 adult individuals initially collected to RNAlater. 214 RNA was extracted with RNAzol and quality was assessed by TapeStation. We used a SMARTer 215 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 8 Ultra Low RNA kit and TruSeq RNA stranded library construction and then sequenced ~10Gb 216 on the HiSeq2500 with 2 x 100bp mode. Library construction and sequencing were performed 217 by Macrogen Europe. Raw reads were used to predict protein coding genes by mapping them 218 to the assembled genome with Hisat2 software and junctions were extracted using regtools. 219 Only junctions with a coverage of 10 were kept. Genomethreader was used to produce spliced 220 alignments from mite proteins on the genome. 221 Structural gene models were produced on the repeat masked genome of WCM, using 222 AUGUSTUS trained for WCM, infused with hints from transcriptome read-junctions and 223 protein alignments (from other mites). These prediction s were further improved with 224 EvidenceModeler (EVM) by integrating augustus predictions (primary prediction), transcript 225 alignments, and protein homology evidence. Gene structures were further refined and 226 updated using PASA , aiming at extending UTRs . Annotation quality and consistency were 227 monitored using BUSCO throughout the workflow. 228 2.2. Experimental evolution 229 The stock population of Aceria tosichella MT-1 was established in 2017 by pooling 26 field -230 derived populations collected from 10 wheat fields across Poland (see Skoracka et al. 2022 for 231 details). Genotypic identity was confirmed via COI barcoding, and colonies were periodically 232 validated to ensure the absence of contamination from non-MT-1 genotypes. Approximately 233 1,000 individuals from each field -derived population were combined to maximise genetic 234 diversity. The resulting stock population was maintained on wheat prior to the 235 commencement of experimental evolution. 236 The experimental evolution was conducted for 45 generations under three host-plant selection 237 regimes, each replicated four times: (i) constant wheat, Triticum aestivum (cT line); (ii) constant 238 barley, Hordeum vulgare (cH line); and (iii) alternating Triticum and Hordeum (aTH line). This 239 design established two regimes with constant biotic conditions and one with fluctuating biotic 240 conditions. 241 Each replicate population was founded using approximately 300 individuals transferred from 242 the stock population onto clean potted plants (20 plants per pot) using an aspirator . 243 Populations were maintained in a growth chamber at 27 °C, with a 16:8 light/dark (L/D) cycle 244 and 60% relative humidity (RH). Every three weeks, approximately 300 mite individuals from 245 each replicate population were transferred to a new pot containing 20 new clean plants 246 according to their respective selection regimes. As A. tosichella development from egg to egg 247 is temperature-dependent, a three-week period at 27 °C corresponds to approximately three 248 generations (Karpicka-Ignatowska et al. 2021). Consequently, the 45-generation experiment 249 spanned 15 such transfer cycles. 250 2.3. Resequencing and genomic analyses 251 2.3.1. Genomic sampling and mapping 252 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 9 For genomic analyses material was sampled from populations in all three evolutionary regimes 253 (cT, cH, and aTH) at generations 15, 30, and 45, as well as from the stock population prior to 254 the start of the experimental evolution. Approximately 1 ,000 individuals per sample were 255 preserved in 70% ethanol and stored at -20°C. Total genomic DNA was extracted using the 256 DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) following the protocol described by 257 Dabert et al. (2008). DNA was quantified using a Qubit fluorometer (Invitrogen) and stored at 258 -20 °C before being dispatched to Macrogen Europe for library preparation and sequencing. 259 Libraries were prepared using the Illumina TruSeq Nano “low-input” protocol (350 bp insert 260 size). Sequencing was performed on the Illumina NovaSeq 6000 platform using an S4 flow cell 261 in 2 x 150 bp mode. Raw reads were trimmed using Trimmomatic ( v.0.39; Bolger et al. 2014); 262 unpaired reads were discarded from analyses. Fastq files were mapped to the assembled 263 genome using BWA-MEM (v0.7.10-r789; Li, 2013) with default settings. Resulting SAM files 264 were converted to BAM format, sorted, and duplicates were marked. To ensure high mapping 265 quality, reads with a mapping quality (MQ) score below 20 were removed using SAMtools 266 (v1.16; Li et al. 2009) and Picard Tools. Mapping s tatistics were assessed using Qualimap 267 (García-Alcalde et al. 2012). 268 2.3.2. SNP Identification and estimation of nucleotide diversity 269 Bam files were converted to pileup format using SAMtools. To avoid false-positive SNPs, indels 270 and their surrounding regions (5bp either side) were identified and filtered using identify -271 genomic-indel-regions.pl and filter -pileup-by-gtf.pl scrpts from the PoPoolation package 272 (v.1.2.2; Kofler et al. 2011 a). The filtered pileup file was used to determine the coverage 273 distribution by sampling every 10,000 line. Minimum and maximum coverage thresholds were 274 set at 60$ \times$ and 240$ \times$, respectively, based on the mean coverage across all 275 samples. These thresholds were applied using a custom Python script to filter the pileup file. 276 The final mpileup file was converted to a sync format using mpileup2sync.jar from 277 PoPoolation2 (Kofler et al. 2011b). Polymorphic SNPs were identified using the snp-frequency-278 diff.pl script (PoPoolation2) with the parameters --min-count 6 and --min-coverage 40. SNPs 279 meeting these criteria were then extracted from the sync file via a custom Python script. 280 To estimate nucleotide diversity the final mpileup file was used. For each sample, the relevant 281 columns were extracted and processed using the Variance-sliding.pl script (Popoolation) to 282 calculate nucleotide diversity ( 𝜋) in 20kb windows (10kb step size) with –min-count 3 –min-283 coverage 15, --pool-size 500. As the data included estimates of 𝜋 from the same population 284 across multiple generations, we performed analysis by repeated measures analysis of variance 285 (ANOVA), which takes into account the non -independence of samples , following the 286 methodology of Parrett et al. (2022). Comparisons were performed using the mean values for 287 each experimentally evolved population. Additionally, Tajima’s D was calculated for the stock 288 population in 20 kb non-overlapping windows using --min-count 2 and --min-coverage 50. 289 2.3.3. Identification of differentiated SNPs 290 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 10 The sync file containing the identified SNPs was first explored with cvtk Python package 291 (Buffalo and Coop 2020). After integrating the genomic data with metadata and the annotation 292 file, the stock population was multiplied and defined as “generation zero” for each selection 293 regime. This step had no future consequences for analyses, but it was required for some steps 294 in data exploration (not shown). Number of fixed variants within samples and associations 295 between coverage and calculated variance/covariance were then calculated following the 296

of Buffalo and Coop (2020). Data arrays were then used for future filtering, retaining 297 only SNPs with minor allele frequency (averaged across all samples) higher than 5%. The 298 frequencies of the major alleles (defined across all samples) were written to a file and analysed 299 using R script. Each SNP having coverage between 60 and 240x in all samples were used to 300 perform three GLMs, comparing pairs of experimental treatments in each time point 301 separately. We compared the counts of the major alleles against counts of minor alleles to 302 determine consistent allele frequency changes between treatments. For all samples +1 was 303 added to minor and major allele counts and qu asibinomial error structure was used for the 304 model. SNPs were defined as differentiated between treatments, if GLM-derived p-value was 305 smaller than 0.001. Custom R scripts were used to visualize results. 306 To identify the functional context of significant SNPs, we computed overlaps between 307 candidate SNPs and annotated genes using the findOverlaps function from the 308 GenomicRanges package. This process was performed separately for each treatment dataset. 309 Genes were compared both against all SNPs in a comparison (background set) and against 310 differentiated SNPs (foreground set). Gene identifiers and their annotations were extracted 311 from the overlapping annotation elements. 312 2.3.4. Functional analyses of differentiated SNPs 313 Gene Ontology (GO) enrichment analysis was performed to characteri se the biological 314 functions associated with genes containing significant SNPs. GO annotations were derived 315 from a custom mapping file retrieved from the OrcAE platform. This annotation was parsed 316 and formatted to create a comprehensive gene-to-GO term mapping. 317 Enrichment analyses were conducted using the topGO package in R, performed separately for 318 each of the three ontology domains: Biological Process (BP), Molecular Function (MF), and 319 Cellular Component (CC). For each pairwise comparison (cT vs. cH, aTH vs. cT, and aTH vs. cH), 320 a binary vector indicating the presence or absence of a gene in the foreground set was used 321 to construcy a topGOdata object. The weighted Fisher's exact test (algorithm = "weight01") 322 was applied to assess statistical overrepresentation. The top 10 significantly enriched GO terms 323 (Fisher p-value < 0.05) were extracted and saved for reporting. 324 To verify the expression of specific genes harbouring differentiated SNPs, we used generated 325 by us RNA-Seq data (see above) combined with data from Gupta et al. (2019). Raw reads were 326 trimmed with Trimmomatic in the same way as described above for resequenced data. STAR 327 (v2.7.6a; Dobin et al. 2013) was then employed to map transcriptomic reads to the reference 328 genome, retaining only reads contain ing junctions that satisfied the filtering criteria ( --329 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 11 outFilterType BySjout). We allowed a minimum block size for spliced alignments to be 8 ( --330 alignSJoverhangMin 8), a maximum mismatch ratio of 0.04 ( -331 outFilterMismatchNoverReadLmax 0.04), and a range for i ntron length and mate-pair gaps 332 between 20 bp and 100 ,000 bp. The resulting BAM files where then visualised alongside the 333 genome annotation using IGV software (Thorvaldsdóttir et al. 2013) to provide qualitative 334 confirmation of gene expression and splice site accuracy. 335 2.4. Frequency changes of candidate SNPs on brome 336 To investigate whether haplotypes carrying two SNPs—identified as significantly differentiated 337 in all three pairwise comparisons —mediate differences in the mites' ability to maintain 338 positive growth on smooth brome ( Bromus inermis ; see Skoracka et al . 2022, Fig. 1C), we 339 monitored frequency shifts of these candidate SNPs over three generations of exclusive 340 maintenance on brome. This experiment was conducted approximately three years following 341 the initial whole-genome resequencing . We used the aTH populations, which had been 342 continuously maintained under alternating wheat and barley conditions. We focused on aTH 343 populations because the SNPs segregated in them at more intermediate frequencies than in 344 cT and cH (see Results, Figure 1). 345 We first extracted DNA from 90 mite individuals collected from three aTH lines (G10, G11, G18; 346 30 individuals per each line) using the Chelex protocol (Bouneb et al. 2014). A 445 bp fragment 347 of Contig_37, containing the two candidate SNPs, was amplified using the primers 5’ -348 GTGCATGCATGTGCTCTACC-3’ and 5’ -GAATGGTGCTCACCAAGCG-3’ and characterised via 349 Sanger sequencing. Subsequently, 3,000 mite individuals from each of the three aTH 350 populations were transferred from wheat to brome using an aspirator. Six replicates were 351 established for each population (i.e. 6 x G10, 6 x G11, 6 x G18) . Colonies were maintained at 352 26 °C, 16:8 h L/D cycle and 60% RH for three weeks, corresponding to approximately three 353 generations. At the end of this period, we collected all surviving individuals and sequenced the 354 same genomic fragment to asses allele frequency changes. To compare allele frequencies 355 across lines and time points, we employed a Generalised Linear Model (GLM) using allele 356 counts as the response variable, with generation and population included as fixed variables. 357 2.5. Code availability 358 Scripts used to analyse data are available on 359 https://github.com/konczal/WCM_EvolveResequence 360 361 3. Results 362 3.1. Reference genome 363 Long-read Nanopore sequencing generated 6.7 million reads with an N50 of ~3.7 kb and a total 364 yield of 11.6 Gb. These data were used to construct a curated de novo reference genome for 365 Aceria tosichella. The final assembly comprises 88 contigs with a total length of 44.84 Mb and 366 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 12 a mean coverage of 194x. Contig sizes range from 105 bp to 9,100,975 bp, with a mean length 367 of 509.6 kb and a median of 15.5 kb. The GC content across contigs spans 25.6–60.0%, with a 368 genome-wide average of 46.3%, consistent with other eriophyoid mites. These steps resulted 369 in a high-quality curated reference genome, providing a robust foundation for all downstream 370 analyses. 371 Genome completeness, assessed using BUSCO arachnida_od12 dataset (n = 1123) revealed 372 59.8% complete genes, including 1.6% duplicated, 8.4% fragmented and 31.9% missing. 373 Parallel BUSCO analyses performed on a genome of another eriophyoid mite , Aculops 374 lycopersici (Greenhalgh et al. 2020) , yielded comparable results ( Table S1 in Supporting 375 Information 2). The AUGUSTUS-EVM–PASA annotation pipeline predicted 8,363 genes, with an 376 average density of 186 genes per Mb. Of these, 4,819 were single-exon genes, and the average 377 exon length was 998 bp. Detailed summary statistics and a comparison with A. lycopersici are 378 provided in the Tab le S2 in Supporting Information 2 . The sequence and full annotation are 379 publicly available via t he ORCAE platform: 380 https://bioinformatics.psb.ugent.be/orcae/overview/Aceto 381 3.2. Genetic differentiation in experimentally evolved populations 382 On average, 97% (range: 88–100%) of reads from the resequenced populations successfully 383 mapped to the reference genome. All samples were sequenced at high coverage, with an 384 average depth of 147× (range: 127–158×; Table S3 in Supporting Information 2), enabling the 385 identification of 625,667 SNPs —approximately one SNP per 10,000 bp. Mean nucleotide 386 diversity in the stock population was 1.0 × 10⁻³, an order of magnitude higher than in any of 387 the evolved populations (range: 2.1 –4.1 × 10⁻⁴; Figure S1 in Supporting Information 1 ). 388 Genome-wide nucleotide diversity did not significantly differ between selection regimes (df = 389 2, F = 3.2, p = 0.10), generations (df = 1, F = 1.7, p = 0.24), or their interaction (df = 2, F = 0.07, 390 p = 0.94). 391 Comparison of allele frequencies revealed significant differentiation in 237 SNPs between cT 392 and cH line ages, 275 SNPs between cT and aTH, and 191 SNPs between cH and aTH (Figure 393 1D). Among the associated Gene Ontology (GO) terms, inositol 1,4,5-trisphosphate binding 394 (GO:0070679) was significantly enriched across all three pairwise comparisons. ABC -type 395 transporter activity emerged as a recurrent functional category in both the aTH vs. cT and aTH 396 vs. cH comparisons, reflected by overlapping enriched GO terms. Additional significantly 397 enriched terms included homophilic cell adhesion via plasma membrane adhesion molecules 398 (aTH vs. cH), double-strand break repair (aTH vs. cT), and xenobiotic metabolic process (cT vs. 399 cH). A complete list of enriched GO terms is provided in the Table 1. 400 Two SNPs were consistently differentiated across all three pairwise comparisons (Figures 2D, 401 2E), while eight SNPs showed differentiation between the alternating and constant 402 environments (aTH vs. cH/cT; Figure S3 in Supporting Information 1 ). Additionally, 50 SNPs 403 were differentiated between the wheat-evolved populations (cT) and both other treatments 404 (aTH/cH; Figure S2 in Supporting Information 1). This relatively high number of differentiated 405 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 13 SNPs in comparisons involving cT appears to be largely driven by a single genomic region 406 (highlighted in red in Figures 1E–G). Of these 50 SNPs, 42 were located on contig_37, with 33 407 clustered within a 120-kbp region (contig_37:2,870,000–2,990,000). 408 3.3. Genomic region with the largest differentiation 409 The genomic region containing the majority of SNPs differentiated between cT and cH/aTH 410 lineages spans 13 protein -coding genes (Fig ures 2G–I; Table 2). In the stock population, this 411 region exhibits elevated nucleotide diversity (1.57 × 10⁻³) compared to the genome -wide 412 median (0.91 × 10⁻³; Figure 2B), along with a slightly less negative Tajima’s D (–2.27 vs. –2.40; 413 Figure 2C), indicating increased genetic variation and a higher proportion of intermediate -414 frequency alleles relative to the rest of the genome. 415 Within this region, many SNPs show a pronounced decline in major allele frequency in cT 416 lineages, a subtle decline in cH lineages, and a slight increase in aTH lineages (Figure 2A). This 417 pattern is consistent across most SNPs differentiated between cT and cH/aTH ( Figure S2 in 418 Supporting Information 1 ). In contrast, SNPs differentiated in the aTH vs. cT/cH comparison 419 (Figure S3 in Supporting Information 1 ) and in the cH vs. cT/aTH comparison ( Figure S4 in 420 Supporting Information 1) exhibit distinct allele frequency dynamics, except for those located 421 within the same region of contig_37. 422 The region also includes the two SNPs differentiated across all three selection regimes (Figures 423 2D, 2E). Although initial automated annotation placed them in an intergenic region, RNA -Seq 424 read mapping revealed that they are actually located within the intronic region of the 5′ UTR 425 of the aceto37g05650 gene (Figure S5 in Supporting Information 1). This gene shows strong 426 homology to Tetranychus urticae gene tetur01g16794 (BLASTp E-value: 1.6e–81), previously 427 identified as encoding the Allatostatin C receptor 2 (AstC -R2) (Veenstra et al. 2012). Another 428 SNP in this region, also differentiated between alternating and constant environments (i.e., 429 aTH vs. cT/cH; Figure 2F), is located within the aceto37g05600 gene, which is annotated as an 430 inositol 1,4,5 -trisphosphate receptor (ITPR1). All 13 genes identified within this genomic 431 region, along with their putative molecular functions, are listed in Table 2. 432 Sanger-sequenced and genotyped two SNPs located in the 5′ UTR of the aceto37g05650 gene 433 in a total of 94 individuals derived from three aTH lines (sampled ~ 3 years after the genomic 434 analyses). Of these, 74 individuals were sampled at the start of the experiment, and 20 were 435 collected after three generations on the refuge host ( smooth brome). Allele frequencies did 436 not differ among the three populations (df = 2, χ² = 3.28, p = 0.19), nor did they differ between 437 the two time points of the experiment (df = 1, χ² = 0.37, p = 0.54). 438 439 4. Discussion 440 Heterogeneous environments —such as those created by crop rotation or seasonal 441 harvesting—are predicted to favo ur the evolution of generalist strategies, defined by the 442 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 14 ability to exploit multiple resources (Futuyma and Moreno 1988; Kassen 2002; Levins 1968; 443 Sexton et al. 2017). The evolution of such generalism has profound ecological and agricultural 444 implications; therefore understanding the underlying molecular mechanisms is essential for 445 informing pest management. However, the genetic architecture that facilitates adaptation to 446 specific hosts versus the maintenance of broad niche breadth remains relatively poorly 447 understood (Shih et al. 2023). 448 Here, we characterise the molecular basis of generalism in Aceria tosichella, a global cereal 449 pest and viral vector. We identified single nucleotide polymorphisms (SNPs) associated with 450 selection for specialisation and generalism, following up on trade -offs previously reported 451 between host-specific performance and broad resource utilisation (Skoracka et al. 2022). We 452 subsequently carried out a follow-up experiment and functional analysis of candidate SNPs 453 underlying the trade-offs we previously reported between improved specialisation to a stable 454 environment and the utilisation a of other plant species. The trade-off was particularly striking 455 for cT populations (evolving on their original host, wheat), which exhibited a negative growth 456 rate on a refuge host, brome, whereas the growth rate of cH populations ( evolving on an 457 alternative host, barley) on brome was intermediate between aTH and c T (Figure 1C). This 458 suggests that the ability to survive on brome was lost in cT lineages due to a genetic trade-off. 459 Such ability was retained in aTH and, to a lesser extent, in cH line ages, likely via the 460 maintenance of ‘ generalist variants’ enabling persistence on refuge hosts. We therefore 461 expected that SNPs which differ between specialist s and generalist, particularly between cT 462 and aTH, would be enriched for functional categories associated with generalism. 463 We initially focused on SNPs showing differentiation in allele frequencies among the three 464 treatments (cH, cT, aTH) and detected hundreds of such loci in each pairwise comparison. 465 However, only two SNPs displayed consistent differentiation across all three comparisons. Both 466 loci were located in close proximity within the intronic regions of the AstC-R2 gene, whose 467 paralog in Drosophila mediates intestinal nutrient signalling (Kubrak et al. 2022). The loss of 468 this receptor impairs lipid and sugar mobili sation during fasting, leading to hypoglyc aemia 469 while substantially increasing starvation resistance (Kubrak et al. 2022). Notably, changes in 470 allele frequencies at these SNPs (Figures 2D, 2E) closely paralleled the variation in population 471 growth rate s measured on the refuge host, smooth brome (Fig ure 1C). We therefore 472 hypothesised that these polymorphisms might be associated with fitness on this plant species. 473 Specifically, the positive population growth rate on brome observed in the aTH lineages might 474 reflect enhanced starvation resistance arising from genetic changes in the AstC-R2 gene. 475 However, results from the subsequent brome-derived starvation experiment did not support 476 this hypothesis. Allele frequencies at AstC-R2 variants remained unchanged among individuals 477 that survived on brome, suggesting weak or absent selection on these variants during 478 exposure to brome. We thus conclude that the frequency changes in AstC-R2 variants in 479 response to treatments involving wheat and/or barley were unlikely to be a major driver of 480 the correlated changes in the ability to utilise brome. Instead, these correlated changes may 481 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 15 stem from other SNPs that shifted significantly in frequency within the cT lineages, particularly 482 when compared to aTH lineages, which exhibited superior performance on brome. 483 Interestingly, most SNPs that exhibited frequency differentiation between cT and aTH—as well 484 as between cT and cH—are located in the genomic region adjacent to AstC-R2, which harbours 485 numerous differentiated SNPs showing similar allele frequency shifts (Figure 2A, Figure S2 in 486 Supporting Information 1). This genomic region is characterised by increased genetic variation 487 in the base population (Figure 2B); however Tajima’s D value for this region was moderate 488 rather than extreme compared to the rest of the genome (Fig ure 2C), suggesting that the 489 observed polymorphism is not driven by balancing selection. 490 Notably, the frequency of major alleles in this region tends to decrease in cT lineages, remain 491 largely unchanged in cH lineages, and slightly increase or remain stable in aTH lineages (Figure 492 S2 in Supporting Information 2), mirroring the pattern of changes in growth rates on brome in 493 a manner similar to that of AstC-R2. This suggests that one or more genes in this region—494 which encompasses 13 protein -coding sequences —may contribute to the generalist 495 phenotypes capable of surviving on brome. One of these SNPs (Fig ure 2F) falls within the 496 coding sequence of the inositol 1,4,5-trisphosphate receptor (ITPR1). The loss of IP3R function 497 in Drosophila also has impo rtant physiological consequences , leading to adult obesity and 498 increased starvation resistance (Subramanian et al. 2013). Knockdown flies exhibit reduced 499 metabolism of long -chain fatty acids and impaired appetite control. Several other notable 500 genes reside within this region. For example, the 1,4-alpha-glucan branching enzyme is 501 essential for glycogen synthesis and cellular iron homeostasis (Huynh et al. 2019). Another is 502 the precursor of cathepsin L, a key cysteine protease involved in fat body catabolism (Yang et 503 al. 2020), which is primarily localised in the midgut, where it participates in the degradation 504 of both food and foreign pathogens (Cristofoletti et al. 2003). There is also a nicotinic 505 acetylcholine receptor subunit, which is a part of the neurotransmitter receptor family 506 targeted by many insecticides (Ihara et al. 2020). Consequently, variation in these genes might 507 contribute to the broadening of the ecological niche breadth observed in aTH line ages. We 508 note that since AstC-R2 is localised within the same region, changes in its variant frequencies 509 resulting from differential survival on brome should be correlated with changes of other SNPs 510 in the region . However, our brome survival experiment was performed three years 511 (approximately 100 generations) after genetic samples were collected, a period which may 512 have broken down linkage disequilibrium (LD) within the region. Therefore, these genes should 513 not be dismissed as candidates for further functional investigation, despite the lack of support 514 for the role of AstC-R2 in the current experiment. 515 There were also numerous differentiated SNPs located across other genomic regions, 516 indicating that the evolution of generalism has a polygenic basis. For example, hexokinase-2 517 (aceto07g03060), a key regulator of energy metabolism in insect (Lin and Xu 2016) , 518 agglutinine (aceto11g12900), which plays roles in immunity and cell recognition, and ephrin 519 type-B receptor (aceto41g02960), involved in nervous system development (Boyle et al. 520 2006), all showed differentiated allele frequencies between the aTH and cT populations. 521 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 16 Similarly, differences between the aTH and cH lineages involved genes such as nuclear 522 receptor corepressor 1 (aceto156g00200), ATP-binding cassette sub -family A member 3 523 (aceto41g00020), sub-family A member 2 (aceto41g00040), and a retinal-specific ATP-binding 524 cassette transporter (aceto41g00060), all of which are associated with xenobiotic 525 detoxification and insecticide resistance (Birnbaum and Abbot 2020). Additionally, variants in 526 GMC oxidoreductase ( aceto46g01390), a gene potentially involved in manipulating plant 527 defences (Lin et al. 2023) were also detected. Together, these patterns support the widely 528 held view that ecologically relevant traits are typically polygenic, often producing only subtle 529 genomic signals (Barghi et al. 2020). 530 Our analysis of Gene Ontology (GO) term enrichment identified key molecular processes and 531 functions shaping the niche breadth in A. tosichella. For example, we detected an enrichment 532 of SNPs differentiated between aTH populations and other treatments among genes 533 associated with ABC-type transporter activity . ABCs are canonical detoxification genes and 534 often exhibit greater regulated plastic responses to novel hosts (Birnbaum and Abbot 2020). 535 Another enriched GO term is xenobiotic transport, which is crucial for processes such as stress-536 related detoxification and has been demonstrated to mediate aphid colonisation of previously 537 resistant soy varieties (Bansal et al. 2014). Additionally, we found enrichment in inositol 1,4,5-538 trisphosphate-gated calcium channel activity , which is involved in regulation of diverse 539 physiological responses (Agrawal et al. 2009). These finding support the view that generalist 540 herbivores require a comprehensive detoxification metaboli c system to adapt to a range of 541 plant species (Van Leeuwen and Dermauw 2016). 542 In conclusion, our analysis identified several SNPs associated with the selection response 543 towards specialisation and generalism, as well as the correlated response related to niche 544 breadth. Significantly diverged SNPs were enriched in highly polymorphic genomic region on 545 contig 37. However, the follow-up experimental testing of these candidate SNPs indicated that 546 they are not directly associated with niche breadth, defined here as the ability to sustain 547 populations on a refuge plant species that was un available during experimental evolution . 548 Other candidate SNPs present within this region warrant further scrutiny in future research. In 549 addition to SNPs within the region, we identified numerous SNPs with diverse molecular 550 functions across other genomic regions, supporting the view that host adaptation is typically 551 polygenic. Finally, functional analysis revealed that diverged SNPs are enriched for several 552 functional categories, particularly those associated with detoxification. This highlights that 553 generalists require an all-purpose detoxification metabolism allowing them to persist across a 554 broad range of plant species and adapt to novel or variable environments. Overall, our results 555 advance the understanding of the molecular basis underlying the evolution of generalist 556 strategies and niche breadth in the cereal pest Aceria tosichella within heterogeneous host 557 environments. 558 559 Acknowledgments 560 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 17 We are grateful to Anna Radwańska and Kamila Zalewska for their contributions to 561 experimental evolution and material collecting, to Wiktoria Szydło and Monika Stefańska for 562 their help assistance in maintaining the colonies, to Sebastian Chmielewski for his help with 563 the DNA isolation protocol, and to Lechosław Kuczyński for his valuable suggestions at the 564 initial stage of drafting the manuscript. We would also like to thank the DANKO Hodowla Roślin 565 company for the Triticum aestivum and Hordeum vulgare seeds, and the CENTNAS Sp. z o.o. 566 company in Krotoszyn, Poland, and the Botanical Garden in Bydgoszcz, Poland, for the Bromus 567 inermis seeds. This study was supported by the National Science Centre (NSC), Poland, 568 research grants no. 2017/27/N/NZ8/00305 and 2019/32/T/NZ8/00151, which were awarded 569 to Alicja Laska. The experiment measuring frequency changes of candidate SNPs in smooth 570 brome was supported by the National Science Centre (NSC), Poland, research grant no. 571 2021/41/B/NZ8/01703 awarded to Anna Skoracka. 572 573

574 Agrawal, N., Padmanabhan, N., & Hasan, G. (2009). Inositol 1, 4, 5 -trisphosphate receptor 575 function in Drosophila insulin producing cells. PLoS ONE , 4(8), e6652. 576 https://doi.org/10.1371/journal.pone.0006652 577 Altieri, M. A. (1999). The ecological role of biodiversity in agroecosystems. Agriculture, 578 Ecosystems & Environment, 74(1-3), 19-31. https://doi.org/10.1016/S0167-579 8809(99)00028-6 580 Andow, D. A. (1991). Vegetational diversity arthropod population. Annual Review of 581 Entomology, 36, 561–586. https://doi.org/10.1146/annurev.en.36.010191.003021 582 Bansal, R., Mian, M. A. R., Mittapalli, O., & Michel, A. P. (2014). RNA -Seq reveals a xenobiotic 583 stress response in the soybean aphid, Aphis glycines, when fed aphid-resistant soybean. 584 BMC Genomics, 15(1), 972. https://doi.org/10.1186/1471-2164-15-972 585 Barghi, N., Hermisson, J., & Schlötterer, C. (2020). Polygenic adaptation: a unifying framework 586 to understand positive selection. Nature Reviews Genetics , 21(12), 769 –781. 587 https://doi.org/10.1038/s41576-020-0250-z 588 Birnbaum, S. S. L., & Abbot, P. (2020). Gene expression and diet breadth in plant -feeding 589 insects: summarizing trends. Trends in Ecology & Evolution , 35(3), 259 –277. 590 https://doi.org/10.1016/j.tree.2019.10.014 591 Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina 592 sequence data. Bioinformatics, 30(15), 2114 –2120. 593 https://doi.org/10.1093/bioinformatics/btu170 594 Bono, L. M., Draghi, J. A., & Turner, P. E. (2020). Evolvability costs of niche expansion. Trends 595 in Genetics, 36(1), 14–23. https://doi.org/10.1016/j.tig.2019.10.003 596 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 18 Bono, L. M., Smith Jr, L. B., Pfennig, D. W., & Burch, C. L. (2017). The emergence of 597 performance trade-offs during local adaptation: insights from experimental evolution. 598 Molecular Ecology, 26(7), 1720–1733. https://doi.org/10.1111/mec.13979 599 Bouneb, M., de Lillo, E., Roversi, P. F., & Simoni, S. (2014). Molecular detection assay of the 600 bud mite Trisetacus juniperinus on Cupressus sempervirens in nurseries of central Italy. 601 Experimental and Applied Acarology , 62(2), 161–170. https://doi.org/10.1007/s10493-602 013-9733-3 603 Boyle, M., Nighorn, A., & Thomas, J. B. (2006). Drosophila Eph receptor guides specific axon 604 branches of mushroom body neurons. Development, 133(9), 1845 –1854. 605 https://doi.org/10.1242/dev.02353 606 Buffalo, V., & Coop, G. (2020). Estimating the genome -wide contribution of selection to 607 temporal allele frequency change. Proceedings of the National Academy of Sciences , 608 117(34), 20672–20680. https://doi.org/10.1073/pnas.1919039117 609 Condon, C., Cooper, B. S., Yeaman, S., & Angilletta Jr, M. J. (2014). Temporal variation favors 610 the evolution of generalists in experimental populations of Drosophila melanogaster. 611 Evolution, 68(3), 720–728. https://doi.org/10.1111/evo.12296 612 Cristofoletti, P. T., Ribeiro, A. F., Deraison, C., Rahbé, Y., & Terra, W. R. (2003). Midgut 613 adaptation and digestive enzyme distribution in a phloem feeding insect, the pea aphid 614 Acyrthosiphon pisum . Journal of Insect Physiology , 49(1), 11 –24. 615 https://doi.org/10.1016/S0022-1910(02)00222-6 616 Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., … Gingeras, T. R. (2013). 617 STAR: Ultrafast universal RNA -seq aligner. Bioinformatics, 29(1), 15 –21. 618 https://doi.org/10.1093/bioinformatics/bts635 619 Egas, M., Dieckmann, U., & Sabelis , M. W. (2004). Evolution restricts the coexistence of 620 specialists and generalists: the role of trade -off structure. The American Naturalist , 621 163(4), 518–531. https://doi.org/10.1086/382599 622 Futuyma, D. J., & Moreno, G. (1988). The evolution of ecological specialization. Annual Review 623 of Ecology and Systematics , 19, 207 –233. 624 https://doi.org/10.1146/annurev.es.19.110188.001231 625 García-Alcalde, F., Okonechnikov, K., Carbonell, J., Cruz, L. M., Götz, S., Tarazona, S., … Conesa, 626 A. (2012). Qualimap: Evaluating next -generation sequencing alignment data. 627 Bioinformatics, 28(20), 2678–2679. https://doi.org/10.1093/bioinformatics/bts503 628 Gompert, Z., Jahner, J. P., Scholl, C. F., Wilson, J. S., Lucas, L. K., Soria-Carrasco, V., … Forister, 629 M. L. (2015). The evolution of novel host use is unlikely to be constrained by trade-offs 630 or a lack of genetic variation. Molecular Ecology , 24(11), 2777 –2793. 631 https://doi.org/10.1111/mec.13199 632 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 19 Greenhalgh, R., Dermauw, W., Glas, J. J., Rombauts, S., Wybouw, N., Thomas, J., … Feyereisen, 633 R. (2020). Genome streamlining in a minute herbivore that manipulates its host plant. 634 eLife, 9, e56689. https://doi.org/10.7554/eLife.56689 635 Gupta, A. K., Scully, E. D., Palmer, N. A., Geib, S. M., Sarath, G., Hein, G. L., & Tatineni, S. (2019). 636 Wheat streak mosaic virus alters the transcriptome of its vector, wheat curl mite (Aceria 637 tosichella Keifer), to enhance mite development and population expansion. Journal of 638 General Virology, 100(5), 889–910. https://doi.org/10.1099/jgv.0.001256 639 Huynh, N., Ou, Q., Cox, P., Lill, R., & King-Jones, K. (2019). Glycogen branching enzyme controls 640 cellular iron homeostasis via Iron Regulatory Protein 1 and mitoNEET. Nature 641 Communications, 10(1), 5463. https://doi.org/10.1038/s41467-019-13237-8 642 Ihara, M., Furutani, S., Shigetou, S., Shimada, S., Niki, K., Komori, Y., … Hikida, M. (2020). 643 Cofactor-enabled functional expression of fruit fly, honeybee, and bumblebee nicotinic 644 receptors reveals picomolar neonicotinoid actions. Proceedings of the National 645 Academy of Sciences , 117(28), 16283 –16291. 646 https://doi.org/10.1073/pnas.2003667117 647 Karpicka-Ignatowska, K., Laska, A., Rector, B. G., Skoracka, A., & Kuczyński, L. (2021). 648 Temperature-dependent development and survival of an invasive genotype of wheat 649 curl mite, Aceria tosichella . Experimental and Applied Acarology , 83(4), 513 –525. 650 https://doi.org/10.1007/s10493-021-00602-w 651 Kassen, R. (2002). The experimental evolution of specialists, generalists, and the maintenance 652 of diversity. Journal of Evolutionary Biology , 15(2), 173 –190. 653 https://doi.org/10.1046/j.1420-9101.2002.00377.x 654 Kennedy, G. G., & Storer, N. P. (2000). Life systems of polyphagous arthropod pests in 655 temporally unstable cropping systems. Annual Review of Entomology , 45(1), 467–493. 656 https://doi.org/10.1146/annurev.ento.45.1.467 657 Ketola, T., Mikonranta, L., Zhang, J. I., Saarinen, K., Örmälä , A.-M., Friman, V.-P., … Laakso, J. 658 (2013). Fluctuating temperature leads to evolution of thermal generalism and 659 preadaptation to novel environments. Evolution, 67(10), 2936 –2944. 660 https://doi.org/10.1111/evo.12148 661 Kofler, R., Orozco -terWengel, P., De Maio, N., Pandey, R. V., Nolte, V., Futschik, A., … 662 Schlötterer, C. (2011a). PoPoolation: a toolbox for population genetic analysis of next 663 generation sequencing data from pooled individuals. PLoS ONE , 6(1), e15925. 664 https://doi.org/10.1371/journal.pone.0015925 665 Kofler, R., Pandey, R. V., & Schlötterer, C. (2011b). PoPoolation2: Identifying differentiation 666 between populations using sequencing of pooled DNA samples (Pool -Seq). 667 Bioinformatics, 27(24), 3435–3436. https://doi.org/10.1093/bioinformatics/btr589 668 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 20 Kubrak, O., Koyama, T., Ahrentløv, N., Jensen, L., Malita, A., Naseem, M. T., ... & Rewitz, K. 669 (2022). The gut hormone Allatostatin C/Somatostatin regulates food intake and 670 metabolic homeostasis under nutrient stress. Nature communications , 13(1), 671 692. https://doi.org/10.1038/s41467-022-28268-x 672 Laska, A., Magalhães, S., Lewandowski, M., Puchalska, E., Karpicka-Ignatowska, K., Radwańska, 673 A., … Skoracka, A. (2021). A sink host allows a specialist herbivore to persist in a seasonal 674 source. Proceedings of the Royal Society B , 288(1958), 20211604. 675 https://doi.org/10.1098/rspb.2021.1604 676 Leonard, A. M., & Lancaster, L. T. (2022). Evolution of resource generalism via generalized 677 stress response confers increased reproductive thermal tolerance in a pest beetle. 678 Biological Journal of the Linnean Society , 137(2), 374 –386. 679 https://doi.org/10.1093/biolinnean/blac082 680 Levins, R. (1968). Evolution in changing environments: some theoretical explorations . 681 Princeton University Press. 682 Li, H. (2013). Aligning sequence reads clone sequences and assembly contings with BWA -683 MEM. arXiv Preprint. https://doi.org/10.48550/arXiv.1303.3997 684 Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., … Durbin, R. (2009). The 685 Sequence Alignment/Map format and SAMtools. Bioinformatics, 25(16), 2078 –2079. 686 https://doi.org/10.1093/bioinformatics/btp352 687 Lin, X.-W., & Xu, W. -H. (2016). Hexokinase is a key regulator of energy metabolism and ROS 688 activity in insect lifespan extension. Aging (Albany NY) , 8(2), 245 –259. 689 https://doi.org/10.18632/aging.100885 690 Lin, Y.-H., Silven, J. J. M., Wybouw, N., Fandino, R. A., Dekker, H. L., Vogel, H., … Haring, M. A. 691 (2023). A salivary GMC oxidoreductase of Manduca sexta re-arranges the green leaf 692 volatile profile of its host plant. Nature Communications , 14(1), 3666. 693 https://doi.org/10.1038/s41467-023-39353-0 694 Litovska, I., van der Plas, F., Buijs, G., Alexandre, N., & Kleijn, D. (202 6). Nature -inclusive 695 farming results in higher arthropod abundance, without compromising agricultural 696 productivity. Agriculture, Ecosystems & Environment , 397, 110033. 697 https://doi.org/10.1016/j.agee.2025.110033 698 Miller, A. D., Umina, P. A., Weeks, A. R., & Hoffmann, A. A. (2012). Population genetics of the 699 wheat curl mite (Aceria tosichella Keifer) in Australia: Implications for the management 700 of wheat pathogens. Bulletin of Entomological Research , 102(2), 199 –212. 701 https://doi.org/10.1017/S0007485311000526 702 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 21 Olazcuaga, L., Baltenweck, R., Leménager, N., Maia-Grondard, A., Claudel, P., Hugueney, P., & 703 Foucaud, J. (2023). Metabolic consequences of various fruit -based diets in a generalist 704 insect species. eLife, 12, e84370. https://doi.org/10.7554/eLife.84370 705 Paredes, D., Rosenheim, J. A., Chaplin-Kramer, R., Winter, S., & Karp, D. S. (2021). Landscape 706 simplification increases vineyard pest outbreaks and insecticide use. Ecology Letters , 707 24(1), 73–83. https://doi.org/10.1111/ele.13622 708 Parrett, J. M., Chmielewski, S., Aydogdu, E., Łukasiewicz, A., Rombauts, S., Szubert-Kruszyńska, 709 A., ... & Radwan, J. (2022). Genomic evidence that a sexually selected trait captures 710 genome-wide variation and facilitates the purging of genetic load. Nature Ecology & 711 Evolution, 6(9), 1330-1342. https://doi.org/10.1038/s41559-022-01816-w 712 Poveda, K., Karp, D. S., Chaplin-Kramer, R., Centrella, M., Luttermoser, T., Perez-Alvarez, R., … 713 Grab, H. (2025). The importance of landscape composition for pest control and crop 714 yield: A global quantitative synthesis. Ecology Letters , 28(11), e70250. 715 https://doi.org/10.1111/ele.70250 716 Remold, S. (2012). Understanding specialism when the jack of all trades can be the master of 717 all. Proceedings of the Royal Society B: Biological Sciences , 279(1749), 4861 –4869. 718 https://doi.org/10.1098/rspb.2012.1990 719 Saarinen, K., Laakso, J., Lindström, L., & Ketola, T. (2018). Adaptation to fluctuations in 720 temperature by nine species of bacteria. Ecology and Evolution , 8(5), 2901 –2910. 721 https://doi.org/10.1002/ece3.3823 722 Sant, D. G., Woods, L. C., Barr, J. J., & McDonald, M. J. (2021). Host diversity slows 723 bacteriophage adaptation by selecting generalists over specialists. Nature Ecology & 724 Evolution, 5(3), 350–359. https://doi.org/10.1038/s41559-020-01364-1 725 Sexton, J. P., Montiel, J., Shay, J. E., Stephens, M. R., & Slatyer, R. A. (2017). Evolution of 726 ecological niche breadth. Annual Review of Ecology, Evolution, and Systematics , 48(1), 727 183–206. https://doi.org/10.1146/annurev-ecolsys-110316-023003 728 Shih, P. -Y., Sugio , A., & Simon, J. -C. (2023). Molecular mechanisms underlying host plant 729 specificity in aphids. Annual Review of Entomology , 68(1), 431 –450. 730 https://doi.org/10.1146/annurev-ento-120220-020526 731 Skoracka, A., Laska, A., Radwan, J., Konczal, M., Lewandowski, M., Puchalska, E., … Kuczyński, 732 L. (2022). Effective specialist or jack of all trades? Experimental evolution of a crop pest 733 in fluctuating and stable environments. Evolutionary Applications, 15(10), 1639–1652. 734 https://doi.org/10.1111/eva.13360 735 Skoracka, A., Lopes, L. F., Alves, M. J., Miller, A., Lewandowski, M., Szydło, W., … Kuczyński, L. 736 (2018a). Genetics of lineage diversification and the evolution of host usage in the 737 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 22 economically important wheat curl mite, Aceria tosichella Keifer, 1969. BMC 738 Evolutionary Biology, 18(1), 122. https://doi.org/10.1111/bij.12213 739 Skoracka, A., Rector, B. G., & Hein, G. L. (2018b). The interface between wheat and the wheat 740 curl mite, Aceria tosichella , the primary vector of globally important viral diseases. 741 Frontiers in Plant Science, 9, 1098. https://doi.org/10.3389/fpls.2018.01098 742 Subramanian, M., Jayakumar, S., Richhariya, S., & Hasan, G. (2013). Loss of IP3 receptor 743 function in neuropeptide secreting neurons leads to obesity in adult Drosophila. BMC 744 Neuroscience, 14(1), 157. https://doi.org/10.1186/1471-2202-14-157 745 Thorvaldsdóttir, H., Robinson, J. T., & Mesirov, J. P. (2013). Integrative Genomics Viewer (IGV): 746 High-performance genomics data visualization and exploration. Briefings in 747 Bioinformatics, 14(2), 178–192. https://doi.org/10.1093/bib/bbs017 748 VanWallendael, A., Soltani, A., Emery, N. C., Peixoto, M. M., Olsen, J., & Lowry, D. B. (2019). A 749 molecular view of plant local adaptation: incorporating stress -response networks. 750 Annual Review of Plant Biology , 70(1), 559 –583. https://doi.org/10.1146/annurev-751 arplant-050718-100114 752 Van Leeuwen, T., & Dermauw, W. (2016). The molecular evolution of xenobiotic metabolism 753 and resistance in chelicerate mites. Annual review of entomology , 61(1), 475 -498. 754 https://doi.org/10.1146/annurev-ento-010715-023907 755 Veenstra, J. A., Rombauts, S., & Grbić, M. (2012). In silico cloning of genes encoding 756 neuropeptides, neurohormones and their putative G -protein coupled receptors in a 757 spider mite. Insect Biochemistry and Molecular Biology , 42(4), 277 –295. 758 https://doi.org/10.1016/j.ibmb.2011.12.009 759 Yang, H., Zhang, R., Zhang, Y., Liu, Q., Li, Y., Gong, J., & Hou, Y. (2020). Cathepsin -L is involved 760 in degradation of fat body and programmed cell death in Bombyx mori . Gene, 760, 761 144998. https://doi.org/10.1016/j.gene.2020.144998 762 763 Data Accessibility and Benefit-Sharing Section 764 The data supporting the findings of this study, as well as the code used to generate the results, 765 are openly available. All raw sequencing data have been deposited in the NCBI BioProject 766 under accession number PRJNA1395124. The assembled reference genome and its annotation 767 are available via the ORCAE platform 768 (https://bioinformatics.psb.ugent.be/orcae/overview/Aceto). Scripts used to analy se the 769 evolve-and-resequence experiment are available at 770 https://github.com/konczal/WCM_EvolveResequence. Benefits from this research accrue 771 from the sharing of our data and results as described above. 772 773 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 23 Author Contributions 774 ALM, MK, AS conceptualized and designed the study; ALM, AS, ML, JRaubic performed 775 research; ALM under supervision of SR conducted genomic data analysis, MK performed 776 analysis of evolve-and-resequence data; ALM, MK, JRadwan, AS, SR critically interpreted the 777 results; MK, AS, JRadwan wrote the manuscript with assistance of ALM, ML, SR; All authors 778 contributed substantially to revisions, read and approved the final version of the manuscript. 779 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 24 Tables 780 Table 1. GO terms enriched among genes differentiated between treatments. 781 GO.ID Term Annotated Significant Expected p-value aTH vs cH Biological Processes GO:0007156 homophilic cell adhesion via plasma membrane adhesion molecules 8 2 0.18 0.012 GO:0016925 protein sumoylation 1 1 0.02 0.022 GO:0006809 nitric oxide biosynthetic process 1 1 0.02 0.022 Molecular Functions GO:0140359 ABC-type transporter activity 10 3 0.23 0.0012 GO:0005220 inositol 1,4,5-trisphosphate-gated calcium channel activity 1 1 0.02 0.0234 GO:1990817 poly(A) RNA polymerase activity 1 1 0.02 0.0234 GO:0070679 inositol 1,4,5 trisphosphate binding 1 1 0.02 0.0234 GO:0004517 nitric-oxide synthase activity 1 1 0.02 0.0234 GO:0019948 SUMO activating enzyme activity 1 1 0.02 0.0234 GO:0016500 protein-hormone receptor activity 2 1 0.05 0.0462 aTH vs cT Biological Processes GO:0006302 double-strand break repair 1 1 0.03 0.025 GO:0006890 retrograde vesicle-mediated transport, Golgi to endoplasmic reticulum 1 1 0.03 0.025 GO:0006351 DNA-templated transcription 75 3 1.88 0.042 GO:0006814 sodium ion transport 2 1 0.05 0.049 GO:0042908 xenobiotic transport 2 1 0.05 0.049 Molecular Functions GO:0005272 sodium channel activity 1 1 0.02 0.023 GO:0005220 inositol 1,4,5-trisphosphate-gated calcium channel activity 1 1 0.02 0.023 GO:0015016 [heparan sulfate]-glucosamine N- sulfotransferase activity 1 1 0.02 0.023 GO:0070679 inositol 1,4,5 trisphosphate binding 1 1 0.02 0.023 GO:0003906 DNA-(apurinic or apyrimidinic site) endonuclease activity 1 1 0.02 0.023 GO:0003924 GTPase activity 14 2 0.33 0.040 GO:0008559 ABC-type xenobiotic transporter activity 2 1 0.05 0.046 Cellular Components GO:0030126 COPI vesicle coat 2 1 0.04 0.039 cT vs cH Biological Processes GO:0006805 xenobiotic metabolic process 1 1 0.01 0.010 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 25 GO:0042908 xenobiotic transport 2 1 0.02 0.021 GO:0007269 neurotransmitter secretion 2 1 0.02 0.021 Molecular Functions GO:0070679 inositol 1,4,5 trisphosphate binding 1 1 0.02 0.017 GO:0005220 inositol 1,4,5-trisphosphate-gated calcium channel activity 1 1 0.02 0.017 GO:0008559 ABC-type xenobiotic transporter activity 2 1 0.03 0.034 GO:0016409 palmitoyltransferase activity 2 1 0.03 0.034 Cellular Components GO:0034751 aryl hydrocarbon receptor complex 1 1 0.02 0.016 GO:0008021 synaptic vesicle 1 1 0.02 0.016 782 Table 2. Genes identified in genomic region differentiated between cT and cH/aTH lines. 783 Gene ID Location Functional description GO terms aceto37g05560 contig_37: 2869366– 2873696 SAGA associated factor 29-like SAGA Complex (GO:0000124); Protein Binding (GO:0005515) aceto37g05570 contig_37: 2874121– 2875947 guanine nucleotide- binding protein-like 1 GTPase Activity (GO:0003924); GTP Binding (GO:0005525) aceto37g05580 contig_37: 2876685– 2879654 nicotinic acetylcholine receptor subunit alpha 1 Transmembrane Receptor Activity (GO:0004888); Ion Channel Activity (GO:0005216); Ion Channel Activity (GO:0005230); Ion Transport (GO:0006811); Integral Component of Membrane (GO:0016021); Ion Transporter Activity (GO:0022848); Voltage-Gated Ion Channel Ac tivity (GO:0034220); Postsynaptic Membrane (GO:0045211) aceto37g05590 contig_37: 2879959– 2880165 Small EDRK-rich factor 2 - aceto37g05600 contig_37: 2881684– 2888145 inositol 1,4,5- trisphosphate receptor type 1 Inositol 1,4,5 -trisphosphate- sensitive calcium -release channel activity (GO:0005220); Calcium channel activity (GO:0005262); Endoplasmic reticulum (GO:0005783); Calcium ion transport (GO:0006816); Membrane (GO:0016020); Calcium ion transmembrane transport (GO:0070588); Inositol 1,4,5 - trisphosphate binding (GO:0070679) .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 26 aceto37g05610 contig_37: 2889333– 2898267 merozoite surface antigen 2 cell wall protein - aceto37g05620 contig_37: 2899617– 2900684 Cathepsin L Protease activity (GO:0006508); Metalloendopeptidase activity (GO:0008234) aceto37g05630 contig_37: 2902367– 2904397 1,4-alpha-glucan- branching enzyme Catalytic activity (GO:0003824); Alcohol dehydrogenase (GO:0003844); Glycosidase activity (GO:0004553); Carbohydrate metabolic process (GO:0005975); Carbohydrate binding (GO:0005978); Polysaccharide binding (GO:0043169) aceto37g05640 contig_37: 2908779– 2910386 coronin-6 Protein Binding (GO:0005515) aceto37g05650 contig_37: 2939191– 2941792 AllostatinC receptor type2-like G-protein coupled receptor activity (GO:0004930); G -protein coupled receptor signaling pathway (GO:0007186); Integral component of membrane (GO:0016021) aceto37g05660 contig_37: 2942610– 2948825 glycine-rich cell wall structural protein 1.8-like - aceto37g05670 contig_37: 2956698– 2959139 helix-loop-helix protein 6- like Protein Binding, Specifically to Protein Domains (GO:0046983) aceto37g05680 contig_37: 2959826– 2960791 dehydrogenase/reductase SDR family member 7-like Oxidoreductase Activity (GO:0016491) 784 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 27 Figures 785 786 787 Figure 1. Identification of SNPs differentiated between evolved populations. 788 A–C: Population growth rates of Aceria tosichella populations experimentally evolved under 789 different environmental conditions: aTH – alternating wheat and barley, cH – constant barley, 790 and cT – constant wheat. Growth rates were measured on wheat (A), barley (B), and smooth 791 brome (C) (adapted from Skoracka et al. 2022). D: Number of SNPs identified as significantly 792 differentiated among the aTH, cH, and cT line ages. E–G: Manhattan plots showing SNP 793 differentiation between evolved populations. The horizontal line indicates the threshold for 794 significantly differentiated SNPs. Violet points represent SNPs differentiated in all three 795 pairwise comparisons; red points indicate SNPs differentiated in comparisons involving cT vs. 796 aTH/cH. 797 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint 28 798 Fig. 2. Characterization of the genomic region with the highest differentiation in 799 experimentally evolved populations. A: Changes in major allele frequencies for all SNPs 800 located between positions 2,870,000 and 2,990,000 bp on contig_37 (the focal genomic 801 region). B: Genome-wide distribution of nucleotide diversity in the stock population. The 802 vertical red line indicates nucleotide diversity in the focal genomic region. C: Genome-wide 803 distribution of Tajima’s D values in the stock population. The vertical red line indicates the 804 Tajima’s D value for the focal region. D–F: Patterns of allele frequency change for three SNPs 805 located within the focal region. D–E: SNPs showing differentiation in all three pairwise 806 comparisons. F: SNP located within the ITPR1 gene. G–I: Zoomed-in views of the regions of 807 interest on Manhattan plots. 808 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted January 8, 2026. ; https://doi.org/10.64898/2026.01.08.698337doi: bioRxiv preprint