Systemic α-synuclein injection triggers selective neuronal pathology as seen in patients with Parkinson’s disease

Parkinson’s disease (PD) is an α-synucleinopathy characterized by the progressive loss of specific neuronal populations. Here, we develop a novel approach to transvascularly deliver proteins of complex quaternary structures, including α-synuclein preformed fibrils (pff). We show that a single systemic administration of α-synuclein pff triggers pathological transformation of endogenous α-synuclein in non-transgenic rats, which leads to neurodegeneration in discrete brain regions. Specifically, pff-exposed animals displayed a progressive deterioration in gastrointestinal and olfactory functions, which corresponded with the presence of cellular pathology in the central and enteric nervous systems. The α-synuclein pathology generated was both time dependent and region specific. Interestingly, the most significant neuropathological changes were observed in those brain regions affected in the early stages of PD. Our data therefore demonstrate for the first time that a single, transvascular administration of α-synuclein pff can lead to selective regional neuropathology resembling the premotor stage of idiopathic PD. Furthermore, this novel delivery approach could also be used to deliver a range of other pathogenic, as well as therapeutic, protein cargos transvascularly to the brain.

Introduction
Parkinson's disease (PD) is a progressive neurological disorder defined by the loss of nigrostriatal dopaminergic neurons and formation of α-synuclein-positive Lewy bodies (LB) [1]. It is well established that neurodegeneration and LB pathology affects only a subset of vulnerable cells in PD. Recent evidence, however, suggests that pathogenic α-synuclein can propagate between cells and seed pathology in a prion-like manner [2,3]. The question therefore arises as to why such a spread only generates pathology in certain neuron populations. Indeed, while neuropathological studies suggest that there is a stereotypic, spatiotemporal pattern of LB deposition in post-mortem PD brains [4,5], studies using a direct intracerebral injection of α-synuclein preformed fibrils (pffs) to transform endogenous α-synuclein into pathological species are dependent on the site of injection [6][7][8][9]. This approach therefore does not capture the pathology of PD nor does it reflect the experimental evidence suggesting that α-synuclein pathology can spread from the periphery, via the enteric (ENS), to the central nervous system (CNS) [10].
Misfolded α-synuclein has been shown to cross the blood-brain barrier following repeated systemic administration, possibly through downregulating endothelial tightjunction proteins, but at a level that does not generate α-synuclein pathology or neurodegeneration [11,12]. We therefore sought to optimize the efficiency of this transvascular protein transduction by taking advantage of a modified rabies virus glycoprotein (RVG9R), which has previously been used to mediate transvascular RNA delivery to the brain [13,14]. Negatively charged nucleotides can be electrostatically complexed with a cationic nonaarginine sequence, while the RVG ligand itself interacts with a neural-specific acetylcholine receptor to induce a transient membrane inversion and enable translocation of the RNA cargo into the cytoplasm [15]. By interacting with anionic amino acids on the surface of proteins, we hypothesized that RVG9R could also be used to enable efficient protein transduction into the CNS following peripheral administration, without it affecting the physiological activities of the cargo protein.
Using this approach, we demonstrate for the first time that systemic delivery of pathological α-synuclein pffs using this method can lead to selective pathology, neurodegeneration, and functional deficits in the rat CNS and ENS. Critically, this occurs in a pattern that strikingly mirrors that reported in early-stage, idiopathic PD patients. This model therefore recapitulates some key features of early-stage PD, and so has significant potential for preclinical testing of putative disease-modifying agents, when such agents are thought to have their maximal chance of therapeutic benefit. Finally, the strategy we used for transvascular delivery of α-synuclein could also be leveraged to deliver other pathogenic and/or therapeutic proteins to the CNS.

Peptide/protein complex formation
To prepare the complex for delivery, peptides and protein cargos were mixed in PBS in a glass vial (Supelco) and left on an orbital shaker (Stuart SSM1, 60 rpm), for 30 min at room temperature before use. The optimal (peptide:cargo protein) molar ratio for application was empirically determined to be 10:1 for monomeric proteins. The molecular weights for the heterodimeric chondroitinase ABC (ChABC) and misfolded α-synuclein pff were assumed to be 120 and 14.5 kDa, respectively, and the peptide/protein complex was mixed at a 10:1 molar ratio.

Animal procedures
Female Sprague-Dawley (SD) rats, weighing around 200 g at the time of injection, were housed with unrestricted access to food and water unless otherwise stated. Animals were housed in groups of four on a 12 hours light/dark cycle. For intravenous delivery, a total volume of 300 μl of the peptide/protein complex, containing 100 μg of the cargo proteins or 3 U ChABC per animal, was administered via the lateral tail vein.
A detailed "Material and methods" section can be found in the Supplementary information.

The RVG9R system mediates efficient protein transduction
To ascertain if the RVG9R system could mediate protein transduction in vitro, monomeric human α-synuclein or the green fluorescent protein (GFP) were added in culture alone or mixed with RVG9R at 1:10 molar ratio. Overnight incubation of RVG9R:α-synuclein or RVG9R:GFP with primary rat or human foetal cortical neurons, respectively, induced cellular update of the cargo proteins (Supplementary Fig. 1a). Adult SD rats were then intravenously administered with either RVG9R:α-synuclein, RVG9R: GFP, or appropriate RVG9R or cargo controls and, after 24 hours, distinctive axonal human-specific α-synuclein and somatic GFP stainings were detectable in various CNS regions ( Supplementary Fig. 1b). To assess the biodistribution of RVG9R-mediated cargo delivery, animals were given a single intravenous RVG9R:GFP administration and, after 24 hours, GFP expression was detectable in the CNS only when complexed with RVG9R ( Supplementary  Fig. 1c). The presence of GFP in the peripheral organs of RVG9R:GFP injected animals could be attributed to the rabies virus-mediated anterograde trans-synaptic from the CNS to the peripheral nerve terminals [16], although further investigation is required to validate whether RVG9R behaves in a similar fashion as the intact virus.
We next assessed whether the RVG9R system can transvascularly deliver biologically active proteins, such as the heterodimeric ChABC enzyme. The formation of chondroitin sulfate proteoglycans-enriched perineuronal nets (PNNs), which are abundantly expressed in the sensorimotor cortex, impedes CNS plasticity, and this can be modulated by enzymatic digestion of chondroitin sulfate chains using ChABC [17]. A single peripheral injection of RVG9R:ChABC treatment progressively reduced Wisteria floribunda agglutinin (WFA)-positive PNNs. Consistently the expression of unsulfated C-0-S stubs, resulting from ChABC digestion, was also elevated, although the level of extracellular cartilage link protein 1 remained unaffected ( Supplementary Fig. 1d, e). This suggested that RVG9Rdelivered ChABC removed chondroitin sulfates without affecting PNN structure.

