Collections of Aedes aegypti

Highland and lowland study sites for each of the Philippine big islands (Luzon, Visayas, and Mindanao) were established (Fig. 1). These highland and lowland sites included, respectively, Baguio city (BG) (16.402°N; 120.596°E) and Quezon city (QC) (14.676°N; 121.044°E) in Luzon; Cebu city (CC) mountains (10.317°N; 123.891°E) and Liloan (LIL), Cebu (10.4121°N; 123.986°E) in the Visayas; and Maramag, Bukidnon (BUK) (8.454°N; 124.632°E) and Cagayan de Oro (CDO) city (8.051°N; 124.923°E) in Mindanao. These sites were selected based on: 1) the same biotype using Corona's modified climate classification (Corporal-Lodangco and Leslie 2017), 2) high incidence of dengue cases (NDRRMC 2019, WHO WPRO 2019), 3) elevation, and 4) range of latitudes.

Ae. aegypti parental eggs were collected from natural populations in the study sites, preferably outside human dwellings, under the roof gutter with the use of modified ovicidal/larvicidal (O/L) traps. Permission from the Department of Health—Bureau of Quarantine Office, Cebu city for shipment of Aedes eggs was obtained beforehand. Aedes eggs were hatched inside transparent plastic cups (TPCs) half-filled with ascorbic acid solution (Lacour et al. 2015). Larvae were carefully pipetted into new TPCs half-filled with distilled water (DW) (∼150 mL) and covered with fine-mesh cloth. Daily feeding of fish food (0.02 g) (Sakura, All Aquariums Co., Ltd., Thailand) and changing of DW every other day to avoid scum formation on the water surface were done until the larvae metamorphosed into adults. Newly emerged adults of Ae. aegypti were identified (Belkin 1962) and sorted out by sex, and were stored into a 1.5 mL microcentrifuge tube individually for DNA extraction.

To test for Mendelian inheritance, first generation (F1) offspring of Ae. aegypti were obtained by coupling the male and female reared adults from parental eggs collected from CC and BUK sites. These adult couples were put inside a TPC covered with fine-mesh cloth and with a filter paper at the bottom. They were fed with a 10% sucrose solution soaked in a cotton ball that was placed on top of the fine-mesh cloth. A 3- to 5-d mating window was given for each couple. After which, females were starved for 24 h, and were fed with chicken blood mixed with ethylenediaminetetraacetic acid (EDTA) placed on the bottom depression of plastic bottles and sealed with a thin parafilm membrane. Each plastic bottle was placed in an upright position on top of the fine-mesh cloth of PTC with coupled mosquitoes, and was frequently added with lukewarm water to simulate body heat to attract female mosquitoes to suck blood. The mosquitoes were then fed with 10% sucrose solution soaked in cotton balls that were put on top of the fine-mesh cloth of PTC. The filter paper at the bottom of the mosquito cup was moistened with DW for females to lay F1 eggs relatively 2–5 d later. F1 eggs were then hatched and reared until adult emergence.

DNA Extraction

DNA samples (n = 30 each site) were extracted from Ae. aegypti natural populations from BG, QC, LIL, and CDO city, and F1 siblings from CC highlands and BUK using Macherey-Nagel insect DNA extraction kit (Cat. No. 740470.50; Macherey-Nagel, Düren, Germany) following the manufacturer's protocol in a level 2 bio-safety cabinet.

Selection of Isolated Microsatellite Loci of Aedes aegypti

Briefly, a total of 327,094 perfect SSRs were searched in the updated reference genome assembly of Ae. aegypti Liverpool strain (Matthews et al. 2018) using the MISA tool (Beier et al. 2017). Only 52,002 SSRs were located in chromosomes 1–3. The annotated 1,513 SSRs were filtered based on motif length (n = 2–4 bp), equal or less than 20 repeats, and located in the intergenic region. Primer BLAST was performed to generate candidate primers at the flanking sequences of the SSRs. Out of the 1,513 sequences, 21 SSRs with di-, tri-, and tetra-nucleotide motifs from each of the three Ae. aegypti chromosomes (total = 63 SSRs) were randomly selected for characterization considering also the logistics of this study. Their primer sequences were sent to International DNA Technologies (IDT) Inc., Singapore for synthesis.

Test on Polymorphic Microsatellite Loci

