Mitogenome analysis included both newly acquired and sequenced samples of A. coluzzii from São Tomé and Príncipe islands, and publicly available data generated previously for this species (Figure 1; Table S1). The former sample set included immature stages of A. coluzzii collected from a total of 73 breeding sites in STP (38 localities on São Tomé and 35 on Príncipe) (Table S1). The latter sample set was comprised of A. coluzzii from the UC-Davis Vector Genetics Laboratory (VGL) archive (21) or from the Ag1000G database, phase 2 (22). These samples originated from STP (São Tomé, N=14; Príncipe, N=15) and 10 countries in West and Central Africa: Angola (N=15), Benin (N=10), Burkina Faso (N=15), Cameroon (N=9), Côte d’Ivoire (N=15), Equatorial Guinea (here after Bioko Island, N=5), Gabon (N=5), Ghana (N=15), Guinea (N=4), Mali (N=14).

Individual mosquito genomic DNA from specimens taken from the VGL sample archive were extracted using a Qiagen Biosprint following our established protocol (23). Species diagnostics was performed using species-specific SNPs included in the DIS assay (24). The methods used for DNA extraction, species identification and genome sequencing of individuals from Ag1000G were as described by The Anopheles gambiae 1000 Genomes Consortium (22).

Individual mosquito DNA from the VGL samples was measured using a dsDNA high sensitivity assay kit on a Qubit instrument (Thermo Fisher Scientific, Waltham, MA, USA). Individual library preparation was done with 10 ng of genomic DNA as input and using KAPA HyperPlus Kit (Roche Sequencing Solutions, Indianapolis, Indiana, USA) following our protocol (25). AMPure SPRI beads (Beckman Coulter Life Sciences, Indianapolis, Indiana, USA) were used for library size selection and clean-up. Sequencing was performed on an Illumina HiSeq 4000 instrument at the UC Davis DNA Technologies Core facility for a pooled library of equal amounts of individual indexed libraries.

After demultiplex and adapter removal, raw-sequencing reads were used to assemble the mitochondria contig using NOVOPlasty version 2.6.7 (26) and 33 as K-mer’s value and default settings. Samples from the Ag1000G were downloaded as binary alignment map (BAM) files from the European Nucleotide Archive (ENA) under accession number PRJEB36277. BAM files were converted to FastQ format using BEDTools (27) which were then used for mitochondria assembly as described above.

De novo assembled mitochondrial genomes were imported into Geneious (2021.1.1) and aligned with a reference mitogenome sequence of Anopheles coluzzii (JADFTP010000004.1) acquired from GenBank. The alignment was visually inspected to confirm that no missing data, singletons or polymorphism/ambiguities were present in the sequences. Genetic diversity indices such as nucleotide diversity (π) and haplotype diversity (Hd) were assessed for each population using DnaSP 6.12.03 (28).

The annotated mtGenome of A. gambiae (NCBI Ref. Seq. NC_002084.1) was used to partition the aligned mitogenomes. We conducted a partition test of heterogeneity to determine if the complete mtGenome was appropriate to use for phylogenetic analysis, as implemented in PartitionFinder 2.7.1 (29). The models of evolution that best fit each partition were determined by PartitionFinder 2.7.1 (Table S1). The A. coluzzii sequence alignment was aligned with mitogenome sequences of A. melas (KT382823.1), A. arabiensis (KT382816.1), A. gambiae (NC_002084.1) (Table S1) to use as additional outgroups for phylogenetic analyses using MAFFT v7.450 (30, 31).

A maximum likelihood (ML) analysis was conducted using IQ-TREE v2.1.3 (32) to estimate a ML tree with 1,000 ultrafast bootstrap replicates and resampled partitions and then sites within resampled partitions (32–36). Likelihood settings followed the general time reversible model (GTR) with a gamma distribution and estimates of the proportion of invariable sites (+I+G). Confidence in the resulting topologies was assessed using rapid bootstrapping and a search for the best-scoring tree with 1,000 replicates. We performed Bayesian inference (BI) analyses using parameters selected by PartitionFinder 2.7.1 (37) and conducted in MrBayes 3.2.7a (38). An MCMC algorithm ran for 100,000,000 generations, sampling one tree every 10,000 generations. Likelihood scores were used to determine the value of burn-in and stationary distributions. When split frequency in the Bayesian analyses reached <0.01, a 50% majority-rule consensus tree was obtained from the remaining trees. Posterior probabilities for clades were compared for congruence between analyses, bootstrap values ML > 70 and BI > 0.90 are presented on the phylograms.