RVG9R reversibly complexes with the C-terminus of α-synuclein
To further understand how RVG9R complexes with proteins, we first studied the interaction between RVG9R and α-synuclein using analytical ultracentrifugation (AUC). Recombinant human α-synuclein was observed to be monomeric but, in the presence of RVG9R, a shift of the sedimentation coefficient was observed, indicative of complex formation (Fig. 1a, b). The interaction was studied in detail by NMR spectroscopy using 15 N-labelled α-synuclein, for which chemical-shift assignments were available [18]. RVG9R was titrated into 15 N-labelled α-synuclein and the chemical-shift changes monitored by 15 N heteronuclear single quantum coherence (HSQC) spectroscopy, until saturation was achieved (Fig. 1c). By far the most significant chemical-shift changes were observed for residues 106-140 at the C-terminus of α-synuclein, which is both anionic and disordered (Fig. 1d). The shift changes for the peaks showing the largest changes were fitted to a one-site model to calculate affinity (Fig. 1e, f). Equilibrium was achieved within 2 min of the addition of RVG9R for each titration point. An average dissociation constant (K d ) of 80 ± 4 µM for RVG9R:α-synuclein was calculated from the fitted K d for residues 120-140 at the C-terminus, leading to a lower limit of k off > 0.004 s −1 , corresponding to a dissociative half-life (ln2/k off ) of <170 s (Fig. 1f). As k off is independent of concentration, this would be expected to hold true in vivo.
A single intravenous administration of RVG9R:αsynuclein pff generates GI and olfactory deficits It was recently shown that full-length human α-synuclein pffs have a highly ordered core with unstructured tails [19]. We therefore hypothesized that this C-terminus region could be accessible to RVG9R complexing and sought to investigate whether RVG9R-mediated delivery of wild-type human pff ( Supplementary Fig. 2a) can induce α-synuclein pathology when delivered non-selectively to the nervous system. Rats were given a single intravenous administration of either RVG9R:pff, pff alone, or PBS, and followed behaviourally over a 6-month period. Animals receiving RVG9R alone were not included as our previous work has shown that repeated administration of the RVG9R control did not generate any phenotypes [14]. Over the assessment period, there were no significant group differences in their locomotor function, motor coordination, and affective behaviour ( Supplementary Fig. 3). There were also no significant group differences in body weight gain or in faecal profiles over the 6-month period. In contrast, solid GI transit was significantly prolonged at 6 months post injection, specifically in the RVG9R:pff cohort, suggestive of subtle GI dysfunction (Fig. 2a). Post-mortem analysis in these animals revealed phosphorylated α-synuclein inclusions and PK-resistant aggregates in the myenteric plexus, but not in the submucosal plexus, of the duodenum (Fig. 2b, Supplementary Fig. 2b). Whilst this was not accompanied by overt cell loss or changes in total α-synuclein expression in the duodenum (which exists in dimeric forms [20]), there was a significant increase in glial fibrillary acidic protein (GFAP) expression in the LMMP, as well as an elevation in GFAP + reactive astrocytes around the myenteric plexus ( Supplementary Fig. 2c, d).
Together with GI dysfunction, hyposmia is another clinical feature of premotor PD. We therefore assessed the ability of the RVG9R:pff exposed animals to use olfactory cues for foraging after overnight fasting. In the buried food test, the RVG9R:pff animals displayed an increasing retrieval latency over time (Fig. 2c). Similarly, odour-cross habituation was disrupted in these animals at 6 months post injection, with a significant group effect on the discrimination index (Fig. 2d, e). Post-mortem histopathological analysis of the olfactory bulb (OB) revealed substantial phosphorylation of α-synuclein in the mitral and granule, but not in the glomerular cell layer, in RVG9R:pff animals ( Fig. 2f). Using the novel pentameric thiophene derivative (pFTAA) that specifically recognizes oligomeric and fibrillar protein species [21], we observed no evidence of abnormal protein folding in these regions (Supplementary Fig. 2e). There was however a significant elevation in the number of dopaminergic cells residing locally in the periglomerular area (Fig. 2f, g), although the level of tyrosine hydroxylase (TH) expression and the concentration of dopamine in the whole OB remained unchanged (Supplementary Fig. 2f-h).