Successful amplification of the 63 SSR markers in the Philippine populations of Ae. aegypti was crucial in determining their presence in the genome of Ae. aegypti local populations. These loci were amplified using pooled DNA samples of male and female Ae. aegypti from each site (n = 12). An M13 sequence (5'- TCC CAG TCA CGA -3') was added to the 5' end of each forward (F) primer of the nine SSR markers. PCR reactions (final volume of 10 µL (containing 5µL of 2x Type-it multiplex PCR master mix [Qiagen Type-It microsatellite PCR kit; Cat. No. 206243], 1 µL each of 10 µM forward and reverse primers, 1 µL pooled DNA sample, 2.25 µL nuclease-free water [NFW]) were run on a Veriti 96-well Thermal Cycler (Cat. No. 4375786: Thermo Fisher Scientific, Applied Biosystems, MA, USA) with 94°C initial denaturation for 10 min, then 35 cycles of 94°C x 30 s, 55°C x 40 s, 72°C x 30 s and final extension for 5 min at 72°C (Slotman et al. 2007). Gel electrophoresis was performed using 1% agarose to determine the quality of microsatellite bands. A mixture of 2 µL loading dye and 2 µL PCR product of each candidate SSR maker was loaded into each well of an agarose gel together with AC4 locus (Slotman et al. 2007) as positive control and NFW as negative control. Bands in agarose gels were viewed through a UV trans-illuminator (Maestrogen, Nevada).

Test for Broadly Applicable Microsatellite Markers

To determine the broad applicability of these molecular markers, outgroups of Ae. aegypti samples from Rayong (Ry), Thailand and Ho Chi Minh (HCM), Vietnam were assessed using previously screened polymorphic microsatellite markers. Pooled DNA samples (n = 30) from each outgroup were amplified following the same concentrations of reagents and PCR profile described earlier. Amplified products were visualized on agarose gels through the UV trans-illuminator (Maestrogen, Nevada). Markers with smeared and faint bands were eliminated.

Characterization of Polymorphic Microsatellite Loci of Aedes aegypti

All 180 DNA extracts of Philippine Ae. aegypti samples together with 60 DNA samples (n = 30 each site) from Ry, Thailand and HCM, Vietnam were amplified following the PCR profile such as denaturation step at 94°C for 5 min, then 38 cycles of 94°C x 30 s, 54°C x 30 s, 72°C x 20 s and final extension for 10 min at 72°C in a Veriti thermal cycler (Cat no. 4375786; Applied Biosystems). Several optimizations (i.e., from 52°C up to 68°C with an increment of 2°C) were performed to determine the optimal annealing temperature (54 °C) of all SSR markers. PCR reactions with final volume of 10 µL (containing 5 µL of 2x Type-it multiplex PCR master mix [Qiagen Type-It microsatellite PCR kit; Cat. No. 206243], 0.2 µL of 10 µM forward [F] primer, 0.8 µL each of 10 µM reverse [R] primer and fluorescent M13, 1 µL pooled DNA sample, 2.2 µL NFW) were prepared. To get maximum peak signals, two important points were observed during the preparation of PCR reaction mixture (Schuelke 2000): 1) equimolar amounts of R primer and fluorescent-labeled M13 primer should be added to the reaction mix, and 2) the concentration of the F primer should be one-fourth that of the concentrations of R primer and M13 primer to facilitate the M13 primer-driven amplification reactions when the F primer was consumed.

PCR products were sent to the Philippine Genome Center (PGC), QC for fragment analyses by using ABI 3730xl DNA Genetic Analyzer (Cat no. 3730XL; Applied Biosystems) (Carlsbad, CA) and GS 500 Rox allele size standard (Cat no. 4340060A; Applied Biosystems). The two highest allele peaks of each genotype were scored using ABI peak scanner software (v.1) (Applera Corp., Norwalk, CT). The height of each peak on the y-axis was zoomed up to 28,000 relative fluorescence units (RFUs), while the allele size on the x-axis was zoomed at 90–300 bp. Capillary DNA sequencing was also performed to confirm the sequences of these microsatellite markers at the PGC, QC. Data sequences were analyzed using Chromatogram Explorer Lite of DNA Baser Assembler software (Heracle BioSoft S.R.L., Romania).

Genetic Analyses

MICROCHECKER software (v. 2.2.3) (Oosterhout et al. 2004) was used to determine the presence of null alleles. For Mendelian codominant inheritance, goodness-of-fit for Hardy-Weinberg equilibrium (H-WE) and chi-square test (P < 0.001) using GENEPOP software (v. 4.7.5) (Raymond and Rousset 1995) for genotypes of F1 offspring samples were performed. Other GENEPOP tests were done such as number of alleles, allele frequencies, estimates of observed (H_O) and expected heterozygosities (H_E), goodness-of-fit probability tests for H-WE to determine if the markers were selectively neutral, Fisher's exact test to detect linkage disequilibrium (LD), and pairwise estimates of Fisher's statistics (F_ST) between populations to assess the robustness of the SSR markers only. Sequential Bonferroni multiple comparison tests for all loci were done in H-WE and LD tests. For the latter, Markov chain parameters were set at 1,000 dememorizations, 100 batches, and 1,000 iterations. Corrections were made on the number of batches and changes were made in increments of 100.

