Mosquito collection

Adult Cx. quinquefasciatus are routinely collected by Harris County Public Health (HCPH) personnel from Houston, Texas, United States. The samples used in this study were collected four times per year in 2016 and 2017 from two operational areas, 415 and 802. Samples from a single time point were collected within a one week time frame. Mosquitoes were stored at −80 ℃ until processing. There were no mosquito control operations by HCPH in area 415 during 2016 and 2017. In contrast, 26 malathion spraying operations were performed on ten separate dates, and ten permethrin spraying operations were performed on five separate dates during the summer of 2016 at area 802. Neither site was subjected to adulticidal spraying in 2017. From area 415, mosquito samples collected in weeks 14, 23, 38, and 46 in 2016 and in weeks 2, 20, 32, and 42 in 2017 were included in this study. From area 802, mosquito samples collected in weeks 14, 23, 39, and 48 in 2016 and in weeks 12, 20, 31, and 42 in 2017 were included.

DNA extraction

For each timepoint and location, 95 individuals were included. Total DNA was extracted and purified from each individual mosquito. The whole mosquito body was homogenized using a Qiagen® TissueLyser (Qiagen, Germantown, Maryland). Mosquito tissues were individually transferred into a well of a 96-well plate, which included a negative control. DNA extraction and purification were performed on the BioSprint® 96 workstation (Qiagen, Hilden, Germany) following their standard protocol. The final DNA product was stored at 4 ℃.

Species diagnostic

Morphological species identification was conducted by personnel from HCPH as the primary filtering for the field-collected adult Cx. quinquefasciatus. To confirm the result of morphological identification, species diagnostic polymerase chain reaction (PCR) was performed on each sample (Crabtree et al., 1995), using the PQ10 and the CP16 primers. Each 20 µL reaction included 10 µL 2× Thermo Scientific™ PCR Master Mix (Thermo Fisher Scientific, Carlsbad, California), 10 ng PQ10 primer, 10 ng CP16 primer, 20 ng template DNA, and nuclease-free water. The PCR program was set up following protocol by Crabtree et al. (1995) and was performed on the Eppendorf™ Mastercycler X50a thermocycler (Eppendorf, Hamburg, Germany). PCR products were visualized on a 2% agarose gel. Those samples that failed to show a clear 698-bp band were not used for microsatellite genotyping and data analysis in the following steps.

Microsatellite Genotyping

Eighteen microsatellite loci for Cx. quinquefasciatus (Fonseca et al., 2004, Smith et al., 2005, Edillo et al., 2007, Hickner et al., 2010) were amplified in multiplex PCR reactions. The forward primers were labeled using one of three fluorescent dyes, FAM, HEX and NED. Primers for six loci were multiplexed in a single reaction, with each primer in equimolar concentrations of 0.2 µM. Each 50 µL amplification reaction included 25 µL 2× QIAGEN® Multiplex PCR Master Mix, 5 µL 10× Primer Mix, 20ng DNA and nuclease-free water (with 3mM MgCl₂ in final concentration). The thermoprofile for the multiplex PCR was: 2 min at 95 ℃, followed by 30 cycles performed for 30 sec at 95 ℃; 90 sec at 57 ℃; and 30 sec at 72 ℃, and a final cycle for 30 min at 72 ℃. The PCR products were submitted to the DNA Analysis Facility on Science Hill at Yale University for fragment analyses. Genotyping was performed using Genemarker (Hulce et al., 2011).

Data Analysis

Microchecker version 2.2.3 was used to check for microsatellite null alleles (Van Oosterhout et al., 2004). Loci CX4, CX10, and CX11 were not included in subsequent analyses, as null alleles were detected in these three loci. Due to low amplification efficiency, data of CX5 from area 802 was also not included in downstream analysis. Microsatellite genotyping data were converted to Genepop format with GenAlEx version 6.503 (Peakall and Smouse, 2012).

Data were loaded into DIYabc version 2.1.0 (Cornuet et al., 2014). This software package implements approximate Bayesian computation inferences about population history (Cornuet et al., 2014). We tested four basic demographic scenarios in our study, which included variations of population size (bottleneck, expansion, decline, or constant)for two time periods. For each scenario, 2,000,000 simulated datasets were generated, with specific mutation models and prior distributions for parameters. The summary statistics of each simulated dataset were compared with those of the observed dataset. More specifically, we used three single-sample statistics, including the mean number of alleles, mean genetic diversity (Nei, 1987), mean allele size variance across loci, and three two-sample statistics, including the fixation index (F_ST) (Weir and Cockerham, 1984), mean index of classification (Rannala and Mountain, 1997, Pascual et al., 2007), and Euclidean distance between every two samples (Goldstein et al., 1995). Through this process, the posterior probability of each scenario was calculated, and the most likely scenario was identified. The effective population size of Cx. quinquefasciatus from each site of the Houston area was obtained by estimating the posterior distribution for the most likely scenario.

We set two weeks as the average generation time of Cx. quinquefasciatus, as the life cycle from egg to adult stage usually takes 10 to 14 days (Manimegalai and Sukanya, 2014). The minimum value for effective population size (N_e) was set as 0 and the maximum value for 1,000,000, except for the bottleneck scenario. The prior range for N_e at earlier time points was between 0 and 1,000,000, while the prior range of Ne at the latter time points was between 0 and 200,000 for the bottleneck scenario. 

We tested three time points from spring in 2016 to winter in 2017 per area to find the effects of the winter season. Genotyping data from weeks 14 and 23 in 2016 and week 2 in 2017 were used for area 415, while genotyping data from weeks 14 and 23 in 2016 and week 12 in 2017 were used for analysis for area 802. First, we performed the parameter estimation with a prior normal distribution, where parameters had a mean value and a standard deviation within a range between extremum values (minimum and maximum). We also performed the parameter estimation with a uniform prior distribution to test if different prior distributions affect the choice of best demographic scenario. In this case, the prior distribution had minimum and maximum values but without the mean value and standard deviation.

Next, we performed N_e estimation from three time points from 2016 and three timepoints from 2017 at area 802 to evaluate the effects of Hurricane Harvey on the N_e of the mosquito population. The mosquito samples were collected before and after Hurricane Harvey, which landed in the Houston area on August 27, 2017. For area 802, genotyping data from weeks 14, 23, and 48 in 2016 were used for pre-hurricane analysis in DIYabc, while genotyping data from weeks 12, 20, and 42 in 2017 were used for post-hurricane analysis.

In addition to DIYabc, we used NeEstimator version 2.1 (Do et al., 2014) to estimate the N_e of Cx. quinquefasciatus from eight sampling time points in 2016 and 2017 at area 802 to further explore the impact of Hurricane Harvey on the N_e of mosquito population in Harris County. The N_e at each time point was estimated with the Linkage Disequilibrium (LD) method. The mating system was defined as random mating.