Structural imaging changes on MRI and their functional correlates
We next examined if there were any macroscale changes in rat brain volume following RVG9R:pff exposure and if so, how these related to the subtle behavioural phenotype of these animals. Comparing the relative volumes of the summary regions of interest (ROIs) for total grey and white matter and the individual rat brain ROIs (n = 79), there were no significant treatment effects after multiple comparison adjustment (Fig. 3a). The relative volumes of the olfactory system were also not affected, despite the presence of olfactory deficits (Fig. 3a). The relative volumes of two individual atlas ROIs in the RVG9R:pff rats, the medial cerebellar peduncle (MCP) and the pyramidal decussation (PDec), were increased (+5) and decreased (−14), respectively (Fig. 3a). Effect sizes (Cohen's d) for each region were +1.3 and −1.3, respectively. Whole-brain, voxel-wise TBM analysis confirmed these data, with no significant volumetric changes between treatment groups after stringent correction for multiple comparisons (family-wise error rate, p < 0.05). There were however clusters of voxels with both apparent changes in volume diffusely in the neocortex, brain stem, and cerebellum, at an exploratory threshold of p < 0.05 (uncorrected, Fig. 3b).
We next explored our MRI dataset for correlations between relative brain ROI volumes in the RVG9R:pff treated animals with their performance in the buried food test. Since this is an exploratory analysis we did not correct the resulting p values for multiple comparisons. This analysis suggested that two ROIs, the neocortex and anterior commissure, were negatively correlated with olfactory impairment, such that smaller relative volumes of these structures were associated with longer retrieval latency (Fig. 3c). A further 2/80 ROIs, the pontine nuclei and the bed nucleus of the stria terminalis (BNST), were positively correlated with performance on the buried food test (Fig. 3c). There were no such correlations with any of these regions in the control group.

RVG9R:pff treatment induces neurodegeneration in the midbrain and brain stem
We next sought to define whether there was any relevant CNS histopathology in brain regions reported to show neuropathology in early-stage PD [4,5], as well as relevant sites affected at trend level in the MRI data. There were significant reductions of TH + and vesicular monoamine  (Fig. 4a, b, Supplementary Fig. 4a, b). These pathologies were likely to have resulted from endogenous rodent α-synuclein aggregation, as no human α-synuclein was detected in the host brain at 6 months after the one-off intravenous administration ( Supplementary Fig. 4c), consistent with the injected pff seeding host pathology. There was also a progressive decline in noradrenergic cell density in the locus ceruleus (LC), cholinergic neurodegeneration in the dorsal motor nucleus of the vagus nerve (DMN-X), as well as in the hypoglossal nerve (DMN-XII) (Fig. 4d-f). These inclusion bodies in the DMN-X were PKresistant, LB-like, amorphous precipitates that were reactive to thioflavin-S ( Supplementary Fig. 5a). Such inclusions were also ubiquitylated and phosphorylated at the serine-129 residue of α-synuclein ( Supplementary Fig. 5b).
RVG9R:pff treatment generates site-specific and time-dependent α-synucleinopathy These data provided justification for further analysis, which revealed a differential pattern of α-synuclein pathology across different brain regions. For instance, LB-like, phosphorylated α-synuclein inclusions were observed in the LC and DMN-X, with pFTAA-positive oligomers found surrounding the inclusion bodies (Fig. 5a). The presence of both voluminous and threadlike Lewy pathologies, similar to those described in the dorsal medulla of idiopathic PD patients [22], could be detected in the DMN-X at 6 months after RVG9R:pff injection (Fig. 5a). The presence of αsynuclein pathology could be observed in other brain regions that have previously been implicated in early-stage PD, such as the anterior olfactory nucleus and the SNc. The selective nature of the pathology observed was also evidenced by the absence of pathology at other sites, such as the thalamus and hippocampus, at 6 months after RVG9R: pff injection (Fig. 5a).
Finally, the extent of α-synuclein aggregation was assessed using western immunoblotting. Triton-, SDS-, and urea-soluble α-synuclein were sequentially extracted from different regions of the CNS and ENS. At 6 months after RVG9R:pff injection, Triton-soluble fractions displayed both total and phosphorylated, monomeric α-synuclein in all brain regions and dimeric α-synuclein in the duodenum (Fig. 5b). The presence of high molecular weight (HMW) α-synuclein in the SDS-soluble fraction, with the absence of α-synuclein in the cerebellum, OB, and in the duodenal mucosa in the urea-soluble fraction, indicates that there is a region-specific susceptibility to intravascularly delivered pathogenic pff (Fig. 5b). Using the cortex as an example, aggregation of α-synuclein in RVG9R:pff rats was also found to be a time-dependent process, as both the size and intensity of HMW α-synuclein was progressively increased in the urea-soluble fraction in animals receiving RVG9R:pff (Fig. 5c, d). The presence of phosphorylated α-synuclein in the Triton-and SDS-soluble fractions of PBS-injected animals was likely associated with the membrane fixation process used to enhance the sensitivity of detecting HMW α-synuclein [23].

Discussion
There is increasing evidence demonstrating that misfolded α-synuclein can spread from one cell to another, possibly through direct synaptic transfer, or via exosomes, tunnelling nanotubes, and LAG3-dependent endocytosis [24][25][26][27]. The relevance of such spread in the pathogenesis of idiopathic PD however remains debated [28]. In this study, we have developed a novel method to systemically deliver pathogenic α-synuclein pff, and demonstrated that a single intravascular injection using this approach can induce a region-selective and time-dependent neuronal pathology resembling that is seen in premotor, idiopathic PD. This strategy therefore provides a potential model system to study the earliest stages of disease pathogenesis in PD. This not only will improve our understanding of the mechanism and functional relevance of α-synuclein pathogenic spread, but will also be beneficial for the evaluation of novel disease-modifying therapeutics. Moreover, therapeutic proteins could also be delivered using this same carrier system, as we have described here using RVG9R:ChABC and previously with RVG9R:p137 in neurotoxic models of PD [14].
While there are many symptomatic treatments for the dopaminergic aspects of PD, no disease-modifying therapies have yet been successfully translated to the clinic. This is in part due to the lack of animal models recapitulating the earliest stages of this disease. Most neurotoxin-based models are acute in nature without generating α-synucleinopathy. Similarly, transgenic models only represent a small proportion of the familial PD cases, and Fig. 2 Transvascular RVG9R:pff delivery generates GI and olfactory impairments. a No group differences were observed in body weight and faecal profiles over time, although there was a treatment effect in GI transit (F 2,55 = 4.559, p = 0.015), reaching statistical significance at 6 months after a single intravenous RVG9R:pff injection (p = 0.021 vs pff alone, p = 0.028 vs PBS). b Abnormally phosphorylated (DABperoxidase, green arrows) and PK-resistant α-synuclein inclusions (red) were detected in the duodenal myenteric plexus (peripherin, green). c There was a significant treatment effect in the latency to retrieve food at 6 months post injection (F 2,54 = 4.483, p = 0.016), post hoc analysis vs pff alone (p = 0.002) and PBS (p = 0.003) groups. d There was a trend for the RVG9R:pff treated rats to be less capable at discriminating the new odour presentation (A3-B1) at 6 months post injection, although this failed to reach significance after post hoc adjustment. e There was a significant reduction on the discrimination index in the RVG9R:pff rats at 6 months post injection (F 4,43 = 2.710, p = 0.042), post hoc analysis vs PBS at 6 months post injection (p = 0.041) and RVG9R:pff at 2 months post injection (p = 0.022). f Extensive α-synuclein phosphorylation (Ser129, red) was found in the mitral and granule cell layers, with more TH + neurons (green) also being observed in the periglomerular region in the RVG9R:pff rats. g There was a significant increase in the TH cell counts per glomerulus in the RVG9R:pff rats at 6 months post injection (F 2,33 = 39.843, p < 0.001), post hoc analysis vs pff only (p < 0.001) and PBS (p < 0.001) groups. *p < 0.05, **p < 0.01, ***p < 0.001 compared with PBS group at the same time point. Data were analyzed using multivariate ANOVA or repeated measure ANOVA in d. n = 6-8 per group at 2/ 4 months post injection, n = 9-12 at 6 months post injection. Three glomeruli were randomly selected for analysis from each animal and averaged in g. Nuclei were counterstained with Hoechst (blue). Scale bars, 30 µm At an exploratory threshold (p < 0.05 uncorrected), RVG9R:pff exposure induced subtle effects on local volume and revealed previously unappreciated clusters of voxels with volumetric changes in the brain stem and cerebellum. c Linear regression revealed significant correlations between the relative volumes of four brain regions (r 2 = 0.50 for neocortex, r 2 = 0.39 for anterior commissure, r 2 = 0.43 for pontine nuclei, and r 2 = 0.36 for BNST) and retrieval latency in the buried food test in RVG9R:pff treated rats. Associated p values (uncorrected) are shown. *p < 0.05, **p < 0.01. Data were analyzed using t-tests for each atlas ROI and then corrected using FDR. Data shown in b are the results of statistical testing for differences in the log scaled Jacobean determinant between RVG9R:pff and PBS-injected rats, at an exploratory threshold of p < 0.05 uncorrected for multiple comparisons. n = 12 per group even in these models, α-synuclein pathology is often not sufficient to induce neurodegeneration at sites relevant to PD [29]. Animal models generated by direct pff inoculation faithfully reproduce many of the pathological features that are seen in the post-mortem brains from idiopathic PD cases [6][7][8][9]. Nevertheless, direct intracerebral pff injection never induces pathology outside the brain, at odds with the clinical picture of idiopathic PD in which α-synuclein pathology is abundant in the ENS and may precede that found centrally [30]. Although a recent study has demonstrated that an enteric inoculation of pff generates pathological features in the GI system that propagates into the CNS [31], such findings have not been universally reported [32][33][34]. Thus, whether ENS abnormalities precede and relate to the onset of motor problems in idiopathic PD is still a matter of debate [35]. Nevertheless, a transvascular delivery of pffs that critically generates pathology only in selective, PDrelated CNS and ENS sites, renders our RVG9R:pff model clinically relevant, given its ability to recapitulate the α-synuclein pathology seen in many patients with earlystage PD.
The molecular mechanism causing such differential susceptibility remains elusive. As the propagation of α-synucleinopathy is dependent on the level of α-synuclein of the recipient cells [6], it has been hypothesized that selective neuronal vulnerability may relate to the differential expressions of endogenous α-synuclein. Indeed, glutamatergic neurons in the hippocampus express higher endogenous α-synuclein levels and are more susceptible to PFF toxicity than GABAergic neurons [36]. Similarly, a subset and PBS (p < 0.001) groups, as well as in the DMN-XII vs pff alone (p < 0.01) and PBS (p < 0.01) groups, at 6 months post lesion. *p < 0.05, **p < 0.01, ***p < 0.001 compared with pff only group at the same time point. Data were analyzed using multivariate ANOVA, n = 6-12 per group at 2 and 4 months post injection, n = 8-12 at 6 months post injection. Scale bars, 500 µm for a and e, and 100 µm for c of Math2 + glutamatergic neurons in the hippocampus are more susceptible than Prox1-containing glutamatergic neurons in the dentate gyrus, which is again related to the differential expression of endogenous α-synuclein [37]. However, this cannot be the whole story, as many brain areas have similar levels of endogenous α-synuclein yet exhibit selective PD pathology [38].
There are some limitations despite the obvious utility of the RVG9R:pff model. The modest level of nigral dopaminergic cell loss (15-20%) at 6 months post injection was inadequate to generate motor impairments. This likely also explains our MRI findings of no clear brain atrophy following adjustment for multiple comparisons. There were, however, some regional volumetric changes, which correlated with olfactory dysfunction, albeit at an exploratory threshold. For example, the negative correlation between cortical volume and olfactory performance in the RVG9R: pff-exposed rats is in good agreement with similar relationships in non-demented older adults with olfactory impairment [39,40] and in PD patients with no co-morbid dementia [41,42], although such changes are not found in all studies [43]. The positive correlations with changes in the pontine nuclei and BNST are slightly harder to interpret, since these regions do not have obvious involvement in olfactory processing. Nevertheless, such changes could be related to neuropathological changes, such as neuronal swelling and gliosis in these regions, associated with systemic RVG9R:pff administration. However, it should be noted that these correlations were exploratory and not corrected for multiple comparisons, and so should be viewed cautiously. Further studies will be required to confirm these relationships and elucidate the cellular correlates driving them. Overall though, our MR imaging endophenotype is consistent with the general lack of structural brain pathology detected by MRI in de novo PD patients [44]. Our data are also consistent with other reports that have injected pff directly into the brain, where the neuropathology and relationship to behavioural dysfunction may become clearer with longer follow-ups coupled to longitudinal in vivo MR imaging.
We did however detect a significant cholinergic neurodegeneration in the DMN-XII at 6 months in the RVG9R: pff rats. Although the presence of ubiquitin-positive degenerating neurites in the DMN-XII has been previously reported in PD patients [45], cholinergic neurodegeneration has not been established in this area. This could be attributed to the intrinsic biological difference between rodent and human cholinergic neurons in the DMN-XII, contributing to their differential susceptibility to α-synuclein pathology. However, a trend of cholinergic neurodegeneration (p = 0.073) in the DMN-XII has recently been reported in patients dying with incidental LB disease, but intriguingly not in idiopathic PD patients [22]. Prominent α-synuclein pathology could also be observed in the cortical regions at 6 month after RVG9R:pff injection. It is possible that transvascular delivery of α-synuclein may lead onto cognitive deficits, which define prodromal PD dementia [46], although we did not thoroughly explore this in our study. Furthermore, we used human pff as opposed to rodent pffs, which may be better to generate pathology, as homologous seeding is more efficient and pathogenic than cross-seeding [7].
Another limitation of the current study is that our analysis in the periphery was restricted to describing the behavioural and cellular deficits in the GI tract in the RVG9R:pff animals, and there may be additional pathological features that we did not examine (e.g. heart) [47,48]. Similarly, our histological examination of the CNS was restricted to the brain regions that have been well established to have pathology in early PD. This was also driven by our MRI analysis. Nonetheless, the datasets presented herein converge to suggest that we have a phenotype at 6 months of age consistent with early-stage PD.
In addition to developing this highly relevant new model of premotor PD, we have also shown that the RVG9R carrier is a versatile system for in vivo protein transduction to complement its previously defined uses in the transvascular delivery of small interfering RNA [13], non-coding RNA [14], and liposome-encapsulated nucleotides [49]. As we have Fig. 5 Site-specific and time-dependent α-synuclein pathology after a single systemic injection. a Representative figures showing the degree of α-synuclein phosphorylation (Ser129, red) and oligomerization (pFTAA, green) with confocal microscopy in various brain regions, as well as the presence of PK-resistant inclusions in bright-field microscopy in DMN-X. Lewy neurite-like structures, as well as LB-like structures with a dense core of phosphorylated α-synuclein surrounded by pFTAA + oligomers, was detectable in the SNc and DMN-X. n = 6-12 per group per time point. b Representative blots showing the development of total (syn42) and phosphorylated (Ser129) α-synuclein aggregation after sequential protein extraction, comparing between various regions of the CNS and ENS. At 6 months post injection, total and phosphorylated α-synuclein was detectable in the urea-soluble fraction in all brain regions except for the cerebellum (CB) and OB, as well as in the duodenal LMMP but not in the mucosa. c The aggregation and phosphorylation of α-synuclein in the CTx was also time dependent in the urea-soluble fraction. d Semi-quantitative analysis on the level of monomeric and HMW α-synuclein in the urea-soluble fraction, at 6 months post injection. There was a significant interaction effect between treatment and α-synuclein conformation (F 8,58 = 5.748, p < 0.001). Post hoc analysis revealed no significant difference in the expression of Syn42 α-synuclein monomers between 2 and 6 months post injection, although there was a significant increase in HMW total α-synuclein aggregation between 2 and 4 months post injection. On the contrary, the expression of Ser129 phosphorylated α-synuclein was significantly elevated between 2 and 4 months post injection, both in terms of monomeric (p = 0.031), and HMW α-synuclein (p < 0.001). Abbreviations include cortex (CTx), basal ganglia (BG), midbrain (Mid), brain stem (BS), and duodenal mucosa (Mus). n = 4 per group per time point. Nuclei were counterstained with Hoechst (blue). Scale bars, 30 µm for confocal microscopy and 60 µm for bright-field microscopy shown here, cargo complexing was simple and achievable within 30 min, in contrast to alternative approaches that require laborious procedures and extensive preparations to create chimeric, CPP-cargo fusion constructs. We have also ascertained that the same delivery system can be adapted to deliver a range of protein cargos including GFP, heterodimeric ChABC, as well as α-synuclein fibrils. In principle the same strategy could also be used to model other neurodegenerative proteinopathies [50].
In summary we have developed a new animal model in which the transvascular delivery of RVG9R:pff mirrors the earliest pathology and behavioural deficits of premotor PD. This opens up new ways to study PD pathogenesis and its treatment, and critically disease-modifying therapies that are likely to have their greatest impact at the onset of the disease process.
Author contributions WLK conceived the study. WLK, KS, JCFK, and RAB designed the experiments. WLK, KS, XH, SY, KH, and ACV performed the experiments. WLK, KS, TCW, ACV, YZ, OT, FIA, and RAB analyzed the results. All authors were involved in manuscript preparation.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons. org/licenses/by/4.0/.


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
Exome‐wide association of deltamethrin resistance in Aedes aegypti from Mexico

Abstract Aedes aegypti is the major vector of a number of arboviruses that cause disease in humans. Without vaccines or pharmaceuticals, pyrethroid insecticides remain the major tool for public health protection. Pyrethroid resistance is now widespread. Replacement substitutions in the voltage‐gated sodium channel (vgsc) that reduce the stability of pyrethroid binding account for most of the resistance, but metabolic mechanisms also inactivate pyrethroids. High‐throughput sequencing and the A. aegypti L5 annotated physical map has allowed interrogation of the exome for genes and single‐nucleotide polymorphisms associated with pyrethroid resistance. We exposed females of A. aegypti from Mexico to a deltamethrin discriminating dose to designate them as resistant (active after 1 h) or susceptible (knocked down with no recovery after 4 h). The vgsc on chromosome 3 had the highest association, followed by genes proximal to vgsc. We identified potential detoxification genes located singly (eg HPX8C) or within clusters in chromosome 2 [three esterase clusters, two of cytochrome P450 monooxygenases (CYP)] and chromosome 3 (one cluster of 16 CYP325 and seven CYP9 genes). Deltamethrin resistance in A. aegypti is associated with mutations in the vgsc gene and a large assortment of genes.

Introduction
Insecticide resistance in Aedes aegypti mosquito populations threatens our ability to control arboviral diseases such as dengue, chikungunya, and Zika. During the last decade, intense global reliance on pyrethroid insecticides in vector control programmes worldwide has resulted in widespread pyrethroid resistance in A. aegypti populations (da-Cunha et al. , 2005;Harris et al. , 2010;Marcombe et al. , 2012;Aponte et al. , 2013;Flores et al. , 2013). Assessment and characterization of pyrethroid resistance mechanisms are important to the establishment of insecticide resistance management plans as insecticidal spraying efforts increase.
Various mechanisms of resistance have been reported in insects, including the insensitivity of the target site, increased metabolism, reduced penetration throughout the cuticle and changes in behaviour that usually result in avoidance of insecticide exposure. Although insecticide resistance is multifactorial, most research has focused on target site insensitivity and metabolic mechanisms separately, whereas reduced penetration and behavioural changes are just starting to be understood (Zalucki and Furlong, 2017;Balabanidou et al. , 2018).
Pyrethroids are neurotoxins that target the voltagegated sodium channel (vgsc ) in insect neurons, resulting in a disrupted action potential that culminates with rapid knockdown and mortality. Specific amino acid substitutions at vgsc result in target insensitivity in resistant insects. Around 57 nonsynonymous mutations in the vgsc have been associated with channel insensitivity (Rinkevich et al. , 2013); interestingly, convergence of knockdown resistance (kdr) mutations across different insect genera validate the important role of this mechanism. In A. aegypti , 11 nonsynonymous substitutions from different areas of the world have been reported (Rinkevich © Saavedra-Rodriguez et al. , 2018). Currently, five have been confirmed to contribute to insecticide resistance by electrophysiology Hirata et al. , 2014;Haddi et al. , 2017). Interestingly, co-existence of kdr mutations at different loci is common in insect species; furthermore, in A. aegypti these kdr mutations seem to be geographically distributed, conferring different levels of resistance (Hirata et al. , 2014;Saavedra-Rodriguez et al. , 2018). The three loci examined in this paper are 410, 1016 and 1534 corresponding to the position of amino acid substitutions. A valine at locus 410 (V410) confers susceptibility, whereas leucine (L410) confers resistance. A valine at locus 1016 (V1016) confers susceptibility, whereas isoleucine (I1016) confers resistance. A phenylalanine at locus 1534 (F1534) confers susceptibility, whereas cysteine (C1534) confers resistance.
Increased metabolism can result from overexpression (by amplification or cis-trans regulation) of insecticide metabolizing enzymes or by enzyme variants that provide higher or insecticide-specific catalytic properties (Hemingway and Karunaratne, 1998). Four enzyme families are associated with insecticide metabolism: cytochrome P450 monooxygenases (CYP), carboxyl/choline esterases (CCE), glutathione-S -transferases (GSTs) and reduction/oxidation (Red/ox) enzymes. Identification of detoxification genes in A. aegypti confirm the high redundancy of genes potentially involved in insecticide metabolism, describing 160 CYP, 44 CCE and 26 GSTs (Strode et al. , 2008). With improved microarray technology, whole transcriptomes had been interrogated to identify candidate genes conferring insecticide resistance in field or in artificially selected A. aegypti strains Marcombe et al. , 2009a;Bingham et al. , 2011;Bariami et al. , 2012;Grisales et al. , 2013;Kasai et al. , 2014). Unsurprisingly, genes and their expression patterns varied among regions of the globe. Nonetheless, a few common genes were successfully confirmed to produce insecticide-metabolizing proteins (Stevenson et al. , 2012).
As next-generation sequencing becomes more affordable, whole genomes or transcriptomes are interrogated for association with specific phenotypic traits. Few studies using this technology in mosquitoes have assessed the transcriptional changes and replacements identified in permethrin-resistant A. aegypti (David et al. , 2014;Faucon et al. , 2015). A recent complete partially annotated physical map of the A. aegypti genome (Matthews et al. , 2018) will allow us for the first time to perform genome-wide association studies. In the present study, we use an exome-enrichment genomic DNA sequencing approach (Juneja et al. , 2015) to identify singlenucleotide polymorphisms (SNPs) associated with deltamethrin resistance. We analysed an A. aegypti collection from Viva Caucel from Yucatán State in southern Mexico that was subjected to permethrin selection in the field from 1998 to 2010. Our hypothesis was that replacements in detoxification genes that have significantly different frequencies between resistant and susceptible mosquitoes are indicative of their potential involvement in deltamethrin resistance or are due to their physical proximity to vgsc on chromosome 3. The expected vgsc replacement polymorphisms were identified in the dataset and also allowed us to identify genome regions surrounding vgsc that are likely to have been involved in a genetic sweep.