Polymorphic Microsatellite Markers

Twenty-six out of 63 primer pairs of SSR markers consistently amplified using pooled Ae. aegypti DNA extracts (n = 12) from different Philippine sampling sites, as depicted on electrophoretic gel bands. Seven out of the 26 loci were located in the first chromosome, nine in the second chromosome, and ten in the third chromosome. Three motifs were observed, namely, dinucleotide (n = 10), trinucleotide (n = 12) and tetra-nucleotide (n = 4). These motifs were consistent in previous studies (Slotman et al. 2007, Chambers et al. 2007, Brown et al. 2011).

Table 1.

Characterization of novel and broadly applicable microsatellite markers of Ae. aegypti using Philippine and outgroup populations

Broadly Applicable Microsatellite Markers

Nine out of the 26 SSRs of Ae. aegypti showed promising bands on pooled DNA extracts of outgroup samples, suggesting their consistent amplification, hence, are good candidates of broadly useful microsatellite markers. Faint and smeared bands of PCR products of 26 polymorphic SSR markers that amplified samples from both outgroups (Thailand and Vietnam) on an agarose gel indicated unsuccessful amplification and were eliminated.

Screening for Null Alleles and Mendelian Inheritance

MICRO-CHECKER software (v. 2.2.3) did not detect significant evidence of null alleles in four (Aaeg1-3D, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C) out of the nine SSR markers among genotyped samples from local populations (Table 1;  Supplemental file 1). GENEPOP (v. 4.7.5) analysis in screening for Mendelian codominant pattern of inheritance revealed that for F1 Ae. aegypti belonging to CC family, four markers (Aaeg1-3D, Aaeg1-3H, Aaeg2-2E, and Aaeg3-4C) adhered to Mendelian codominant pattern (X² test, df = 18, P > 0.0011) after incorporating sequential Bonferroni correction for multiple comparison. For BUK F1 mosquito family, four markers (Aaeg1-3D, Aaeg1-3H, Aaeg2-2A, and Aaeg3-4C) also conformed to Mendelian pattern of inheritance (X² test, df = 18, P > 0.0011).

Screening for H-WE, LD, and F_ST Estimates

After sequential Bonferroni correction for multiple comparison, Table 1 revealed that overall, four SSR markers (Aaeg1-3D, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C) conformed to H-WE (P > 0.0011). Aedes aegypti populations in BG and QC conformed to H-WE (P > 0.0011) in five loci (Aaeg1-3D, Aaeg1-3H, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C). LIL mosquito population had six loci (Aaeg1-3D, Aaeg1-3H, Aaeg2-2A, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C) that conformed to H-WE (P > 0.0011). CDO population had five loci (Aaeg1-3D, Aaeg2-2A, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C) that adhered to H-WE (P > 0.0011). LD between two loci after sequential Bonferroni multiple comparison (P > 0.00028) in all mosquito populations were not observed. The 9-locus dataset produced similar F_ST estimates as the 7-locus dataset that included the two best markers (Aaeg1-3D and Aaeg3-4C) and the rest of the SSR markers (Aaeg1-3H, Aaeg2-2A, Aaeg2-2B, Aaeg2-4D, and Aaeg3-3I) but differed in one F_ST estimate only between Ry and BG population pairs based on the 7-locus dataset that included the two good markers but not co-dominant (Aaeg2-2E and Aaeg3-2A) and the rest of the SSR markers mentioned.

Table 2.

Primer sequences (F—forward, R—reverse), chromosome position based on the reference genome (Matthews et al. 2018), allele size range, total number of alleles observed, and mean H_E and H_O of the nine polymorphic SSR markers of Philippine Aedes aegypti (Sample size of all loci = 240)

Polymorphism and Heterozygosities (H_E and H_O)

Table 2 shows that nine loci were polymorphic although they varied in the number of alleles per locus (2–20 alleles) in each population. These sequence data were submitted to the GenBank database with their accession numbers (Table 2). In all populations, Aaeg3-3I locus was the most polymorphic with 37 alleles; Aaeg2-4D, the least polymorphic with only three alleles. H_E estimates were moderate to high (0.41–0.94) (Table 2). No singletons were detected in Aaeg2-2E, Aaeg2-4D, and Aaeg3-4C, whereas one singleton in Aaeg1-3D, Aaeg1-3H, and Aaeg2-2A, and ten in Aaeg3-3I.