Haplotypes were identified and corrected using Poppr 2.0 and DnaSP 6.12.03 (39). An initial TCS haplotype network was constructed using POPArt 1.7 (40). Mainland Africa samples that differed by >25 mutations from individuals from São Tomé and Príncipe were removed and a second haplotype network constructed. A third haplotype network including only individuals from São Tomé and Príncipe Islands was constructed to identify the number of haplotypes and their geographic distribution.

Patterns of migration between mainland Africa and STP were investigated using a Bayesian approach compiled in Migrate-n v. 4.4.3 (41). Parameters were estimated under a full migration model that allowed gene flow to occur among STP and their most recent common ancestral populations. Migrate-n analyses were conducted using default parameters with two independent runs sampling every 200 steps for 200,000 recorded steps and a burn-in of 20,000. The migration models were compared based on their marginal likelihood and probability using thermodynamic integration with Bezier approximation as implemented in Migrate-n (42). Visualization of the migration pattern inferred by Migrate-n analysis was constructed using the Migest Package in R (43).

In order to create a mitogenome-based timeline for the introduction of A. coluzzii onto STP, as reflected in their mitogenome sequences, additional phylogenetic trees were estimated. For this purpose, a multi-locus coalescent model was applied to the data using the software BEAST2 [BEAST 2.5; (44)]. Intraspecific divergence times were concurrently estimated using a substitution rate of 1.2 x 10^-8 mutations per site per year and assuming a generation time of three weeks (45). BEAST2 analyses were performed in triplicate at the University of California – Davis Genome Center High Performance Computing cluster with Markov chains run for 100 million generations or until convergence, with 10 million generations of each run discarded as burn in and chains sampled every 10,000 generations. Based on the results from PartitionFinder 2.7.1 the HKY substitution model was used in combination with gamma site-specific rate variation and a proportion of invariant sites parameters with a relaxed log normal molecular clock. A strict clock was selected over a relaxed clock comparing the likelihood of the results of 10,000,000 generation test runs. Tracer v1.7.1 (46) was used to assess convergence. LogCombiner (47) was used to resample 10000 subtrees from BEAST2 analyses and then TreeAnnotator (44) was used to generate Maximum Clade Credibility (MCC) trees (44).

A total of 456 mitochondrial genomes were used to investigate the origin of A. coluzzii in São Tomé and Príncipe islands. New complete A. coluzzii mitochondrial genome sequences were generated for 158 individuals from São Tomé and 158 individuals from Príncipe. After multi-alignment, variable sites were identified in each sequence. São Tomé and Príncipe presented the lowest number of polymorphisms and nucleotide diversity estimates, despite being the largest sampling dataset (Table 1). Haplotype diversity ranged from 0.134 to 1, in which the Príncipe population had the lowest value (Table 1).

This combined phylogenetic analysis (Figure S1A) represents 414 terminals from the genus Anopheles, representing four species (A. melas, A. arabiensis, A. gambiae and A. coluzzii) and 12 western African (WA) sites (Angola, Benin, Bioko, Burkina Faso, Cameroon, Cotê d’Ivoire, Gabon, Ghana, Guinea, Mali, São Tomé, and Príncipe), and includes at least one representative of every unique haplotype (Figure S1B).

All relationships within the mitogenome phylogeny were resolved with significant support using ML and Bayesian analyses (ML > 70 and BI > 0.90; Figure 2). Anopheles arabiensis is recovered as the earliest branching lineage and is a sister clade to A. gambiae + A. coluzzii (Figure 2). The A. coluzzii reference sequence (JADFTP010000004.1) is recovered as the earliest branching lineage within A. coluzzii and is sister to five WA clades. Individuals from São Tomé and Príncipe are recovered as a monophyletic clade and are sister to individuals from only the northwestern region of West Africa (NW: Benin, Burkina Faso, Cotê d’Ivoire, Ghana, Guinea, and Mali) with strong support (Figures 2 and S1).

Individuals from Central Africa (Angola, Bioko, Cameroon, Gabon) form an unresolved cluster from haplotypes within the network, with northern West African (NW) haplotypes appearing as the closest relatives of haplotypes present in STP (Figure S1). No shared haplotype between any west African country and STP were identified (Figures 3A and S1). São Tomé and Príncipe haplotypes are separated from the closest mainland haplotypes by a least 11 mutations (Figure S1). We identified 5 unique haplotypes present in São Tomé (ST-1, 2, 3, 5 and 6), two unique haplotypes in Príncipe (PR-1and 2) and three shared haplotypes between the two islands (Figure 3A). We found no geographic substructure of haplotypes within either island. Migrate-n analysis indicates similar migration estimates of A. coluzzii among the west African populations in the countries of Benin, Burkina Faso, Cotê d’Ivoire, Ghana, Guinea, and Mali, and from these populations to those in São Tomé (Figure 4). This analysis also indicates migration from São Tomé to Príncipe, with an insignificant chance of migration in the opposite direction (from Príncipe to São Tomé). Migration from both islands to any of continental populations also appears highly unlikely (Figure 4).

