High-resolution genetic mapping reveals cis-regulatory and copy number variation in loci associated with cytochrome P450-mediated detoxification in a generalist arthropod pest

Chemical control strategies are driving the evolution of pesticide resistance in pest populations. Understanding the genetic mechanisms of these evolutionary processes is of crucial importance to develop sustainable resistance management strategies. The acaricide pyflubumide is one of the most recently developed mitochondrial complex II inhibitors with a new mode of action that specifically targets spider mite pests. In this study, we characterize the molecular basis of pyflubumide resistance in a highly resistant population of the spider mite Tetranychus urticae. Classical genetic crosses indicated that pyflubumide resistance was incompletely recessive and controlled by more than one gene. To identify resistance loci, we crossed the resistant population to a highly susceptible T. urticae inbred strain and propagated resulting populations with and without pyflubumide exposure for multiple generations in an experimental evolution set-up. High-resolution genetic mapping by a bulked segregant analysis approach led to the identification of three quantitative trait loci (QTL) linked to pyflubumide resistance. Two QTLs were found on the first chromosome and centered on the cytochrome P450 CYP392A16 and a cluster of CYP392E6-8 genes. Comparative transcriptomics revealed a consistent overexpression of CYP392A16 and CYP392E8 in the experimental populations that were selected for pyflubumide resistance. We further corroborated the involvement of CYP392A16 in resistance by in vitro functional expression and metabolism studies. Collectively, these experiments uncovered that CYP392A16 N-demethylates the toxic carboxamide form of pyflubumide to a non-toxic compound. A third QTL coincided with cytochrome P450 reductase (CPR), a vital component of cytochrome P450 metabolism. We show here that the resistant population harbors three gene copies of CPR and that this copy number variation is associated with higher mRNA abundance. Together, we provide evidence for detoxification of pyflubumide by cytochrome P450s that is likely synergized by gene amplification of CPR.

action [34]. As revealed by recent studies, T. urticae is highly amenable to BSA genetic mapping, with

113
The development of pesticides with novel modes of action is crucial to maintain control over 114 pest populations that have developed multi-resistance. Pyflubumide is a carboxanilide acaricide and 115 inhibits mitochondrial complex II, or succinate dehydrogenase. This recently developed compound is 116 highly selective to spider mite pest species [39,40] and acts as a pro-acaricide that requires bioactivation 117 into a toxic derivative. Pyflubumide is converted within the spider mite body into its active deacylated 118 metabolite (carboxamide form, NNI-0711-NH) that strongly inhibits the mitochondrial complex II 119 through binding to the quinone-binding pocket [39,40]. In the current study, we characterize the genetic 120 basis of resistance to this novel acaricide using two resistant T. urticae strains (JPR-R1 and JPR-R2).

121
Target-site resistance is not involved in high-level pyflubumide resistance in these strains, and 122 synergism/antagonism bioassays strongly suggested the involvement of cytochrome P450s [41]. In this 123 study, we employed high-resolution BSA genetic mapping and uncovered three QTL associated with 124 pyflubumide resistance. In parallel, transcriptomic analyses were conducted on these experimental 125 populations to further characterize the molecular mechanisms of pyflubumide resistance. In vitro 126 functional characterization of the cytochrome P450 CYP392A16, a QTL candidate, suggested a role for 127 resistance to pyflubumide and its active (deacylated) carboxamide metabolite. Additionally, another 128 QTL interval was centered on cytochrome P450 reductase (CPR), which is required for cytochrome 129 P450 activity, and is amplified in the resistant populations. Together, our results indicate that the 130 mechanism of pyflubumide detoxification involves cytochrome P450 activity enhanced by gene 131 amplification at the CPR locus.

135
Mode of inheritance of pyflubumide resistance 136 Previously, we selected a field-collected strain for higher resistance levels to pyflubumide and 137 obtained two highly resistant strains, JPR-R1 and JPR-R2 [41]. In this study, we first confirmed the high 138 levels of pyflubumide resistance and uncovered LC 50 values higher than 1,000 mg/L of pyflubumide for 139 both strains (Table S1). Reciprocal crosses between resistant (JPR-R1 and JPR-R2) and susceptible (JP-140 S and Wasatch) strains revealed that pyflubumide resistance is not maternally inherited and has an 141 incomplete recessive mode of inheritance (Fig 1 and Table S1). The dose-response relationships of the 142 two back-crossed populations were significantly different from the expected curves for a recessive 143 monogenic mode of inheritance (p< 0.001) [42]. These results indicate that pyflubumide resistance in 144 both JPR-R1 and JPR-R2 was determined by multiple loci.