Characteristics of Microsatellite Loci in Outgroup Populations

Polymorphism and Heterozygosities

All nine loci were polymorphic. H_E estimates were moderate (0.46) to high (0.89) in populations of Ae. aegypti from the Philippines and outgroups (Ry, Thailand and HCM, Vietnam). Excess of homo-zygotes (or heterozygote deficiencies) were observed in Aaeg1-3H, Aaeg2-2A, Aaeg2-2B, Aaeg2-4D, and Aaeg3-3I, implying that homo-zygotes were greater than expected under H-WE (Selkoe and Toonen 2006). Several reasons include the presence of null alleles (Callen et al. 1993, Ardren et al. 1999), effects of inbreeding (Callen et al. 1993, Selkoe and Toonen 2006), and selection for or against a specific allele (Selkoe and Toonen 2006).

Null Alleles and Mendelian Inheritance

Absence of null alleles were consistently observed in four markers (Aaeg1-3D, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C) in both Ae. aegypti populations from local and outgroups. Null alleles are caused by loss-of-function mutation (Chapuis and Estoup 2006) and may cause biases in estimating genotypic and allelic frequencies in population genetic studies (Oosterhout et al. 2006), heterozygote deficiencies (Callen et al. 1993, Ardren et al. 1999), errors in genetic parentage analysis (Dakin and Avise, 2004), and misleading analysis in population differentiation (Chapuis and Estoup 2006). Thus, it is crucial for a microsatellite marker to lack null alleles.

On the other hand, three markers (Aaeg1-3D, Aaeg1-3H, and Aaeg3-4C) consistently conformed to Mendelian codominant pattern in both F1 mosquito families from CC and BUK, and may be useful in assessing genetic connectivity of Ae. aegypti populations. Departures from Mendelian fashion may be due to nondisjunction during meiosis, presence of null alleles (Oosterhout et al. 2006), and mutations during sexual reproduction (Tesson et al. 2013).

Table 3.

Goodness-of-fit probability tests for Hardy-Weinberg Equilibrium of the nine microsatellite markers in local and outgroup populations of Aedes aegypti. P-value = 0.0011 incorporating sequential Bonferroni correction for multiple comparison tests; (*) deviated from H-WE

Fig. 2.

Location of the previously reported and the new SSR loci in the genome of Aedes aegypti using the published genome assembly (Matthews et al. 2018) as the reference. Each location is marked by the name of the locus and its coordinate (in bp), in blue font for published markers (Slotman et al. [2007], Brown et al. [2011] and Lovin et al. [2009]) or in red font for the new markers (this study). The coordinates indicate the location of the first 5' nucleotide position of the forward (F) PCR primer for the published markers or the first 5' nucleotide position of the repeat region for the new markers (from Table 2, column 3).

Selectively Neutral and Randomly Associated SSRs

Four markers (Aaeg1-3D, Aaeg2-2E, Aaeg3-2A, and Aaeg3-4C) consistently conformed to H-WE in both local and outgroup populations of Ae. aegypti (Table 3). Conformity to H-WE implies that allele frequencies in a population with random mating remain constant from one generation to the next in a large population in the absence of evolutionary forces such as mutation, migration, and natural selection. Violations on any or all of these assumptions cause deviations from H-WE (Hartl 2020). Moreover, deficiency of heterozygotes generated by null alleles may also cause deviations from H-WE (Lanzaro et al. 1998, Lehmann et al. 2003). Moreover, random association between two loci were consistently observed in both local and outgroup populations, suggesting close physical genetic linkage was unlikely. Factors that influence LD include selection, mating system, genetic drift, mutation rate, population structure (Hartl 2020), inbreeding, and null alleles (Callen et al. 1993, Lehmann et al. 1998). LD is not expected if heterozygote deficiency is caused by null alleles (Lehmann et al. 2003).

