Knockdown resistance (kdr) genotypes and collection information for Aedes aegytpi from Iquitos, Peru (2000 - 2017)

This study describes the evolution of knockdown resistance (kdr) haplotypes in Aedes aegypti in response to pyrethroid insecticide use over the course of 18 years in Iquitos, Peru. Based on the duration and intensiveness of sampling (~10,000 samples), this is the most thorough study of kdr population genetics in Ae. aegypti to date within a city. We provide evidence for the direct connection between programmatic citywide pyrethroid spraying and the increase in frequency of specific kdr haplotypes by identifying two evolutionary events in the population. The relatively high selection coefficients, even under infrequent insecticide pressure, emphasize how quickly Ae. aegypti populations can evolve. In our examination of the literature on mosquitoes and other insect pests, we could find no cases where a pest evolved so quickly to so few exposures to low or non-residual insecticide applications. The observed rapid increase in frequency of resistance alleles might have been aided by the incomplete dominance of resistance-conferring alleles over corresponding susceptibility alleles. In addition to dramatic temporal shifts, spatial suppression experiments reveal that genetic heterogeneity existed not only at the citywide scale, but also on a very fine scale within the city.

Methods
Entomological Surveys

Control History. During the past two decades, researchers have preserved mosquito specimens collected throughout the city of Iquitos in a repository (e.g. Cavany et al., 2020; Cromwell et al., 2017; Getis et al., 2003; Gunning et al., 2018; LaCon et al., 2014; Lenhart et al., 2020; Morrison et al., 2004a; Morrison et al., 2004b; Morrison et al., 2006; Morrison et al., 2008; Reiner et al., 2019; Schneider et al., 2004; Tun-Lin et al., 2009). This repository holds samples that were collected prior to pyrethroid application (2000 – 2002), during citywide pyrethroid use (2002 – 2014), and after pyrethroids were discontinued by the Ministry of Health (2014 – 2017). Prior to 2002, no citywide insecticide spraying targeted to control Ae. aegypti occurred. Once targeted control began, multiple sub-classes of pyrethroid insecticides were sprayed by the Iquitos Ministry of Health inside homes from 2002 – 2014. These sub-classes included deltamethrin, cypermethrin, alpha-cypermethrin, lambda-cyhalothrin, and alpha-cypermethrin + pyriproxyfen (Suppl. Tbl. 1). It was uncommon for residents to spray their own homes (A. Morrison, personal communication). In 2014, pyrethroids were discontinued in favor of malathion due to the development of phenotypic pyrethroid resistance (Gunning et al., 2018). During the years of pyrethroid applications, spraying occurred at an average of 3.25 treatments per year, with spraying occurring within a one-month period in most years. Because the applications were Ultra Low Volume space sprays with no residual, most generations of the mosquitoes in any year were not exposed to insecticide and mosquitoes were collected independent of whether sprays had been conducted (Ritchie et al. 2021).

Temporal Collections. Aedes aegypti were collected and stored at -80˚C by NAMRU-6 and University of California at Davis personnel since the late 1990s. Specimens dating back to the year 2000 were available for study. Mosquitoes were collected by backpack aspirator (Clark et al., 1994) prior to June 2009 and by Prokopack Aspirator following June 2009 (Vasquez-Prokopec et al., 2009; Reiner et al., 2019). Each mosquito in the repository was identified to species, sex, collection date, and collection site. Each collection site, typically an individual household, was associated with GPS coordinates (Fig. 1).

Spatial Collections. Intense suppression experiments based on pyrethroid spraying were conducted in 2013 and 2014 (Gunning et al., 2018) to test the predictions of a detailed Ae. aegypti population dynamics model (Magori et al., 2009). In brief, two areas of the city were identified as having relatively high densities of Ae. aegypti and were configured spatially in a way that allowed for a central spray sector with an outer buffer sector to act as an experimental control region (Fig. 2). To limit the impact of migration on resistance allele frequency, site dimensions were selected to be 3 – 5 times larger than the expected Ae. aegypti lifetime flight distance of approximately 150 m (Harrington et al., 2005).

