Population genetics of Anopheles koliensis through Papua New Guinea: New cryptic species and landscape topography effects on genetic connectivity

Abstract New Guinea is a topographically and biogeographically complex region that supports unique endemic fauna. Studies describing the population connectivity of species through this region are scarce. We present a population and landscape genetic study on the endemic malaria‐transmitting mosquito, Anopheles koliensis (Owen). Using mitochondrial and nuclear sequence data, as well as microsatellites, we show the evidence of geographically discrete population structure within Papua New Guinea (PNG). We also confirm the existence of three rDNA ITS2 genotypes within this mosquito and assess reproductive isolation between individuals carrying different genotypes. Microsatellites reveal the clearest population structure and show four clear population units. Microsatellite markers also reveal probable reproductive isolation between sympatric populations in northern PNG with different ITS2 genotypes, suggesting that these populations may represent distinct cryptic species. Excluding individuals belonging to the newly identified putative cryptic species (ITS2 genotype 3), we modeled the genetic differences between A. koliensis populations through PNG as a function of terrain and find that dispersal is most likely along routes with low topographic relief. Overall, these results show that A. koliensis is made up of geographically and genetically discrete populations in Papua New Guinea with landscape topography being important in restricting dispersal.

. Because of this complexity, New Guinea harbors a diverse array of locally endemic species (taxa that are endemic to certain parts of the island; Allison, 2007). Previously identified patterns of biological connectivity reveal that many New Guinean species show strong affinities to both Australia and Asia. The population structure and local endemicity observed in some species suggest that the geological history of the region has been important in the diversification of species. In many cases, vicariance is thought to play a more important role than dispersal in the current distributions of species in New Guinea and parts of Australasia (Heads, 2013). This theory is supported by the presence of biological breaks that occur in locations associated with historically separate landmasses in New Guinea as well as to the north and the south of the Central Range. This mountain range traverses the island from east to west, presenting a northsouth barrier to dispersal for most species (Heads, 2013;Macqueen, Goldizen, Austin, & Seddon, 2011).
While the biogeography of New Guinea has been explored relatively extensively, there is less understood regarding the population genetics and phylogeography of individual species throughout the region. The Punctulatus Group of mosquitoes currently comprise 13 cryptic species that are widely distributed throughout New Guinea and northern Australia and show overlapping distributions. Some species of this group transmit malaria, and they provide a useful system with which to study the genetic diversity of organisms relative to the landscape features of this region. As larvae, they show highly specific aquatic requirements, however, as adults they can be on the wing and their ability to disperse is not well understood and may vary between species and sex (Ambrose et al., 2014). A study of two closely related isomorphic mosquito species in this group-Anopheles hinesorum Schmidt and the coastally restricted Anopheles farauti Laveran reveal similarities in the location of genetic breaks in New Guinea (Ambrose et al., 2012). In stark contrast, a third species in this group-Anopheles punctulatus showed little apparent population structure between these regions despite the use of fast-evolving microsatellite markers (Seah, Ambrose, Cooper, & Beebe, 2013), suggesting it has undergone a recent range bottleneck and expansion throughout New Guinea and has recently traversed the Central Range.
Here, we focus on a fourth cryptic member and important malaria vector of the Punctulatus Group, Anopheles koliensis (Owen; Cooper et al., 2009;Slooff, 1964). This species is found throughout New Guinea, New Britain, and until recently on several islands of the Solomon Archipelago Spencer, Spencer, & Venters, 1974;Taylor, 1975). In Papua New Guinea (PNG-eastern New Guinea), A. koliensis appears to be primarily a lowland species, inhabiting the river valley flood plains mostly below 300 m throughout the northwestern and southeastern lowland regions of PNG-it is rarely found in the distinct wet/dry monsoonal climate area of southwest PNG (Beebe, Russell, Burkot, & Cooper, 2015;. Surveys through PNG suggest oviposition occurs in natural ground pools and swamps, as well as in human-disturbed or human-modified sites, including vehicle wheel tracks and drains (Cooper, Waterson, Frances, Beebe, & Sweeney, 2002;Slooff, 1964).
This mosquito was also found to have a positive association with human habitation .
Due to both overlapping and variable adult morphology of the Punctulatus Group, molecular tools including genomic DNA probes and PCR are required for species identification (Beebe, Foley, Cooper, Bryan, & Saul, 1996;Beebe & Saul, 1995). The most utilized diagnostic tool is a PCR, restriction fragment length polymorphism analysis (RFLP) of the ribosomal DNA (rDNA) internal transcribed spacer 2 (ITS2). This diagnostic marker clearly discriminates A. koliensis from all other members in the Punctulatus Group based on agarose gel electrophoresis (Beebe & Saul, 1995).
An earlier study through northeast PNG (Madang and East Sepik region) using the PCR-RFLP species diagnostic tool with a sensitive acrylamide gel size separation suggested A. koliensis comprised three ITS2-RFLP variants (Benet et al., 2004). These ITS2 RFLP variants were designated genotype W, found only around Wosera in northeast PNG's Sepik region; genotype M, found only in the Madang region and MW, found throughout their study sites in northern PNG (Madang and Wosera regions). The authors suggested these genotypes may be evidence for distinct subspecies as some variation in time of night blood feeding was found (Benet et al., 2004). Variation in night feeding time may be important for the capacity of this species to develop behavioral resistance to the long-lasting insecticide-treated bed nets now deployed in the region (Russell, Beebe, Cooper, Lobo, & Burkot, 2013).
In this study, we further investigate the population genetics of A. koliensis in Papua New Guinea. We initially aim to verify the existence of the three previously described rDNA ITS2 RFLP variants within A. koliensis and assess reproductive isolation between individuals with these genotypes. We use additional nuclear loci (DNA sequence and microsatellite) as well as a mitochondrial locus to achieve this. We use this sequence and microsatellite data to describe the population genetic structure of A. koliensis throughout PNG and hypothesize that the complex geological history and landscape topography of New Guinea have shaped population structure in this mosquito. Finally, we use recently developed methods to assess whether landscape topography plays a role in maintaining this genetic structure, providing the first contribution to the literature on landscape genetics of a species in New Guinea.

| Mosquito collection, identification, and genotyping
Mosquitoes were collected by human-landing catches, CO 2 -baited light traps, or through larval collections with the larvae being bred out to adults. Using these sampling techniques, A. koliensis was collected throughout PNG (Tables 1 and 2, Figure 1) and stored frozen, in alcohol, or desiccated on silica gel. Genomic DNA was extracted using a salt extraction method (Beebe, Ellis, Cooper, & Saul, 1999), and samples were identified to species by a speciesspecific PCR-RFLP of the ITS2 (Beebe & Saul, 1995). Only samples identified as A. koliensis were analyzed further. To confirm the presence of the M, MW, and W genotypes that were identified by using acrylamide gel electrophoresis RFLP separation (Benet et al., 2004), we used a 3% agarose gel and slightly longer ITS2 primers outlined in Alquezar, Hemmerter, Cooper, and Beebe (2010) for improved resolution of the three genotypes (Alquezar et al., 2010).

| DNA sequencing and analysis of the ITS2, COI, and rpS9
Due to the multicopy nature of the rDNA, multiple variant copies (paralogues) of the ITS2 gene are often present within individual genomes and PCR products. These intragenomic sequence variants often cannot be directly sequenced using Sanger DNA sequencing as paralogues containing indels can collapse the chromatogram. For this reason, three individuals from sympatric populations from sites 1,330, 1,109, and 1,136 (see Table 1 for site localities) were PCR-amplified for the ITS2 and cloned into bacteria with 3-5 randomly selected colonies sequenced following the methods of Alquezar et al. (2010).
A 527-bp segment of the mtDNA COI and a 430-bp region of nuclear ribosomal protein S9 (rpS9) were also amplified and directly sequenced following methods of Ambrose et al. (2012). For the rpS9, pseudohaplotypes of the nuclear locus were inferred using the program PHASE implemented in DNAsp v5 (Librado & Rozas, 2009). Additionally, DNAsp v5 was used to estimate haplotype diversity and nucleotide diversity, as well as to test for neutrality using Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) for each locus. Recombination was assessed using the program RPD3 (Martin et al., 2010). Haplotype networks were constructed using TCS 1.21 (Clement, Posada, & Crandall, 2000) under a 95% connection limit. Genetic diversity for mtDNA and rpS9 was estimated in a number of ways including haplotype diversity, haplotype richness and evenness (Kimura, 1964), and nucleotide diversity. Haplotype evenness provides an estimate of the number of equally frequent alleles that would result in the estimated haplotype diversity of the sampled population. Diversity estimates were obtained using the R packages gstudio (Dyer, 2014), pegas (Paradis, 2010), and adegenet (Jombart, 2008).