145
Response to pyflubumide selection pressure in an experimental evolution set-up

146
To further characterize the genetic architecture of pyflubumide resistance, we crossed a single 147 haploid JPR-R1 male to diploid females of the susceptible inbred WasX strain. Mortality at 10, 50, 100,

148
and 500 mg/L of pyflubumide was significantly different between JPR-R1 and WasX (X 2 = 41.631, df= 149 1, p < 0.0001) (Fig 2). The segregating population resulting from the WasX x JPR-R1 cross was allowed 150 to expand and develop for two months (approximately three-four generations), after which 11 pairs (22 151 populations) were set-up and allowed to propagate for approximately 35 generations. In this pairwise 152 set-up, one sister population was exposed to pyflubumide on sprayed plants (a selected population),

153
whereas the other was maintained on unsprayed plants (an unselected, control population). After a 154 prolonged selection with increasing doses of pyflubumide, all populations were phenotyped, revealing 155 that mortality at four pyflubumide doses was significantly lower for pyflubumide-selected populations 156 as compared to the control populations (Fig 2)

171
Genetic mapping uncovers multiple QTL associated with pyflubumide selection

172
To locate the genomic regions that underlie pyflubumide resistance in the pyflubumide-selected 173 populations, we employed a bulked segregant analysis (BSA) approach using high-quality SNP loci 174 [25,29]. Genome-wide JPR-R1 allele frequencies between pyflubumide-selected and susceptible for QTL detection with replicated paired selected and control populations [29], and with an FDR of 0.05,

178
we identified two QTL on chromosome 1 (hereafter QTL-1 and QTL-2) and one on chromosome 2

179
(hereafter QTL-3) in the pyflubumide-selected populations ( Fig 4C and Fig S1). In all cases, alleles 180 from the resistant JPR-R1 parent were selected, and for all pyflubumide-selected populations, near-181 fixation of the JPR-R1 haplotype was observed at the peaks of the three QTL regions (Fig S1).

182
Population R6 exhibited the lowest JPR-R1 allele frequency near the peaks of QTL-2 and QTL-3 (Fig   183  S1), a finding that is consistent with its lower resistance levels (Fig 2). We also noted two additional 184 genomic regions where parental JPR-R1 allele frequencies were elevated in a number, but not all, of the 185 replicated pyflubumide-selected populations (hereafter referred to as subsidiary peaks 1 and 2) (Fig 4C).