We identified a well-supported phylogeny generated by BEAST2 that is concordant with the ML and BI mitogenomic phylogenetic trees. The final BEAST2 trees support a recent divergence of the São Tomé and Príncipe populations from ancestral West African populations (Figures 3B and 5). Divergence time estimates reveal that São Tomé and Príncipe populations began differentiating from West African populations around 500 years before present (Figure 5). Our BEAST2 trees suggest that Príncipe populations represent a contemporary divergence from São Tomé. Furthermore, the BEAST2 trees support at least two introductions of A. coluzzii into São Tomé and at least four introductions from there to Príncipe (Figure 5).

Our phylogenetic analyses and haplotype networks using whole mitochondrial sequences support a West African (Mali, Burkina Faso, Côte d’Ivoire, Guinea, Ghana, Benin) ancestry of A. coluzzii on the islands. Similarly, a population genetic study has presented the smallest F_ST values between STP populations and West African populations from Ghana, Senegal and Côte d’Ivoire using a region of the mitochondrial NADH dehydrogenase subunit 5 (ND5) gene (19). Discordant results were found using nuclear sequences that suggested the population of A. coluzzii in Angola and Central Africa (Bioko, Cameroon, and Gabon) as the genetically closest populations to STP (19, 21). The biology of mitochondrial and nuclear DNA (nDNA) includes fundamental differences such as ploidy, degree of recombination, number of introns, effective population size and mutation rate (48, 49). For A. coluzzii, it is a comparison between the 15,000 base pair mitochondrial genome and 278 million base pair nuclear genome (50). However, mtDNA has been extensively used to evaluate ancestry and demographic changes in populations (51, 52).

Our analyses suggest that A. coluzzii was introduced onto the island of São Tomé from West Africa on at most two occasions. Historical information suggests this introduction occurred during the first major wave of human migration associated with the Portuguese colonization and subsequent slave trade supporting sugar production around 1500-1600_AD (53). The introduction of malaria into São Tomé reported in 1493 (54) supports the presence of the vector at around that time. Our results also suggest that A. coluzzii was only introduced to Príncipe Island from São Tomé four or possibly five times. Of critical importance, we also found no evidence of contemporary genetic exchange between STP and mainland Africa, or between São Tomé and Príncipe which supports the suitability of STP for early field trials of genetically engineered mosquitoes for malaria elimination from the islands.

Within the island of São Tomé, we identified five unique haplotypes (ST-1, 2, 3, 5 and 6) which appear to be diverged from a single ancestral lineage, designated lineage ST-1, (Figure 3A). On Príncipe we identified 2 unique haplotypes that likely originated on São Tomé from one to two introductions (Figure 3A). We also observed three shared haplotypes between the islands that, combined with our bidirectional migration analysis, likely diverged on São Tomé and were subsequently introduced to Príncipe (Figures 3A;  4). The presence of multiple haplotypes on each island suggests that at some point populations were isolated. However, the lack of geographic substructure in haplotype distributions suggests that contemporary gene flow among subpopulations is occurring within each island. The small number of haplotypes (Figure 3A) and absence haplotype divergence among subpopulations within each island agrees with a very recent colonization of STP, resulting in the observed population bottleneck (founder effect) of A. coluzzii in STP (21).

As geographic distance from a source population increases, overall genetic diversity decreases, resulting in fewer observed mtDNA haplotypes. This can be the result of colonization events on islands causing population bottlenecks (founder effect) (21, 55). The founder effect is amplified as the number of islands within a chain increases, because colonization typically occurs in a stepwise pattern moving away from the island first colonized, resulting in multiple population bottlenecks. The presence of a greater number of haplotypes (8 vs 5) on the larger of the two islands (São Tomé vs Príncipe) fits within island biogeography theory and lends further support to São Tomé as the initial island of A. coluzzii introduction (56). The results reported here provide additional evidence that populations of the malaria vector A. coluzzii are genetically isolated from mainland conspecific populations. In addition, our results suggest that introductions of this species into STP are rare events and that there is little or no contemporary dispersal of this species into the islands from mainland populations.