| Microsatellite development and amplification
Primers for microsatellite analysis were obtained from 454 pyrosequencing data of the genomic DNA of A. koliensis following the methods of Ambrose et al. (2014). A total of 24 microsatellite markers were initially screened. Primers were selected based on their ability to consistently amplify a clean product (single locus) from all populations. No markers were excluded due to low variability, and details of primers used can be found in Table S1. Each locus was amplified by PCR using fluorescently labeled forward primers.

| Population parameters and statistics
Alleles for each marker were scored manually in the program between loci. To estimate population genetic structure between all pairs of sites (used subsequently for landscape genetic analyses), we used Hedrick's GST' (Hedrick, 2005), a measure which is not biased by within-population variation as implemented in MMOD (Winter, 2012).

| Population structure: Bayesian clustering
In order to assess broadscale population structure in New Guinea, we used a Bayesian clustering method as well as discriminant analy-

| Population structure: multivariate approach (DAPC)
To complement the Bayesian analysis run in STRUCTURE (outlined below), we also performed a discriminant analysis of principal components (DAPCs) using the R package adegenet (Jombart, 2008;Jombart et al., 2010). This is a multivariate analysis designed for the use on genetic data such as microsatellites. It does not require that specific population genetic assumptions be met or try to fit data to a predefined model. Advantages over principal component analysis (and principal coordinate analysis) include that it is optimized to maximize variation among rather than within groups. We used the online server (Jombart, 2008) to determine the optimal number of principal components (PCs) to retain. We then performed DAPC on the five groups identified using STRUCTURE, retaining 60 PCs and 3 discriminant axes.