186
Finally, a BSA mapping approach was also performed using informative SNP loci identified via the RNAseq data. This analysis further confirmed the three QTL, as well as the two subsidiary peaks ( [25,27,33]. Therefore, we examined ~75 kb genomic intervals that bracket the 199 average of the BSA peaks to detect potential causal candidate genes (Fig 5 and Supplementary Data S2).

200
For QTL-1, the cytochrome P450 CYP392A16 was the only gene in the interval with a predicted function 201 that could be readily linked with xenobiotic metabolism of arthropods. CYP392A16 was within 10 kb of 202 the averaged BSA peak (Fig 5A), and was one of three genes within the ~75 kb genomic interval with parental strains had the same copy number of CYP392A16 (a single copy), ruling out copy number 208 variation (CNV) as underpinning the higher transcript levels of CYP392A16 ( Fig S4).

209
A cluster of three cytochrome P450s, CYP392E6-8, was located very close to the BSA peak of 210 QTL-2 ( Fig 5B), and no other genes with obvious roles in arthropod detoxification were located in the 211 interval. CYP392E8, which is located only 2 kb from the QTL-2 peak, was the only gene in the region 212 that exhibited a significant increase in transcription in pyflubumide-selected populations (a fold change 213 of ~2) (Fig 5 and Supplementary Data S2). Although substantial CNV for this cytochrome P450 gene cluster has been reported in other T. urticae strains [29], all CYP genes in the interval were single copy 215 as revealed by genomic read coverage in the parental and segregating populations ( Fig S4).

216
No cytochrome P450s or genes encoding other detoxification enzymes were apparent at the 217 QTL-3 peak region. However, NADPH cytochrome P450 reductase (tetur18g03390, CPR), the 218 necessary electron donor for cytochrome P450 metabolism, was within 20 kb of the averaged BSA peak 219 for QTL-3 ( Fig 5). Previous work with T. urticae revealed an association between a CPR coding 220 sequence variant and cytochrome P450-mediated detoxification of pesticides [29,33]; however, short-  (Fig S3). This ratio was mirrored in the RNAseq read data where up-regulation of 236 CPR transcription was ~3.66 in pyflubumide-selected versus control populations (Fig 5 and Fig S3).

237
Collectively, this indicates that CNV of CPR resulted in higher mRNA levels in JPR-R1 and 238 pyflubumide-selected populations. It was also striking that in the control (unselected) populations, the 239 JPR-R1 allele frequencies dipped to near zero coincident with the CPR-associated QTL-3 peak ( Fig 4A   240 and Fig S1), potentially reflecting a negative pleiotropic effect on fitness that becomes apparent in absence of pyflubumide selection (see Discussion). This pattern was absent or less apparent for the other 242 QTL peaks (for instance, a superficially similar pattern for QTL-1 is likely explained by linkage drag 243 from a proximal locus at ~10 Mb on chromosome 1 that is unrelated to pyflubumide selection).

244
Consistent with previous work [41], our BSA genetic mapping with JPR-R1 excluded target-245 site insensitivity as a contributing mechanism to pyflubumide resistance. Only one of the genes encoding 246 a mitochondrial complex II subunit, the known molecular target(s) of pyflubumide, was near a QTL 247 peak (tetur08g03210, which is located ~640 kb from the BSA peak for QTL-2). Although the two 248 subsidiary BSA peaks did not reach the significance threshold for QTL detection (Fig 4C), we 249 nonetheless analyzed their genic content. For subsidiary peak 1, no obvious candidate detoxification 250 genes were found within the peak ~75 kb region, nor were any genes differentially expressed. For 251 subsidiary peak 2, two CYP genes, CYP392B2 and a predicted CYP pseudogene (tetur02g06640) in the 252 London reference strain, were within the peak ~75 kb region, with the latter exhibiting lower RNA 253 abundance in the pyflubumide-selected populations (a fold change of -2) (Supplementary Data S2).

254
Several non-synonymous single-nucleotide differences were observed in both CYP genes between JPR-255 R1 and WasX. Despite the lack of statistical significance for subsidiary peak 2, CYP392B2 is of interest 256 as it was previously associated with pyflubumide resistance [43] (and see Discussion).

258
To test if CYP392A16 contributes to pyflubumide resistance in JPR-R1, we functionally 259 expressed the gene to determine if the recombinant protein was active against pyflubumide and/or its pointed to N-demethylation as the likely mechanism of the CYP392A16 reaction (Fig 7) [44].

276
Demethylation is expected to result in a mass difference of -14 Da between the NNI-0711-NH and 277 reaction metabolite product ions. This was reflected in almost all substrate and reaction product ions 278 formed by collision-induced dissociation (Fig 7). More specifically, the product ion MS data indicated 279 that demethylation involves a pyrazole methyl group. We arrived at this conclusion because pyrazole-280 containing product ions originating from the reaction product are 14 Da lower than their corresponding 281 substrate product ions of NNI-0711-NH (Fig 7).

283
Pesticide resistance is a rapidly evolving trait within the arthropod phylum [19,45]. Despite its 284 economic importance and its potential to contribute to our understanding of general evolutionary 285 processes, several questions about the genetic architectures underlying pesticide resistance are still 286 outstanding, especially in the case of polygenic resistance. Pyflubumide is a recently developed 287 carboxanilide acaricide that is highly effective against plant-feeding spider mites [39][40][41]46]. Based on 288 synergism/antagonism assays of previous work [41], the resistant T. urticae JPR-R1 and JPR-R2 strains 289 most likely rely on cytochrome P450-mediated detoxification to overcome pyflubumide toxicity. The 290 two strains also exhibit highly similar constitutive and plastic transcriptional patterns on a genome-wide 291 level [41]. Here, we uncovered that pyflubumide resistance of both JPR-R1 and JPR-R2 is determined 292 by multiple loci. Considering the two strains share a common ancestral genetic background with already 293 moderate levels of pyflubmide resistance [41], these findings further suggest that JPR-R1 and JPR-R2 294 likely share (some) factors underlying pyflubumide resistance. Using JPR-R1 and susceptible WasX as parents, BSA genetic mapping revealed three major QTL that are associated with pyflubumide resistance 296 based on the segregating genomes of 11 pyflubumide-selected and control population pairs.

297
Interestingly, a BSA genetic mapping approach based on RNAseq data (instead of genomic DNA) 298 uncovered highly similar allele frequency differences, mirroring the QTL pattern obtained with DNA.