The 2013 study site covered approximately 750 m x 450 m and contained 1,163 houses. Baseline samples were collected in January 2013. Systematic sampling began on 22 April 2013 and continued for 16 weeks until 8 August 2013. From 29 April 2013 to 3 June 2013, six weekly non-residual, indoor ultra-low volume (ULV) cypermethrin treatments were applied in the treatment sector.

The 2014 study site was larger and covered an approximate 600 m x 600 m area and contained 2,166 houses. Systematic sampling was conducted over a longer period of 44 weeks. ULV spraying of cypermethrin was performed from 28 April 2014 through 2 June 2014 in a similar manner as in the 2013 study. In addition to the study spray in 2014, a citywide spray was conducted in response to a dengue outbreak in February 2014, during which homes in both the experimental and buffer sectors were sprayed with pyrethroids. Throughout both suppression experiments, mosquitoes were collected and stored as described above.

DNA Extractions and Quantification

Male mosquitoes were chosen for genetic analysis throughout this study because female mosquitoes were typically the focus of virology and epidemiological studies and, therefore more males were available in the repository. Using females would also have brought the risk of genomic contamination from male mosquitoes (via insemination) and from humans (via human blood feeding). Males are expected to share similar allele frequencies with females because the VGSC is not sex-linked.

Whole mosquitoes were transferred from Iquitos, Peru to Raleigh, North Carolina, USA with permits from Peruvian and US authorities. Samples were stored at -80˚C prior to genomic DNA (gDNA) isolation and at -20˚C after gDNA isolation. Genomic DNA was extracted from whole male Ae. aegypti by one of two methods: Qiagen DNeasy blood and tissue kit (cat: 69582) or Canadian Center for DNA Barcoding protocol. In brief, for the Qiagen DNeasy kit protocol, whole male mosquitoes were homogenized and incubated in lysis buffer and proteinase K overnight at 55˚C. Following incubation and removal of chitinous material, RNase A treatment was performed to remove RNA contamination for both isolation methods. Then, the standard Qiagen protocol of washes was followed. Final samples were eluted two times in 150 µl warm dH2O (Invitrogen Cat #: 10977-015). A modified Canadian Center for DNA Barcoding (2020) protocol was also used for some mosquito DNA isolations to reduce costs while maintaining quality genomic DNA extractions. Samples were homogenized, incubated, and RNase A treated as described above before the lysate was passed through the filter of an AcroPrep™ PALL2 plate (Cat #: PALL 5053) to bind the gDNA. The filter was washed with Protein Wash Buffer to remove remaining proteins and then washed with cold Wash Buffer to remove additional contaminates. The filter was allowed to dry to ensure that no ethanol remained to interfere with DNA yield. Finally, two washes of 75 µl warm dH2O (Invitrogen Cat #: 10977-015) were performed to elute a final volume of 175 µl gDNA.

Quantification of gDNA was performed using a Quant-iT PicoGreen dsDNA assay (ThermoFisher Scientific - P11496) and samples were read on a Synergy H1 Hybrid Plate Reader (BioTek Instruments, Inc.) in the Genomic Sciences Laboratory at North Carolina State University (GSL).

Genotyping

Allele-specific quantitative PCR and melting curve analysis (AS-PCR) was used to genotype all mosquitoes in duplicate for each of the mutations most commonly found in Central and South America (V1016I and F1534C). If the two reactions were not scored identically, the sample was discarded from further analysis. Mismatches were rare and typically due to non-amplification of a sample or because certain criteria for scoring were not met; i.e., melting peak did not cross threshold. A smaller number (n=92) of individuals were additionally genotyped at the V410L locus to verify the strong linkage disequilibrium that has been previously reported between it and locus V1016I (Saavedra-Rodriguez et al., 2018). Each mosquito genotyped at the V410L locus was also genotyped twice to ensure accuracy.

Genotyping of V1016I. AS-PCR for the V1016I locus was based on the method reported by Saavedra-Rodriguez et al. (2007) and modifications to the I1016R primer made by the Entomology Branch at the Centers for Disease Control and Prevention (CDC), Atlanta, USA (A. Lenhart, personal communication). The PCR volume was reduced to 10 µl per reaction and contained 2.5 µl of dH2O, 0.5 µl of each primer at 10 µM (V1016F, I1016F, I1016R), 5 µl of PerfeCTa SYBR Green Supermix (Quanta – 95054-02K), and 1 µl of template. The primer sequences for V1016F and I1016F used are reported in Saavedra-Rodriquez et al. (2007), but the primer sequence for I1016R was modified to: 5’ - TGA TGA ACC SGA ATT GGA CAA AAG C – 3’ (CDC, personal communication). Samples were genotyped on a BioRad CFX384 Real-Time PCR machine in the GSL, with the following thermal conditions: step 1 - 95˚C for 3 minutes, step 2 – 95˚C for 10 seconds, step 3 – 60˚C for 10 seconds, step 4 – 72˚C for 10 seconds, step 5 - Go to step 2, 39 times, step 6 – 95˚C for 10 seconds, step 7 - Melting Curve 65˚C - 95˚C, increment 0.2˚C per 10 sec plus a plate read.

Genotyping of F1534C. AS-PCR for the F1534C locus was performed following the method reported by Yanola et al. (2011) with the following modifications. The PCR volume was reduced to a total of 10 µl per reaction and contained 5 µl dH2O, 0.2 µl C1534F primer, 0.4 µl F1534F primer, 0.4 µl F1534R primer, 3.0 µl of PerfeCTa SYBR Green Supermix (Quanta – 95054-02K), and 1 µl of template. Samples were genotyped on a BioRad CFX384 Real-Time PCR machine in the GSL, with the following thermal conditions: step 1 - 95˚C for 2 minutes, step 2 – 95˚C for 30 seconds, step 3 – 60˚C for 30 seconds, step 4 – 72˚C for 30 seconds, step 5 - Go to step 2, 34 times, step 6 – 72˚C for 2 minutes, step 7 - Melting Curve 65˚C - 95˚C, increment 0.2˚C per 10 sec plus a plate read.

Genotyping of V410L. The AS-PCR for the V410L locus was based on a protocol developed by K. Saavedra and shared via the Entomology Branch, CDC (A. Lenhart, personal communication). The total volume for each reaction was reduced to 10 µl: 3.8 µl dH2O, 0.05 µl Val410 primer (50 µM) 5’ – GCG GGC AGG GCG GCG GGG GCG GGG CCA TCT TCT TGG GTT CGT TCT ACC GTG – 3’, 0.05 µl Leu410 primer (50 µM) 5’ – GCG GGC ATC TTC TTG GGT TCG TTC TAC CAT T – 3’, 0.1 µl Rev410 primer (50 µM) 5’ – TTC TTC CTC GGC GGC CTC TT – 3’, 5.0 µl PerfeCTa SYBR Green SuperMix (Quanta – 95054-02K), and 1 µl template. Thermal conditions were performed on the BioRad CFX384 Real-Time PCR machine in the GSL, with the following thermal conditions: 95˚C for 3:00, 40 cycles of (95˚C for 0:10, 60˚C for 0:10, 72˚C for 0:30), 95˚C for 0:10, melting curve 65˚C – 95˚C increasing in increments of 0.2˚C per 10 sec plus a plate read.

Analysis for AS-PCR. Melt Curve peak calls were determined using the CFX Maestro Software (Bio-Rad - 12004110), verified by eye, and exported to a customized C++ script to quickly convert melt curve peak calls to genotypes for each sample. Melt curve genotypes were then read into a customized R script (R Core Team 2019) for allele frequency determination and other statistical analysis.

本研究描述了秘鲁伊基托斯市18年间，埃及伊蚊（Aedes aegypti）的击倒抗性（knockdown resistance, kdr）单倍型随拟除虫菊酯类杀虫剂使用的演化历程。基于长达18年、累计约10000份样本的采样时长与强度，本研究是目前为止针对城市内埃及伊蚊kdr种群遗传学最为全面的研究。我们通过鉴定该种群中的两次演化事件，证明了全市范围的程序化拟除虫菊酯喷洒与特定kdr单倍型频率升高之间存在直接关联。即便在杀虫剂使用频率较低的情况下，较高的选择系数也凸显了埃及伊蚊种群的演化速度之快。在我们对蚊虫及其他昆虫害虫的相关文献调研中，尚未发现任何一种害虫能在如此少次暴露于低残留或无残留杀虫剂施用的情况下，以如此快的速度演化出抗性。本次观测到的抗性等位基因频率快速上升，可能得益于抗性相关等位基因相较于敏感等位基因的不完全显性。除了显著的时间维度变化外，空间抑制实验还揭示，遗传异质性不仅存在于全市尺度，也存在于城市内部的极小尺度范围内。

### 方法
#### 昆虫学调查
##### 防控历史
在过去二十年中，研究人员将整个伊基托斯市采集的蚊虫标本保存于标本库中（如Cavany等人，2020；Cromwell等人，2017；Getis等人，2003；Gunning等人，2018；LaCon等人，2014；Lenhart等人，2020；Morrison等人，2004a；Morrison等人，2004b；Morrison等人，2006；Morrison等人，2008；Reiner等人，2019；Schneider等人，2004；Tun-Lin等人，2009）。该标本库涵盖了拟除虫菊酯施用前（2000–2002年）、全市使用拟除虫菊酯期间（2002–2014年）以及秘鲁卫生部停用拟除虫菊酯后（2014–2017年）采集的样本。2002年之前，当地未开展针对埃及伊蚊的全市范围杀虫剂喷洒防控工作。2002年至2014年间，伊基托斯市卫生部开始针对埃及伊蚊开展防控，在居民住宅内喷洒多类拟除虫菊酯类杀虫剂，包括溴氰菊酯、氯氰菊酯、高效氯氰菊酯、氯氟氰菊酯以及高效氯氰菊酯+吡丙醚（补充表1）。居民自行在家中喷洒杀虫剂的情况较为罕见（A. Morrison, 个人通信）。2014年，由于出现表型拟除虫菊酯抗性，卫生部停用拟除虫菊酯，改用马拉硫磷（Gunning等人，2018）。在拟除虫菊酯施用期间，年均施药次数约为3.25次，多数年份的施药活动集中在一个月内。由于采用的是超低容量（ultra low volume, ULV）空间喷雾且无残留效果，当年大多数世代的蚊虫均未接触杀虫剂，且蚊虫样本的采集与施药活动无关（Ritchie等人，2021）。

##### 时间序列采样
自20世纪90年代末以来，美国海军医学研究所第6分队（NAMRU-6）与加州大学戴维斯分校的工作人员已将采集到的埃及伊蚊保存于-80℃环境中。本研究可获取2000年以来的标本。2009年6月前，蚊虫采集采用背负式吸虫器（Clark等人，1994）；2009年6月后，改用Prokopack吸虫器（Vasquez-Prokopec等人，2009；Reiner等人，2019）。标本库中的每只蚊虫均已鉴定至物种、性别、采集日期与采集地点。每个采集点通常为单个住宅，且已关联GPS坐标（图1）。

##### 空间采样
2013年与2014年，我们开展了基于拟除虫菊酯喷洒的高强度抑制实验（Gunning等人，2018），以验证针对埃及伊蚊种群动态的详细模型预测结果（Magori等人，2009）。简言之，我们选定了两个埃及伊蚊密度较高的城区，将其空间布局设置为：中央为施药试验区，外侧设置缓冲区作为实验对照区域（图2）。为限制迁移对抗性等位基因频率的影响，实验区域的尺寸被设定为埃及伊蚊预期终身飞行距离（约150米，Harrington等人，2005）的3至5倍。

2013年的实验区域面积约为750米×450米，包含1163户住宅。基线样本采集于2013年1月。系统采样始于2013年4月22日，持续16周至2013年8月8日。2013年4月29日至6月3日期间，在施药试验区每周开展6次无残留室内超低容量氯氰菊酯处理。

2014年的实验区域面积更大，约为600米×600米，包含2166户住宅。系统采样持续了44周。2014年4月28日至6月2日期间，参照2013年的实验方案，在试验区开展了氯氰菊酯超低容量喷雾。除本次实验施药外，2014年2月因登革热暴发开展了全市范围的喷洒作业，试验区与缓冲区的住宅均被喷洒了拟除虫菊酯。两次抑制实验期间，蚊虫样本的采集与保存方式与前文所述一致。

#### DNA提取与定量
本研究全程选择雄性蚊虫进行遗传分析，原因在于雌性蚊虫通常为病毒学与流行病学研究的对象，因此标本库中雄性样本数量更多。若使用雌性样本，则可能面临雄蚊虫（通过授精）以及人类（通过吸血）带来的基因组污染风险。由于电压门控钠离子通道（voltage-gated sodium channel, VGSC）基因不伴性遗传，雄性与雌性的等位基因频率预期一致。

完整的蚊虫样本经秘鲁与美国当局许可后，从秘鲁伊基托斯转运至美国北卡罗来纳州罗利市。样本在基因组DNA（genomic DNA, gDNA）提取前保存于-80℃，提取后保存于-20℃。基因组DNA从雄性埃及伊蚊全虫中提取，采用两种方法之一：Qiagen DNeasy血液与组织试剂盒（货号：69582），或加拿大DNA条形码中心的实验方案。简言之，对于Qiagen DNeasy试剂盒方案，将雄性全虫匀浆后，于55℃裂解缓冲液与蛋白酶K中孵育过夜。孵育并去除几丁质物质后，两种提取方法均采用RNase A处理以去除RNA污染。随后遵循Qiagen标准洗涤流程。最终样本用150 μl温热无菌去离子水（Invitrogen货号：10977-015）洗脱两次。为降低成本同时保证基因组DNA提取质量，部分蚊虫DNA提取采用了改良的加拿大DNA条形码中心（2020）实验方案。该方案中，样本的匀浆、孵育与RNase A处理流程与前文一致，随后将裂解液通过AcroPrep™ PALL2过滤板（货号：PALL 5053）以结合基因组DNA。过滤板先用蛋白洗涤缓冲液洗涤以去除残留蛋白，再用冷洗涤缓冲液洗涤以去除其他污染物。随后让过滤板自然晾干，确保无乙醇残留以避免影响DNA得率。最后用75 μl温热无菌去离子水洗涤两次，最终获得175 μl基因组DNA洗脱液。

基因组DNA的定量采用Quant-iT PicoGreen双链DNA检测试剂盒（ThermoFisher Scientific - P11496），样本检测在北卡罗来纳州立大学基因组科学实验室（GSL）的Synergy H1多功能酶标仪（BioTek Instruments, Inc.）上完成。

#### 基因分型
本研究采用等位基因特异性定量PCR结合熔解曲线分析（allele-specific quantitative PCR and melting curve analysis, AS-PCR），对中南美洲最常见的两个突变位点V1016I与F1534C进行重复两次的基因分型。若两次反应的分型结果不一致，则该样本被排除后续分析。分型不匹配的情况较为罕见，通常源于样本未成功扩增，或未满足分型判定标准（例如熔解峰未达到阈值）。另有少量样本（n=92）额外在V410L位点进行基因分型，以验证此前报道的该位点与V1016I位点之间存在强连锁不平衡（linkage disequilibrium, LD）（Saavedra-Rodriguez等人，2018）。所有在V410L位点分型的样本同样进行了两次重复实验以保证准确性。

##### V1016I位点基因分型
V1016I位点的AS-PCR方法基于Saavedra-Rodriguez等人（2007）报道的方案，并经美国亚特兰大疾控中心（CDC）昆虫学分支对I1016R引物的修改（A. Lenhart, 个人通信）。每个反应的PCR体系体积缩减至10 μl，包含2.5 μl无菌去离子水、0.5 μl 10 μM的各引物（V1016F、I1016F、I1016R）、5 μl PerfeCTa SYBR Green Supermix（Quanta – 95054-02K）以及1 μl模板。V1016F与I1016F的引物序列见Saavedra-Rodriquez等人（2007），但I1016R的引物序列已修改为：5’ - TGA TGA ACC SGA ATT GGA CAA AAG C – 3’（CDC, 个人通信）。样本在GSL实验室的BioRad CFX384实时PCR仪上进行分型，热循环程序如下：步骤1：95℃孵育3分钟；步骤2：95℃孵育10秒；步骤3：60℃孵育10秒；步骤4：72℃孵育10秒；步骤5：重复步骤2至4，共39个循环；步骤6：95℃孵育10秒；步骤7：熔解曲线分析，65℃至95℃，每10秒递增0.2℃并读取平板荧光信号。

##### F1534C位点基因分型
F1534C位点的AS-PCR实验参照Yanola等人（2011）报道的方法，并做如下修改：每个反应的总体系体积缩减至10 μl，包含5 μl无菌去离子水、0.2 μl C1534F引物、0.4 μl F1534F引物、0.4 μl F1534R引物、3.0 μl PerfeCTa SYBR Green Supermix（Quanta – 95054-02K）以及1 μl模板。样本在GSL实验室的BioRad CFX384实时PCR仪上进行分型，热循环程序如下：步骤1：95℃孵育2分钟；步骤2：95℃孵育30秒；步骤3：60℃孵育30秒；步骤4：72℃孵育30秒；步骤5：重复步骤2至4，共34个循环；步骤6：72℃孵育2分钟；步骤7：熔解曲线分析，65℃至95℃，每10秒递增0.2℃并读取平板荧光信号。

##### V410L位点基因分型
V410L位点的AS-PCR方法基于K. Saavedra开发的方案，并由CDC昆虫学分支共享（A. Lenhart, 个人通信）。每个反应的总体系体积缩减至10 μl：3.8 μl无菌去离子水、0.05 μl 50 μM Val410引物（5’ – GCG GGC AGG GCG GCG GGG GCG GGG CCA TCT TCT TGG GTT CGT TCT ACC GTG – 3’）、0.05 μl 50 μM Leu410引物（5’ – GCG GGC ATC TTC TTG GGT TCG TTC TAC CAT T – 3’）、0.1 μl 50 μM Rev410引物（5’ – TTC TTC CTC GGC GGC CTC TT – 3’）、5.0 μl PerfeCTa SYBR Green SuperMix（Quanta – 95054-02K）以及1 μl模板。热循环程序在GSL实验室的BioRad CFX384实时PCR仪上完成，具体如下：95℃孵育3分钟，40个循环（95℃孵育10秒，60℃孵育10秒，72℃孵育30秒），95℃孵育10秒，熔解曲线分析65℃至95℃，每10秒递增0.2℃并读取平板荧光信号。

##### AS-PCR数据分析
熔解曲线的峰型判定采用CFX Maestro软件（Bio-Rad - 12004110），并经人工目视验证，随后导出至自定义C++脚本，以快速将熔解曲线峰型结果转换为各样本的基因型。基因型数据随后导入自定义R脚本（R Core Team 2019），用于等位基因频率计算与其他统计分析。