| Landscape genetic analysis
We conducted a landscape genetic analysis to investigate the possible effects of landscape topography on genetic differentiation in A. koliensis. The most eastern population in the northern Papuan Peninsula in Figure 1 (site 1 in Table 1) was omitted from all landscape analyses due to this population consisting of a single collection site that showed a distinct site-specific genetic signature. To describe topographic complexity across the study area, we obtained elevation data (90 × 90 m horizontal resolution) for PNG from the Shuttle Radar Topography Mission (SRTM; Jarvis, Reuter, Nelson, & Guevara, 2008) and used these elevation data to produce a layer describing landscape topography using the "terrain ruggedness index" (Wilson, O'Connell, Brown, Guinan, & Grehan, 2007). This index assigns greater values to areas that have higher (e.g., hills) or lower (e.g., valleys) elevation than their surroundings, and thus would be expected to positively correlate with dispersal cost (Row et al., 2015).
We assessed four different hypotheses for connectivity among the northern, southern, and combined populations. The first hypothesis was the null hypothesis, and under this hypothesis, F I G U R E 1 Topographic map of Papua New Guinea. Sites on map are Anopheles koliensis sampling sites used in this study (detailed in Table 1) and are color-coded to show the distributions of the three rDNA ITS2-RFLP genotypes identified pairwise measures of genetic variation were fit to a constant. The second hypothesis described the connectivity between a set of populations as a function of the geographic distance between them (Wright, 1943). The third hypothesis used terrain ruggedness to account for spatial variation in dispersal cost across the landscape and least-cost path distances to model connectivity among populations.
The fourth hypothesis also used terrain ruggedness and used commute distances to model connectivity among populations (Chandra, Raghavan, Ruzzo, Smolensky, & Tiwari, 1996 were subsequently used to examine the relative importance of surface topography on connectivity. We used an information-theoretic approach (Burnham & Anderson, 2003) to assess the relative support for each of the four connectivity hypotheses under the northern, southern, and combined populations. Linear mixed-effects models were fit to the pairwise measures of genetic differentiation using the predictor variables for each connectivity hypothesis (as described above in the optimization process). The corrected Akaike information criterion (AICc) statistic was used to assess model performance due to low sample sizes (Burnham & Anderson, 2003). Following standard methodology (Burnham & Anderson, 2003), AICc statistics were calculated for each model, and the corresponding δAICc statistics and Akaike model weights (w i ) were used to assess the relative support for each hypothesis (calculated using the MuMIn R package (Bartoń, 2017)).
We also calculated marginal and conditional R 2 statistics to provide a more intuitive description of model fit (Nakagawa & Schielzeth, 2013). To assess uncertainty in the relative support for each hypothesis (following Dudaniec et al., 2016;Dudaniec et al., 2016), we conducted a bootstrap resampling analysis (10,000 replicates). In each bootstrap replicate, a subset of 75% of the populations in a population set (i.e., north, south, or combined) were randomly selected, the linear mixed-effects models corresponding to each of the connectivity hypotheses (as described above) were refitted to the subset of populations, and the AICc statistics for the refitted models were calculated. After completing all of the iterations, we finally calculated the average rank and percentage of times that each hypothesis was found to have the most support among the bootstrap replicates.

| Ribosomal DNA ITS2 genotype identification
All individuals included in this study were identified as A. koliensis by PCR-RFLP of the ITS2 followed by restriction digest using Msp I (Beebe & Saul, 1995). Three subtle ITS2 restriction profile variants could be identified in the agarose gel that reflected those previously found using 10% acrylamide gels (Benet et al., 2004; see Figure 2).
Of the 345 samples assessed from 33 sites throughout PNG, a com- The ITS2 sequences were generated from plasmid clones drawn from a subset of individuals taken from sites where the genotypes could be regarded as being in sympatry in northern PNG (G1-G2 [site 13] and G1-G3); however, no sites were identified with all three genotypes present. As expected, ITS2 sequences clustered into three genetic groupings (see haplotype network in Figure 2b) that correlate with the three ITS2 genotypes identified by PCR-RFLP analysis.

| Population genetics: COI and rpS9 sequence analysis
The indicate an excess of singletons. This is reflected in the haplotype networks for both loci which are displayed in Figure 3. Given the lack of common haplotypes, this is unlikely to be due to a recent population expansion or selective sweep and probably reflects a very large and stable population supporting very high diversity.
Site-specific spatial graphs showing mtDNA COI haplotype diversity (Hd) are presented in Figure S1.
The haplotype networks in Figure 3 suggest the rpS9 has more shared haplotypes between populations than the COI, and little discernible structure. Both mtDNA and nuclear loci show a high level of haplotype diversity with many unique haplotypes present in single individuals. At the rpS9 locus, a common haplotype was found in all geographic regions sampled, whereas at the COI locus there are shared haplotypes only between the two populations from northern PNG with these four haplotypes being shared between the Sepik region and Madang/Lae region.

| Population genetics: microsatellite analysis
Evidence of null alleles was found in some populations for some loci (Table S2). In STRUCTURE analyses, mean LnP(K) improved rapidly for K = 2 to K = 4 and plateaued at K = 5 and analyses run in CLUMPAK support K = 5 as the most likely K value. The barplot for K = 5 shows additional informative population structure over K = 4 and is pre- New Guinean population appears most closely related to the Sepik population.

| Effects of landscape topography on connectivity
Landscape topography had a large effect on connectivity among the A. koliensis populations at the broadscale (Figure 5a). The model that F I G U R E 2 Three rDNA ITS2 restriction digest variants were found within Anopheles koliensis in PNG (G1, G2, and G3; upper panel). When a subset of these variants were cloned and sequenced, they revealed intraindividual paralogues and three distinct sequence variant lineages with no shared sequences between lineages. The haplotype network (lower panel) shows the genetic relationship of these cloned sequences with genotypes with G1 common through PNG with G2 and G3 showing nonoverlapping geographic restriction in northern PNG best explained genetic variation among all of the populations was the model that accommodated landscape topography and modeled dispersal using least-cost path distances (Table 3). The relative support for this model eclipsed all others (w i = 1), and this model was The connectivity between the northern populations was also affected by landscape topography (Figure 5b). Similar to the broadscale connectivity patterns, the best support connectivity model was the model that accommodated landscape topography between populations and described dispersal using least-cost distances (w i = 0.99; Table 3). This model was able to explain an adequate proportion of the variation in genetic differences among the sampled populations (0.43R 2 m), and was also reasonably robust F I G U R E 3 Haplotype network for the mtDNA COI and nuDNA rpS9 DNA sequences for Anopheles koliensis. The size of the circle reflects number of individuals sharing a particular sequence, and the connections represent single mutational steps. Colors designate regions in PNG were identified as genetically distinct by the microsatellites There was little evidence to suggest that landscape topography had any effect on connectivity among the southern A. koliensis populations. The best supported model was the null model (w i = 0.86; Table 3). However, the support for this model was so, the null model received the best support (observed in 9.94% of the bootstrap replicates). But when the connectivity models were fitted to a subset of populations that did not include the supposed outlier populations, the models were able to adequately describe the genetic differences among the populations, and so, the connectivity models which included landscape topography had much more support (e.g., terrain ruggedness with least-cost distances had the greatest support in 40.01% of the bootstrap replicates).

| D ISCUSS I ON
The mosquito A. koliensis transmits malaria and is a member of a Punctulatus Group of 13 cryptic species of which 11 can still be distinguished by the original ITS2 PCR-RFLP method (Beebe et al., 2015;Beebe & Saul, 1995). In assessing A. koliensis through PNG, we identified three ITS2 RFLP genotypes complementing a study on this mosquito in northwest PNG (Benet et al., 2004). A common genotype (G1 or Madang Wosera in Benet et al. (2004)) exists throughout PNG, with a second less common genotype (G2 or Madang in Benet et al.) sympatric with G1 through the Madang/ Lae region. A third genotype (G3 or Wosera in Benet et al., 2004) was also sympatric with G1 but only through PNG's Sepik region in northwest PNG. Cloning and sequencing of individuals with these ITS2 genotypes revealed intraindividual ITS2 sequence paralogues but no shared sequences between individuals with different genotypes, despite genotypes occurring at sympatric collection sites (G1-G2 and G1-G3). The analysis of microsatellite data suggests that individuals of ITS2 genotype G3 may be reproductively isolated from individuals other with the other ITS2 genotype. Using landscape genetic analysis of microsatellite data, we found that elevation restricts the dispersal of A. koliensis between populations in Papua New Guinea.

| Population structure of A. koliensis in PNG
Population structure is evident between north and south PNG in COI but not rpS9. The high haplotype diversity in both mtDNA COI and F I G U R E 5 Landscape genetic analysis. The top panels show resistance to gene flow maps estimated using landscape topology for (a) all populations and (b) the northern populations. Colors show the spatial distribution of resistance to gene flow, points correspond to sampled populations, and lines represent potential dispersal routes between sampled populations using the resistance data. The bottom panels show the modeled relationship between landscape topology, measured as the terrain roughness index, and resistance to gene flow for (c) all populations and (d) the northern populations. Maps and data for the southern populations are not shown because the landscape genetic models fit exclusively to these populations failed to explain an adequate proportion of their genetic variation rpS9, and an excess of singleton sequences suggest that the species has a long history in the region. The higher effective number of haplotypes and π for the COI and rpS9 sequence data per site (Table 2) together suggest northern PNG may be an older and more stable population than southern PNG. We suggest that the northern New Guinean population may have founded the southern population incurring an apparent genetic bottleneck caused by the Central Range.

The single population in eastern PNG, north of the Owen Stanley
Range (site CP144), appeared genetically distinct at the level of the COI in that it did not share haplotypes with the other populations.
However, haplotypes of individuals from this site appear in different parts of the haplotype network and the lack of haplotype sharing may be an artifact of the small sample size from CP144 (n = 11).
The analysis of the 11 microsatellites provided enhanced detail in regard to the spatial separation of populations through PNG.
The population in southern PNG seen in the COI marker were now clearly identified, and the northern Papuan Peninsula population could be discriminated. Within northern PNG, microsatellites appeared to pull apart the Sepik region populations to the west from the Madang/Lae populations in the east-an area that could be regarded as continuous.

| Effects of landscape topography on genetic connectivity
Overall, the broadscale differences among the A. koliensis populations could be better explained when considering the topography of the region as it is supported by both our model selection analysis and the bootstrap analysis for the combined set of populations. The least-cost path model for mosquito dispersal suggests that movement between populations occurs along distinct routes, rather than a diffusive model for dispersal where individuals take multiple routes between each pair of populations. This phenomenon has also been observed in montane amphibian populations (Kershenbaum et al., 2014;Zancolli, Rodel, Steffan-Dewenter, & Storfer, 2014).
Focusing only on the northern PNG populations, the best supported hypothesis for fine-scale genetic variation among these populations also included topography with a least-cost distance model for dispersal. This result suggests that at the fine scale, individuals in the northern populations may also exploit specific routes for dispersal. The bootstrap analyses corroborate this outcome, but they suggest that geographic distance may also play a minor role in gene flow. This is may be due to the limited number of northern populations that were sampled in the present study, and additional data covering more of the northern populations could potentially reduce this uncertainty.
Unlike the northern PNG populations, the best supported hypothesis for fine-scale genetic variation among the southern populations was the null hypothesis. The relative support for the null hypothesis was far greater than any other hypothesis when fitting models to all of the sampled southern populations. However, the bootstrap analysis revealed that this strong level of support was only present when fitting models to all of the sampled southern populations and that Note: Data show results for maximum-likelihood population effects (MLPE) describing genetic differences among combined, northern, and south populations. For a given model, AICc represents its corrected Akaike information criterion, δAICc is the difference between its AICc statistic and that of the best supported model in the set, and w i is its Akaike weight which denotes the probability that it is the best in the set. The marginal (R 2 m) and conditional (R 2 c) statistics are also reported. To account for uncertainty, the results from the bootstrap analysis are also reported. These show the mean rank of the models and the percentage of times that each model was the best in its set. the relative support for the null hypothesis was much lower when specific populations were omitted. Therefore, it seems likely that geographic distance and landscape topography have some effect on the spatial patterns of gene flow among southern populations.
However, to understand the effect of landscape topography on gene flow among southern PNG A. koliensis populations, more extensive geographic sampling from this region is needed. These results further highlight the importance of quantifying uncertainty among competing hypotheses in landscape genetics (Dudaniec et al., 2016).

| Reproductive isolation between sympatric rDNA variants and divergent biting behaviors in Anopheles mosquitoes
Microsatellite data suggest that individuals of ITS2 genotype G3 were reproductively isolated from individuals of the more common G1 and G2 ITS2 genotypes. Individuals of ITS2 genotype G3 appear to be spatially restricted to the western Sepik region of northern PNG, while individuals of ITS2 genotype G2 seem to be restricted to the Madang/ Lae region, which supports the observed gene flow restriction between the Sepik and Madang/Lae populations. The ITS2, being part of the rDNA gene family, is tandemly arranged multicopy gene family in metazoans, and the evolutionary machinery maintaining sequence fidelity between copies in this array is not well described and does not follow traditional Mendelian rules of inheritance (Bower, Cooper, & Beebe, 2009;Eickbush & Eickbush, 2007;Nei & Rooney, 2005). In Anopheles mosquitoes, the rDNA is usually positioned on the sex chromosomes adjacent to the centromeres (Kumar & Rai, 1990) and the relatively rapidly evolving ITS2 spacer has been a useful marker for detecting early genetic discontinuity between populations (Alquezar et al., 2010;Beebe, 2018;Coleman, 2009;Muller, Philippi, Dandekar, Schultz, & Wolf, 2007). The positional effect of the rDNA being adjacent to the centromere would be to reduce recombination (Nachman & Churchill, 1996). In Anopheles, elevated levels of genetic divergence have been observed in regions proximal to the X chromosome in the face of gene flow across other parts of the genome (Reidenbach et al., 2012;Weetman, Wilding, Steen, Pinto, & Donnelly, 2012). If alleles for traits pertaining to time of night-biting behavior occur near or within this low recombination landscape, the appearance of intraspecific behavioral differences may associate with rDNA divergence.
Interestingly, the Benet et al. (2004) study observed variation in nightbiting behavior between the A. koliensis genotypes in northern New Guinea. Their study found evidence that the Madang (M) variant (G2 in this study) starts to blood feed later in the night, where MW (G1) and W (G3) were actively seeking a host as early as 6 p.m. It would be reasonable to hypothesize that genes for circadian or other rhythmic actives may be positioned near the centromere on the X chromosome.
Individuals of ITS2 genotype G3 are spatially restricted to the Sepik region of northern PNG and form a clearly separate group in the microsatellite clustering analyses, despite being sampled from the same sites as individuals of ITS2 RFLP genotype G1. The distinct microsatellite profile of these genotypes despite their overlapping ranges provides evidence that A. koliensis may well be more than one species in PNG. Individuals carrying the G1 and G2 genotypes show no evidence of genetic structure at other loci despite appearing to have different host feeding initiation times (Benet et al., 2004).
We found that the geographic structure observed through PNG is best explained by landscape topography, with slope (or elevation) presenting as a significant factor. This makes sense given that the species exhibits a predominantly a lowland distribution . Although the individuals of ITS2 G1 and G2 may have differences in time of feeding, they could not be separated by fastevolving microsatellites. Their distinct ITS2 genotypes may however represent early stages of genetic discontinuity which may, in time, lead to reproductive isolation across other parts of the genome.

ACK N OWLED G M ENT
The authors would like to thank Darren Waterson and Steven U19AI089686-03).

AUTH O R CO NTR I B UTI O N S
LA completed laboratory work and much of the genetic analyses and manuscript writing; JOH performed the landscape genetic analyses and wrote these sections in the manuscript; CR assisted with the genetic analyses and contributed to writing of the manuscript; WX performed a subset of the molecular genotyping, nDNA, and some mtDNA sequencing, while SLF initiated a substantial amount of the genotyping and mtDNA COI sequencing; RDC was responsible for the field collections; and NWB designed the project, contributed to field collections, and assisted in developing the manuscript.