Resistance phenotype
Deltamethrin resistance was phenotyped by exposing 390 A. aegypti adults from a field resistant strain collected in Viva Caucel to bottles coated with 3.0 µg of deltamethrin during 1 h. This concentration discriminated three phenotypes: (1) 149 mosquitoes were 'knockdown resistant' (38%), (2) 148 were 'susceptible' (38%) and (3) 93 'recovered' (24%) (Fig. 1). In this study, we compared pooled mosquitoes from the 'knockdown resistant' and 'susceptible' groups ( Fig.  1). We excluded the 'recovered' group to remove the uncertainty of assigning the incorrect phenotype. Also, we assumed that more contrasting phenotypes would allow us to get more contrasting SNPs underlying deltamethrin resistance.
A Benjamini-Hochberg correction (Benjamini and Hochberg, 1995) for false discovery rate (FDR) was applied to chromosomes 1, 2 and 3 separately (α = 0.01). A total of 14 532 SNPs (0.206% of the targeted exome) surpassed the FDR thresholds of 3.39, 3.52 and 3.22 for chromosomes 1, 2 and 3 respectively. These SNPs occurred in the 3′-untranslated region (UTR; n = 321), 5′-UTR (n = 443), coding sequence (n = 10 359; of which 2614 were nonsynonymous and 7745 were synonymous substitutions), intergenic regions (n = 40), noncoding RNA (n = 118) and introns (n = 3251). Fig. 2 shows the distribution of SNPs with −log 10 (prob) values above the FDR along chromosomes 1, 2 and 3. We provide a complete list of significant SNPs in Tables S2, S3 and S4. The approach used in this paper assumes that significant SNPs are independently distributed (in linkage equilibrium). This assumption is invalid, but we have no understanding of linkage disequilibrium blocks or haplotype structure in this or other mosquito species. Instead, we focus our results and discussion on those genes with the largest −log10(prob) values.
Our results showed that V410L located at IIS6 in vgsc had a high association [−log 10 (prob) = 17] with deltamethrin resistance. This mutation was recently associated with pyrethroid resistance in A. aegypti (Haddi et al. , 2017;Saavedra-Rodriguez et al. , 2018). To determine whether the library sequences accurately estimated the population of V410L allele frequencies, we performed a melting curve PCR SNP assay to genotype the individuals used to construct the four libraries (Saavedra-Rodriguez et al. , 2018). Fig. 3 shows the comparison of V410L allele frequencies obtained by library sequencing and by melting curve PCR SNP genotyping. Allele frequencies were not significantly different between replicates in any of the comparisons. Individual genotyping shows that the frequency of the resistant allele L410 was 0.88 and 0.94 in the knockdown-resistant individuals used in the replicates and 0.28 for the susceptible individuals used in the replicates (Fig. 3). In the sequencing libraries, L410 had a −log 10 (prob) = 17; frequencies in the knockdown-resistant libraries for replicates 1 and 2 were 0.88 and 0.73 respectively. In the susceptible individuals, frequencies were 0.10 and 0.3 in replicates 1 and 2 respectively. Heterogeneity χ 2 tests showed that the allele frequencies were similar regardless of the detection technique used, suggesting Testing for a genetic sweep by vgsc Because 45% of the 168 SNPs with high −log 10 (prob) > 10 were located in close proximity to vgsc (Fig. 2C), we analysed a 1.4 Mbp region upstream/downstream of this gene. By observing SNPs with −log 10 (prob) > 10 in this region, the sweep extended from a TATA-box binding protein upstream to vgsc to a vanin-2 protein downstream to vgsc . Genes in this sweep included a ubiquitin protease (OTU domain-containing protein 7B), a cytoskeleton binding protein (moesin/ezrin/radixin), a protein that forms a structural framework for protein-protein interactions (leucine-rich repeat-containing protein 47), a ribosomal biogenesis protein (large subunit GTPase) and two vanins involved in nitrogen compound metabolic processes (Fig. 4). Since none of these protein-coding genes have been functionally or genetically associated with insecticide resistance, we tested the possibility of a genetic sweep by rapid positive selection at vgsc rather than being selected due to their own activity as resistance loci. Fig. 4 shows the expected heterozygosity H exp of SNPs at vgsc and surrounding genes in knockdown-resistant and susceptible mosquitoes, where H exp = 1 − ∑ n i=1 p i 2 and n is the number of alternate nucleotides at a site. Heterozygosity at position L410 was low in knockdown-resistant mosquitoes and high in knockdown-susceptible mosquitoes. This is expected to occur when beneficial mutations (at vgsc ) are rare; but once a beneficial mutation has occurred it increases in frequency rapidly, thereby drastically reducing genetic variation in the population. To determine whether this pattern was genome wide, we tested for this pattern throughout all three chromosomes. Table 1 and Fig. S1 display the result of three pairwise t -tests, one for each chromosome testing whether the mean of the difference between heterozygosities in resistant and susceptible mosquitoes is different from zero. For each pairwise comparison of SNPs with −log 10 (prob) greater than the FDR we subtracted the expected heterozygosity of susceptible mosquitoes, Het (dead), from the heterozygosity of resistant mosquitoes, Het (alive). When Het (alive) > Het (dead) then the difference is positive, but when Het (dead) > Het (alive) the difference is negative. We would expect genes in surviving mosquitoes to have been swept and, therefore, have a lower heterozygosity. Conversely, we would expect genes in dead mosquitoes (homozygous susceptible) not to have been swept or swept at a lower rate (heterozygotes) and, therefore, have a higher heterozygosity. Along chromosomes 1 and 3 the mean difference was negative and significant in the paired t -test (Table 1), consistent with genetic sweeps.
However, along chromosome 2 the mean difference was positive and significant in the paired t -test (Table 1). Most importantly, in all three chromosomes there are distinct peaks where Het (alive) > Het (dead). This suggests that genes along all three chromosomes have not been swept and that the pattern in Fig. 4 is not genome wide.