Several studies on the population structure of Ae. aegypti (Brown et al. 2011, Fonzi et al. 2015, Sayson et al. 2015, Hopperstad et al. 2019) used the seven markers developed by Chambers et al. (2007), 13 markers by Slotman et al. (2007), four markers by Brown et al. (2011), and four re-designed markers by Lovin et al. (2009). Slotman et al. (2007) and Brown et al. (2011) characterized these markers by assessing deviations from H-WE, screening for null alleles, and LD, although Brown et al. (2011) utilized their four new markers with eight previously developed markers to assess the population structure of Ae. aegypti worldwide. Moreover, Chambers et al. (2007) included LD, classification of microsatellite repeats, and polymorphism. Lovin et al. (2009) assessed polymorphism, deviation from H-WE, genetic map location, and multiplex-ready groups of markers that were used to determine the population genetics of Ae. aegypti populations in Haiti. However, these studies (Chambers et al. 2007, Slotman et al. 2007, Lovin et al. 2009, Brown et al. 2011) did not determine the codominant Mendelian inheritance. Only Lovin et al. (2009) and this current study provided the chromosomal location of microsatellite markers based on the linkage map positions of the previously defined genetic loci, and on Ae. aegypti reference genome (Matthews et al. 2018), respectively. Moreover, we located the chromosomal locations of the current SSR markers relative to the published markers (Chambers et al. 2007, Slotman et al. 2007, Lovin et al. 2009, Brown et al. 2011) using the reference genome (Matthews et al. 2018) (Fig. 2). Thus, this current study characterized comprehensively novel SSR markers utilizing natural populations of Philippine Ae. aegypti, for the first time, by analyzing H-WE, LD, codominant Mendelian inheritance in F1 siblings, screening for null alleles, heterozygosities, chromosomal locations relative to other markers and broad applicability in other Asian countries.

Isolation of polymorphic microsatellite markers in Culex pipiens (Keyghobadi et al. 2004) has also been reported with H_E and H_O analysis only, which was utilized in a separate study by Keyghobadi et al. (2006) to determine the population structure of Cx. pipiens in Hawaii. Edillo et al. (2007) characterized microsatellite markers of Cx. pipiens by analyzing deviations from H-WE, LD, codominant Mendelian inheritance, heterozygosities (H_E and H_O), presence of null alleles, and broader applicability in other countries. Edillo et al. (2009) utilized these good markers (Edillo et al. 2007) in a separate study on the population structure of Cx. pipiens populations in eastern North America. In Anopheles (Lather et al. 2015), genetic analysis of polymorphic microsatellite markers included deviations from H-WE, LD, codominant Mendelian pattern, presence of null alleles, heterozygosities, and cross-amplification in other species. In other taxa, Masran and Majid (2018) developed seven microsatellite markers from Cimex hemipterus (Fabricius) (Hemiptera: Cimicidae) by analyzing results from H-WE, LD, and screening for null alleles on their natural populations. Recent SSR characterization isolated from fin whales (Tardy et al. 2020) and earthworm (Liu et al. 2020) assessed H-WE, LD, allelic diversity, heterozygosities, and cross-amplification with other species. These studies (Keyghobadi et al. 2004, Edillo et al. 2007, Slotman et al. 2007, Lather et al. 2015, Masran and Majid 2018, Liu et al. 2020, Tardy et al. 2020) did not necessarily compare between newly characterized microsatellite markers with a set of previously good published markers and those datasets of published markers with problematic loci like that of Brown et al. (2011) for Ae. aegypti.

Other techniques include genotyping-by-sequencing (GBS) (Elshire et al. 2011) and restriction site-associated DNA sequencing (RADseq) (Davey and Blaxter 2010). The latter has been used to develop markers of Ae. aegypti (Evans et al. 2015, Rasic et al. 2015). Sunde et al. (2020) noted that between RADSeq and microsatellite analysis, RADSeq is a better technique because it produces a higher number of loci.

In conclusion, this study determined two novel and broadly useful microsatellite markers of Philippine Ae. aegypti such as Aaeg1-3D of chromosome 1 and Aaeg3-4C of chromosome 3, that can be added to markers of this species previously developed. Moreover, two markers (Aaeg2-2E of chromosome 2 and Aaeg3-2A of chromosome 3) may still be of use since they met all criteria except that they did not conform to Mendelian codominant inheritance. These latter markers can be used with caution, i.e., these should be utilized together with Aaeg1-3D and Aaeg3-4C and other good microsatellite markers (Chambers et al. 2007, Slotman et al. 2007, Lovin et al. 2009, Brown et al. 2011). In the Philippines, since dengue outbreaks are recurring (NDRRMC 2019, WHO WPRO 2019) and currently only a few studies on the genetic structure of Ae. aegypti were conducted (Fonzi et al. 2015, Sayson et al. 2015, Carvajal et al. 2020), the need to conduct more studies on population structure of this species in the country, in particular, and in the tropics, in general, is very relevant. The additional markers developed in this study will be very useful in studies on the population structure of this species in the tropics, as in our separate study of the Philippine Ae. aegypti to promote effective and timely vector control strategies for public health. In the future, these research initiatives may lead to implementation of urgent actions to safeguard public health against dengue, Zika, chikungunya, and yellow fever that are transmitted by this vector species.