299
Less biological tissue is needed to obtain sufficient RNA, as compared to DNA, for high-quality 300 sequencing, and our findings highlight the feasibility of using RNAseq data for BSA genetic mapping, 301 a consideration that is relevant for microarthropods for which biological material is often limiting.

302
CYP392A16 was located at the center of QTL-1, was over-expressed in pyflubumide-selected 303 populations and the JPR-R1 parent, and metabolized the toxic pyflubumide carboxamide in vitro.

304
Incubation assays with membrane fractions of recombinant protein further uncovered that CYP392A16 305 demethylates the active carboxamide into a derivative that is not toxic to spider mites [46], confirming 306 that CYP392A16 catalyzes a detoxification reaction. The CYP gene cluster of CYP392E6-8 was close 307 to the peak of QTL-2 and CYP392E8 exhibited higher transcription in pyflubumide-selected populations 308 and JPR-R1. In an earlier study, BSA genetic mapping uncovered a broad QTL (referred to as spiro-

309
QTL 2) that included the CYP392E6-8 gene cluster when replicated populations were selected for 310 resistance to spirodiclofen, a pesticide that targets lipid synthesis [29]. This raises the question of  [37,47]. Using the annotated genome of the London reference strain, QTL-2 is predicted 319 to additionally hold several CYP gene fragments and pseudogenes [29,48]. As the CYP multi-gene 320 family is characterized by highly dynamic birth and death rates [49,50], it is possible that these genes 321 could be functional in JPR-R1. Further work is needed to fully characterize the functionality and enzymatic abilities of the multiple cytochrome P450s within QTL-2 and to understand why this genomic 323 region has been a locus of selection for both spirodiclofen and pyflubumide resistance.

324
The causal genetic variants underpinning the over-expression of detoxification genes in resistant  [43].

378
These CYP genes were not differentially expressed between the pyflubumide-selected and control 379 populations and were not included in any of the three well-supported QTL of this study. However, 380 subsidiary peak 2 of the current study largely coincided with this QTL (its peak region also included pesticides [19].

388
In summary, our work uncovered a complex genetic basis for pyflubumide resistance in the 389 arthropod pest T. urticae that likely involves multiple genes with roles in cytochrome P450-mediated 390 metabolism. Our in vitro findings suggest that CYP392A16 detoxifies the active form of the proacaricide 391 pyflubumide, and that cis-regulatory variation at this gene contributes prominently to pyflubumide 392 resistance. We also uncovered a potential novel molecular mechanism of toxicokinetic resistance -gene 393 amplification of CPR -that may act in a synergistic manner.

397
JPR-R1, and JPR-R2. JPR-R1 and JPR-R2 were selected for high levels of pyflubumide resistance from 398 the JP-R strain [41]. WasX was derived from the highly inbred Wasatch strain by two sequential rounds  the RNAseq read data was performed using the regularized log transformed counts.

484
BSA genetic mapping was performed using methods adapted from previous studies [28,29,33], 485 and specifically, using the defaults of the haplodiploid option of Kurlovs et al. [25]. Briefly, after quality 486 control, JPR-R1 was specified as the haplodiploid male parent, and which alleles got passed on from its 498 urticae genome [29]. Finally, the allele frequency differences were averaged across the replicated pairs.

515
The final estimate of CPR copy number was obtained with efficiency (corrected Ct -Ct) .

517
We adopted a previously successful functional expression strategy for CYP392A16 [80].

518
Briefly, competent E. coli BL21STAR cells were co-transformed with the pCW_ompA-CYP392A16 519 and pACYQ_TuCPR plasmids [80]. Transformed cells were grown at 37 °C in Terrific Broth under 520 ampicillin and chloramphenicol selection until the optical density at 595 nm was higher than 0.95 cm -1 .

521
The heme precursor δ-aminolevulinic acid was added to a final concentration of 1 mM. Induction was aliquots at -80 °C. P450 content was measured by CO-difference spectra in reduced CYP392A16 527 membrane samples [81]. CPR activity was estimated by measurements of NADPH-dependent reduction 528 of cytochrome c at 550 nm [82]. CYP activity was tested using the chemi-luminescent substrate 529 Luciferin-ME EGE [80]. Total protein concentrations of CYP392A16 membranes were determined 530 using the Bradford assay [83].  of control (S1-S11) and pyflubumide-selected (R1-R11) populations (depicted in blue and orange, 595 respectively). Error bars represent standard errors. between the eleven paired populations (R1-S1 to R11-S11), as assessed in a sliding window analysis.