Insecticide metabolism: detoxification genes
We identified 268 insecticide detoxification genes  in the AaeL5 assembly map and questioned which significant SNPs in our libraries were included in this categorical group. Genes with insecticide detoxification function surpassed the FDR cut-off calculated for each chromosome. Fig. 5 shows the −log 10 (prob) calculated for the SNPs comprised in 2, 57 and 46 detoxification genes dispersed across chromosomes 1, 2 and 3 respectively. Overall, 62% of the SNPs are classify as CYP genes, 20% as CCE, 16% as Red/ox and 4% as GSTs. In chromosome 1, only 13 SNPs were significant, with 11 SNPs in HPX8C and two SNPs in an oxidase peroxidase (HPX6) (Fig. 5A). Most SNPs were annotated as synonymous substitutions, except for an S766A in HPX8C [−log 10 (prob) = 3.76]. Both genes are currently annotated as chorion peroxidases. Several peroxidase isozymes have been described in Drosophila melanogaster as cyclooxygenase-like protein for both the actin filament bundle formation and the cortical actin strengthening during Drosophila oogenesis.
Because previous studies focused mostly on gene expression, we compared our list of genes with significant SNPs with 22 studies that had measured gene expression in mosquitoes exposed to sublethal doses of insecticides (Poupardin et al. , 2008;Riaz et al. , 2009;David et al. , 2010). We also included published articles that compared insecticide-resistant versus susceptible strains Marcombe et al. , 2009b;Saavedra-Rodriguez et al. , 2012;David et al. , 2014;Reid et al. , 2014;Estep et al. , 2017) and, most recently, allelic variants (David et al. , 2014;Dusfour et al. , 2015;Faucon et al. , 2015;. Fig. S3 shows the list of significant detoxification genes in our study and the expression that has been reported in previous studies. For example, in our study, CYP6N12 (LOC5571529) was significantly associated with deltamethrin resistance; this gene was significantly down/overexpressed in 12 studies, becoming a clear candidate for further study.
Furthermore, Campbell et al. (2019) analysed a compilation of high-throughput sequencing data from three pyrethroid treatment groups and interrogated them to identify functional categories that were overrepresented among all coding genes. Of eight functional categories, three were overrepresented, as calculated by a hypergeometric test, which tests the probability of a gene being represented more often than expected by chance. The three overrepresented groups were metabolism, transport and signalling. The detox gene subset would encompass both the metabolic category and the redox category, which was not overrepresented.

Discussion
With increased insecticide spraying to protect public health, the detection and characterization of resistance loci are a crucial requisite for adequate monitoring of the rise of insecticide resistance that could hamper vector control efforts. Currently, kdr is a major mechanism of pyrethroid resistance; around 11 amino acid replacements in vgsc have been associated with this trait in A. aegypti , and quantification of their frequencies in field populations is becoming a common practice to monitor pyrethroid resistance. Unfortunately, continued use of pyrethroids in the field is leading to the fixation of kdr mutations and possibly the selection of novel or compensatory mutations in vgsc and other genes associated with pyrethroid resistance. Identifying these novel or compensatory mutations is necessary to comprehensively monitor pyrethroid resistance in field populations. In this study, 52-80% of vgsc was sequenced across libraries (out of 7794 nucleotides covered by enrichment probes covering exons 5-30 and introns 4-29) and 19 SNPs were significantly associated with deltamethrin resistance. Notably, three nonsynonymous mutations in vgsc had significant association. Two are novel mutations, and their role in pyrethroid resistance is unknown. Contrastingly, V410L was highly associated with resistance in our study, bringing further evidence in its role in pyrethroid resistance (Haddi et al. , 2017). We recently tracked V410L in a population from Mexico collected in 2002, and frequencies have increased in the last decade in strong linkage disequilibrium to mutations I1016 and C1534 (Saavedra-Rodriguez et al. , 2018). Interestingly, I1016 and C1534 were not significantly associated with deltamethrin resistance in our study. To understand this observation, we individually genotyped mosquitoes included in our libraries and calculated allele frequencies of I1016 and C1534. I1016 had identical frequencies to those found in L410; however, coverage (>25 reads) was not reached at this site. The improvement of library preparation and DNA pool optimization will address the identification of this locus in association with pyrethroid resistance. Contrastingly, the frequency of C1534 was high (>0.9) in both susceptible and resistant libraries, so that −log 10 (prob) values were not high enough to surpass the chromosome FDR threshold. In Mexico, strong selection has pushed C1534 to fixation in many field populations (Vera-Maloof et al ., 2015). Therefore, the use of C1534 as a marker for pyrethroid resistance is no longer informative in Mexico, even though the role of C1534 in pyrethroid resistance has been validated by electrophysiology studies Haddi et al. , 2017). Furthermore, C1534 appears to be necessary for other mutations (I1016) to occur and increase in frequency in the population (Vera-Maloof et al ., 2015). Interestingly, we identified a possible genetic sweep driven by selection of vgsc nonsynonymous mutations associated with pyrethroid resistance. Similar genetic sweeps have been reported in D. melanogaster driven by selection of acetylcholinesterase mutations associated with carbamate and organophosphate resistance and by DDT resistance mediated by an insertion of a transposon in the 5′ regulatory region of CYP6g1 (Garud et al. , 2015). Similarly, strong signals of recent positive selection at several genes that are known to have a role in resistance have been identified in Anopheles , including vgsc , a cluster of GST genes including Gste2 (implicated in the metabolism of DDT and pyrethroids), and Cyp6p , a cluster of genes that include Cyp6p3 , which is upregulated in permethrin and bendiocarb-resistant mosquitoes (Anopheles gambiae 2017 Genomes Consortium, 2017). In A. aegypti , evidence of a genetic sweep (or hitchhiking effect) was demonstrated in a region of chromosome 1 that contains organophosphate-resistance-conferring genes (Yan et al. , 1998). The genetic sweep was reflected by low DNA polymorphisms and gene deletions for loci surrounding the EST-4 locus gene region. In contrast, a hitchhiking effect associated with DDT was not evident for the para (vgsc ) locus, in which heterozygosity was similar to other regions in the genome. This result was consistent with the hypothesis that, in the years since DDT was abandoned, the populations have had time to re-equilibrate. Luckily, with the A. aegypti physical map completed, we can test the gain/loss of genomic variability in populations exposed to different pressures of insecticides.
The role of metabolism in conferring pyrethroid resistance is unquestionable. Even if kdr mutations at vgsc protect mosquitoes from knockdown, pyrethroids are metabolized before reaching vgsc or after dissociation from vgsc by different enzymes to less toxic and more hydrophilic metabolites for excretion. The metabolism of pyrethroids occurs by ester hydrolysis (by CCE) and oxidation at methyl, methylene, alkenyl or aryl substituents (CYP), and ≥80 metabolites have been identified from cis and trans -permethrin (Casida et al. , 1983). To date, fewer than 10 genes have been functionally characterized and confirmed to metabolize pyrethroids in A. aegypti , including members of the CYP9 family (Stevenson et al. , 2012), CYP9M6 and CYP6BB2 , whereas other enzymes, such as CYP6Z8 and the Red/ ox (aldehyde dehydrogenases), catalyse the oxidation of permethrin intermediates (Chandor-Proust et al. , 2013;Lumjuan et al. , 2014). Noticeably, genes CYP6Z8, CYP9J26-28 and CYP9M6 belong to clusters of tandem genes. In our study, a large cluster of CYP325 genes at the p arm of chromosome 3 was strongly associated with deltamethrin. Although, the physiological function of the CYP325 in insects still remains unclear, some CYP325s (eg CYP325-E, K and G) are conserved among mosquito species, some only expanded in the Culex and Aedes genomes (CYP325X and Y) and some are species specific (Yan et al. , 2012). In our review, at least 13 CYP325 genes (Fig. S3) have been associated with changes of expression (two to threefold) in insecticide-resistant A. aegypti. In one study, CYP325N1 was 313-fold more expressed in a pyrethroid-resistant strain from Puerto Rico (Estep et al. , 2017). In other insects, only CYP325A3 was reported in a permethrin-resistant strain of Anopheles gambiae (David et al. , 2005); however, no clear orthologues were found in Aedes or Drosophila . Further research of the CYP325 cluster might help elucidate models for mechanisms of cluster expression and adaptation to specific insecticides, if any.
Similarly, identification and isolation of specific esterases that metabolize pyrethroids is ongoing, mainly because esterases have unspecific substrates that complicate protein isolation. In this study, several clustered CCE genes in chromosomes 2 (CCEaeA, CCEaeB and CCEaeD) and 3 (CCEunk and CCEglt) showed high association with deltamethrin resistance. CCEaeA family members are commonly reported in expression studies (Fig. S3). Genes from this cluster showed strong amplification in two deltamethrin-resistant strains from Thailand but were not present in two pyrethroid-  (Faucon et al. , 2015). Similarly, CCEae3A and 6A have been associated with temephos resistance in an A. aegypti strain for Thailand that was also resistant to pyrethroids . Although amplification of these CCE clusters are a possible mechanism to overexpress enzymes that metabolize insecticides, our study and others (Faucon et al. , 2015; have also identified allele variants within these genes. The roles of these nonsynonymous mutations in catalytic activity is of great interest. For example, amino acid substitutions increased diazinon hydrolysis by esterase E7 in Lucilia cuprina (Newcomb et al. , 1997) and M. domestica (Claudianos et al. , 1999), whereas a different substitution in the same gene was responsible for malathion resistance in L. cuprina (Campbell et al. , 1998). Characterization of CCE genes and their allelic variants will increase our knowledge of possible routes of metabolism.
We have mapped genomic regions associated with deltamethrin resistance in mosquitoes from Mexico. Our results validate the major role of vgsc in pyrethroid resistance, but fixation of replacement mutations in vgsc suggest that kdr is becoming well established in field populations. Additional mechanisms of deltamethrin resistance might play increasingly important roles in the future, following genomic changes in field populations that are routinely pressured with insecticides. It is therefore necessary to identify genes and possible mechanisms being selected by pyrethroids in general. So far, we have identified several detoxification genes that are highly variable in resistant mosquitoes. Research to further associate allelic variants with deltamethrin resistance is ongoing.

Mosquito samples, bioassays and DNA isolation
The A. aegypti population Viva Caucel from Yucatán State in southern Mexico (longitude −89.71827, latitude 20.99827) was collected in 2011 by Universidad Autónoma de Yucatán. This mosquito population was under field selection by a type-1 pyrethroid (permethrin) in vector control programmes from 1999 to 2010. We calculated the deltamethrin concentration to kill 50% of the mosquitoes (LC 50 ) by bottle bioassays to establish the level of resistance in the Viva Caucel population (F 3 ) relative to the susceptible New Orleans (NO). Around 20 non-blood-fed females (2-3 days old) of the Viva Caucel F 3 and NO were exposed to seven deltamethrin (Chem Service, West Chester, PA, USA) concentrations for 1 h, in three replicates. Mortality was scored following a recovery time of 24 h. The deltamethrin LC 50 at 24 h was 47.6-fold higher in the Viva Caucel population than in the NO susceptible strain (10.49 µg vs. 0.22 µg respectively) (Fig. S2). The NO strain was not used further in this study.
To discriminate the deltamethrin-resistant and susceptible individuals within the Viva Caucel population, 390 adult mosquitoes were exposed to 3 µg of deltamethrin for 1 h. This concentration is below the LC 50 calculated at 24 h post-exposure; however, 3 µg allowed us to discriminate the three phenotypes at 4 h of observation (instead of 24 h) while ensuring good DNA quality in 'susceptible' mosquitoes. Additionally, 4 h is a reasonable time for recovery, since longer periods increase the chances of predation and dehydration in knocked-down mosquitoes in the field.
The bioassay was performed by aspirating groups of 50 non-blood-fed female mosquitoes (3-4 days old) into deltamethrin-coated bottles for 1 h. After this time, active mosquitoes were transferred to 1 pint cardboard cups and placed into an incubator (28 °C and 70% humidity) for 4 h. These mosquitoes constituted the 'knockdown-resistant' group. Knocked-down mosquitoes were then transferred to a second cardboard cup. Four hours later, newly recovered mosquitoes were aspirated, frozen and labelled as 'recovered'. These were not used in the current study. The knocked-down mosquitoes that remained inactive at 4 h post-treatment were scored as 'susceptible'. DNA was isolated from individual mosquitoes that were categorized as either resistant or susceptible (Black and DuTeau, 1997) and resuspended in 150 µl of TE buffer (10 mM Tris hydrochloride, 1 mM ethylenediaminetetraacetic acid pH 8.0).

Sample pooling and library preparation
We constructed four genomic DNA libraries. A library was generated from the DNA of 25 knockdown-resistant females (Fig. 1), and this was replicated in a second set of 25 mosquitoes. The same was done with susceptible females (susceptible replicates 1 and 2) (Fig. 1). Before pooling, DNA from each individual mosquito was quantified using the Quant-IT Pico Green kit (Invitrogen, Eugene, OR, USA), and ~40 ng from each individual DNA sample was used for a final DNA pool of 1 µg. Pooled DNA was sheared and fragmented by sonication to obtain fragments between 300 and 500 bp (Covaris Ltd, Brighton, UK). We prepared one library for each of the four DNA pools following the Low Sample (LS) protocol from the Illumina TrueSeqDNA PCR-Free Sample preparation guide (Illumina, San Diego, CA, USA).

Whole exome hybrid capture
Because of the large 1.38 Gbp size of the A. aegypti genome (Nene et al. , 2007), we used an exome capture method. To identify the proportion of L5 genes covered by the capture probes, the analogous L5 probe sequences were blasted (blastn) against the L5 reference transcriptome. The returned genes were cross-checked against the total number of L5 genes (18 081). Thus, 77% (13 941/18 081) of L5 genes were accounted for among the probe set. We performed an exome-capture hybridization to enrich for coding sequences using custom SeqCap EZ Developer probes (Roche NimbleGen Inc., Madison, WI, USA). Probes covered protein coding sequences (not including UTRs) in the AaegL1.3 genebuild using the exonic coordinates specified previously (Juneja et al. , 2015). In total, 26.7 Mb of the genome (2%) containing 15 735 genes was targeted for enrichment. TruSeq libraries were hybridized to the probes using the xGen®Lock®Down recommendations (Integrated DNA Technologies

Analysis pipeline
The raw sequence files (*.fastq) for each pair-ended gDNA library were aligned to a custom reference physical map generated from the assembly AaegL5 (Matthews et al. , 2018) using the package gsnap (Genomic Short-read Nucleotide Alignment; Wu and Nacu, 2010), which performed an SNP-tolerant alignment of the resulting (39-42) × 10 6 paired reads/library. Additional processing details for our in-house fortran pipeline are provided in Dickson et al. (2017). Briefly, we removed SNPs that were significantly different between replicates. Next, we grouped replicates that share the same sites together to build contingency tables and the heterogeneity χ 2 was calculated with (n − 1) nucleotides degrees of freedom. The probability derived from that analysis was −log 10 transformed to provide a −log 10 (prob) value. A Benjamini-Hochberg correction for FDR (Benjamini and Hochberg, 1995) was applied for chromosomes 1, 2 and 3 separately (α = 0.01) using sas.

Validation of nonsynonymous substitution V410L in vgsc
We determined the individual genotype of the mosquitoes used in our four libraries (n = 100) at two nonsynonymous substitutions in the vgsc . We performed melting curve analyses of allele-specific PCR products to genotype at locus V1016I (Saavedra-Rodriguez et al. , 2007) and a recently reported mutation at locus V410L (Haddi et al. , 2017) following an allele-specific PCR described in Saavedra-Rodriguez et al. (2018).

Supporting Information
Additional supporting information may be found in the online version of this article at the publisher's web site: Figure S1. Distribution of the differences in heterozygosity in knockdown resistant versus susceptible mosquitoes (H exp alive (HETA) H exp dead (HETD)) at each significant SNP in chromosomes 1, 2 and 3. Figure S2. Deltamethrin dose response curve for female Ae. aegypti from Viva Caucel and the reference strain from New Orleans. Exposure was performed using bottles coated with increasing concentrations of deltamethrin. Exposure time was 1 h and mortality was recorded at 24 h. Figure S3. List of significant detoxification genes in our study and the expression that has been reported in previous studies. Upregulated genes are highlighted in orange, downregulated are highlighted in blue, downregulated or upregulated in the same study are highlighted in grey. A number represents the number of polymorphisms found at a particular SNP in that gene. Table S1. Number of reads aligned AaeL5 across chromosome 1, 2 and 3 for each replicated libraries. The number of common reads between Table S2. List of Individual SNP, physical position of each SNP, gene description, mutation, frequencies of alternant nucleotide in resistant and susceptible groups, log 10 (prob), nucleotide, type of SNP, substitution, amino acid residue and codon position of each SNP for chromosome 1. Table S3. List of Individual SNP, physical position of each SNP, gene description, mutation, frequencies of alternant nucleotide in resistant and susceptible groups, log 10 (prob), nucleotide, type of SNP, substitution, amino acid residue and codon position of each SNP for chromosome 2. Table S4. List of Individual SNP, physical position of each SNP, gene description, mutation, frequencies of alternant nucleotide in resistant and susceptible groups, log 10 (prob), nucleotide, type of SNP, substitution, amino acid residue and codon position of each SNP for chromosome 3. Table S5. List of non synonymous mutations in deltamethrin resistant Ae. aegypti . The genes, description, SNPID, residue number, frequencies in resistant and susceptible and log 10 (prob) is shown for chromosomes 1, 2 and 3.