Species assignment methods commonly work by comparing a query sequence to a reference database (Ratnasingham and Hebert, 2007; Benson et al., 2018). The completeness and quality of the reference database heavily influence the accuracy of the assignment method. The reference database we constructed consists of well-curated samples sequenced with the panel, in silico extracted reference genomes, and in silico extracted whole genome short read data. As we expect that the reference database will expand to include additional species and populations over time, we assign a version number to the database. The reference database described here, which we call NNv1, contains 186 samples, representing 62 species spread over 4 subgenera. The dataset from Makunin et al., 2022 forms the backbone of the reference panel. In addition, six species in the An. gambiae complex have been included from publicly available whole genome data (The Anopheles gambiae 1000 Genomes Consortium, 2021; Fontaine et al., 2015) in order to increase the resolution in a group of hard to distinguish species. To maintain the advantages of multi-locus assignment, we required samples to have at least 10 targets amplified to be included in the reference database. Ideally, the database would contain several specimens per species to represent within species variation. This is particularly important for species with a wide geographical range.

The amplicon panel is designed to improve accuracy of species assignments over morphological or single-marker methods. However, the species labels supplied by our sample partners were in most cases obtained by the latter methods, hence it was necessary to reconfirm them. All label information for the reference database NNv1 is listed in Supplementary file 1, and the assignment of species labels is discussed in more detail in Appendix 1, but we present an outline of the principles we used here. For most samples sequenced by the amplicon panel, two molecular barcodes (COI and ITS2) were also sequenced and compared to the sequences available in BOLD (Ratnasingham and Hebert, 2007) and NCBI (Benson et al., 2018). Here, we used the barcode information as well as the pairwise k-mer distances between samples in the reference database NNv1 to generate a consensus species label for each sample. A few samples show inconsistencies between their partner labels, molecular barcodes and amplicon assignments; these samples were also removed or flagged as overly diverged in Makunin et al., 2022. For these samples, the distances to samples of the same species label is much larger than the distance to some other samples in the database, suggesting they were mislabelled. In some cases, there is additional evidence from molecular barcodes and it was possible to relabel them to the species they match. Other samples are clearly different from samples with the same species label, but do not clearly belong to any other species in the database. These are labelled as the species-group they belong to appended by ‘_sp1’, ‘_sp2’, etc. Hopefully, by extending the reference database in the future, we can get a better understanding of which species such samples represent. For some closely related species, the different species labels were supported neither by the pairwise-distances nor by the barcodes. Samples from these species are assigned to a group of closely related species and their consensus labels end in ‘_c’ to notify that they are part of a complex of species that cannot be distinguished by the nearest neighbour method. None of the in silico extracted samples showed inconsistencies with pair-wise k-mer distances and hence all retained their published labels.

The species-groups are defined based on the pairwise k-mer distances between samples in the reference database NNv1 and make use of the consensus species labels discussed above. The k-mer distance between two samples, s₁ and s₂, is defined as $d_{s,k}\left(s_1,s_2\right)=\frac{1}{\left|T_1\cap T_2\right|}{\sum }_{t\in T_1\cap T_2}^{}{\sum }_{q_1\in Q_{1,t}}^{}{\sum }_{q_2\in Q_{2,t}}^{}\frac{d_k\left(q_1,q_2\right)}{\left|Q_{1,t}\right|\left|Q_{2,t}\right|}$, where T_i is the set of targets amplified in sample i and Q_i,t is the set of unique haplotypes of sample i at target t and by |S| we mean the number of elements in set S. The k-mer distance d_k between two sequences is defined in the Methods section on k-mers. So in words, for a given target the k-mer distance between two samples is the mean k-mer distance between all pairs of haplotypes from the two samples; for example if sample s₁ is homozygous and s₂ heterozygous at target t, then the k-mer distance at target t between these samples is the mean of the k-mer distance of the haplotype from s₁ compared to the first haplotype from s₂ and the k-mer distance of the haplotype from s₁ compared to the second haplotype from s₂. If a sample has more than two alleles at a single target, we take into account all haplotypes according to the above definition. The k-mer distance between two samples is defined as the average of the k-mer distance between these samples at all the targets that were amplified in both samples.

Figure 1 shows the pairwise k-mer distances between all samples in the reference database. The samples are ordered roughly by phylogeny (as in the tree in Figure 2) and this results in a visible structure in the distance plot. One can observe dark triangles below the diagonal, reflecting that samples of the same or closely related species have a smaller 8-mer distance to each other than to other samples. This effect repeats itself on different scales, mirroring the multi-level structure in the phylogenetic tree.

Figure 1 with 1 supplement see all

Download asset Open asset

Figure 2

Download asset Open asset

We identify clusters of samples in the reference database NNv1 and we will refer to those as species-groups. It is important to emphasise that the way we use the term ‘species-groups’ does not refer to a taxonomic classification, but to the clusters of samples defined by thresholds on the 8-mer distances between samples. Ideally, a threshold would partition the samples in the reference database NNv1 into species-groups, such that the 8-mer distance between each pair of samples within the group is smaller than the threshold and the 8-mer distance from any sample within the group to any sample outside the group is larger than the threshold. The motivating assumption is that the samples in the same species-group share more recent ancestry with each other than with samples that are members of different species-groups. While several thresholds satisfy the partitioning condition, they are most useful for the purpose of species classification when they generate species-groups that (approximately) correspond to classification at an established taxonomic level. Not surprisingly, there are no thresholds that perfectly satisfy the partitioning condition and are completely concordant with current taxonomy, but several thresholds give groups that are close to taxonomic entities and given the paucity of molecular data for some of these species, it may be that the taxonomic entities are less phylogenetically accurate than the k-mer based groupings. In exploring different thresholds (Figure 1—figure supplement 1), we selected three levels that best matched described taxonomic entities. The species-groups at each of the threshold levels we have selected are listed in Supplementary file 1.

For very low thresholds, each sample will form its own species-group, which is not informative. For a threshold of 0.1 on the sample 8-mer distance, most species-groups satisfy or nearly satisfy the partitioning condition, and each species-group contains a single species or multiple species from a known species complex or group, for example An. gambiae and An. coluzzii are grouped together at the 0.1 threshold. We refer to the species-groups at the 0.1 level as the fine level species-groups and they are most useful for species assignment, because they provide the highest resolution. A threshold of 0.3 merges samples representing species from some well-known complexes, like the entire An. gambiae complex, into a single species-group. Similarly, many of the species for which we currently only have a single representative in the reference database NNv1, are merged into larger species-groups. We refer to this level of species-groups as the intermediate level, which provides insight into the degree of similarity between different fine level species-groups. Additionally, the intermediate level species-groups can be informative when we sequence a sample whose species is not represented in the database because the species assignment at the intermediate level places the sample within its most closely related species-group in the database. Thresholds higher than 0.3 tend to violate the partitioning condition to a greater extent, but it is desirable to include a coarse level classification to get an approximate taxonomic position of unrepresented species and of more diverged species. We selected 0.51 as a threshold for the coarse level classification because this gives reasonably clear species-groups that roughly correspond to taxonomic series. However, it is not perfect, as it groups together the Myzomyia and Neocellia series in the Cellia subgenus and it splits the Neomyzomyia series into three distinct species-groups. A similar effect was observed in Makunin et al., 2022, where the Neomyzomyia series did not form a monophyletic clade.

At the fine level, most species-groups contain all samples from a single species, but there are some exceptions. For some species, for example An. nili and An. hyrcanus, the samples are split into multiple fine-level species-groups, because they appear much more distinct from each other than you would expect in a single species (see Appendix 1 and Figure 1—figure supplement 1 for more detailed discussion). In our assignment, these will be treated as distinct species, highlighting to entomologists and taxonomists that further work to refine species in these groups is needed. Conversely, there are also species-groups that contain samples from multiple different species. These do not represent a single species, but instead they represent a complex of closely related species. Species inside species complexes often share a lot of genetic variation and the k-mer distance based method that we discuss here does not have sufficient resolution to reliably distinguish between them. The species-groups at the fine level that contain more than a single species are the An. marshallii group (contains An. hancocki, An. brohieri, and An. demeilloni), An. gambiae/coluzzii, the An. sundaicus complex (contains An. sundaicus and An. epiroticus) and the An. coustani group (contains An. coustani, An. ziemanni, An. tenebrosus, and An. paludis).

The first step of the hierarchical assignment method performs nearest neighbour assignments to samples in the reference database at the three different levels of species-groups introduced above. The assignment is initially done independently at each target, for each sample, computing assignment proportions for the species-groups at the chosen level, normalised such that they sum up to one over all species-groups. The resulting per-target assignment proportions are then averaged over all targets, resulting in the overall sample assignment proportions at the chosen level. If the sample assignment proportion is at least 0.8 for one species-group, the sample is classified as a member of that group. If the classification threshold of 0.8 is not met at the chosen level, the sample remains unassigned at that level.

For example, to assign a sample s at the coarse level, we translate its target sequence corresponding to target 1 to an 8-mer count table, denoted as q_s,1. Next, we compute its 8-mer distance, $d_8\left(q_{s,1},q_{r,1}\right)$ , to every target sequence in the reference database corresponding to target 1, that is to every $q_{r,1}$ . The nearest neighbours of the test sequence q_s,1 are those sequences in the database that minimise the 8-mer distance between themselves and q_s,1. In other words, the nearest neighbour sequences of q_s,1 are the target 1 sequences in the database with the highest percentage of matching 8-mers. The nearest neighbour sequence of q_s,1 can be a sequence that occurs in a single sample in the reference database, or the nearest neighbour sequences can be the same sequence occurring in multiple samples in the database, or the nearest neighbour sequences can be distinct sequences that have the same distance to q_s,1.

Now we bring in the species-groups. For each species-group, we record the frequency of nearest neighbour sequences. This can be thought of as an ‘allele-frequency’ when we classify each target sequence as either a ‘nearest neighbour allele’ or not a ‘nearest neighbour allele’. But just like allele frequency, it takes into account the zygosity of the samples. The nearest neighbour frequencies are normalised, such that they are equal to one when summed over all species-groups. These quantities are the per-target assignment proportions.

This procedure is repeated for every amplified target in the test sample. Finally the per-target assignment proportions are averaged over all successfully amplified targets to give the overall assignment proportions for sample s at the coarse level. If there is a species-group with an assignment proportion of at least 0.8, the sample is classified as a member of this group, otherwise it remains unassigned at this level. Assignments to the intermediate and fine levels are made in the same fashion, starting from the same nearest neighbour assignments, but based on the relevant species-group memberships.

The per-target assignment proportions are based on the frequency of nearest neighbour sequences in the species-groups and not simply on the count of nearest neighbour sequences. The use of frequencies corrects for the different sizes of the species-groups. Suppose a nearest neighbour sequence occurs in 2 out of 10 samples of species-group A and in 1 out of 2 samples in species-group B (and assume all samples are homozygous). By using counts, we would attribute 2/ (2+1)=0.67 and 1/ (2+1)=0.33 assignment proportion to species-group A and B, respectively. But if we did this, species-group A would only have a higher assignment proportion because it contains more samples. What we are really interested in, is how similar the target sequence of the test sample is to the target sequences in the species-groups. So by using the nearest neighbour frequencies, the assignment proportions are (2/10)/(2/10+1/2)=0.29 and (1/2)/(2/10+1/2)=0.71 for species-group A and B, respectively.

The nearest neighbour frequencies observed for a given target of a given sample at a given assignment level are normalised to obtain the per-target assignment proportions. This normalisation ensures that every target is weighted equally. Without the normalisation, the weight of a target would be determined by a combination of the overall frequency of the nearest neighbour sequences and their distribution amongst the species-groups, whilst we are more interested in the distribution than the total frequency.

Targets that did not amplify in the sample are simply ignored. There are different reasons why a certain target does not get amplified in a sample. It might be that the primer binding sites for the target are too diverged or altogether absent in the test sample’s genome. Or the target might not be amplified due to technical reasons like poor DNA quality. In the former case, we are implicitly using the information contained in the missingness, because we restrict our attention to the targets that did amplify in the sample, which should be the same targets that amplified in the samples of the same species contained in the reference database. In the latter case, as long as the missingness is randomly affecting the targets, ignoring missing targets does not bias the assignment proportions.

So far, we have assumed that the test sample was homozygous at each target. To generalise to the heterozygous case, we compute the per-target assignment proportions separately for both target sequences and average them for the final per-target assignment proportions. It does occasionally happen that a sample has more than two different target sequences. This can be due to errors in the PCR amplification or sequencing, contamination by other samples, or a certain target region might be duplicated in the genome of some species. In the NNv1 reference database on average 1.2% of targets per sample have more than two different sequences; in the query datasets discussed later in this article this percentage ranged from 0.3% to 2.3%. We deemed this small enough to simply extend the per-target assignment proportion computation to include targets with more than two different sequences, by taking the average assignment proportions over different sequences as for the heterozygous case.

Figure 2 shows two examples of the nearest neighbour assignment of a test sample. The test samples are Amou-3–2 and Agam-35, an An. moucheti and An. gambiae individual, respectively. For Amou-3–2 we see that for most targets, the nearest neighbour sequence is found in all four An. moucheti samples in the reference database NNv1. For some targets, the nearest neighbour sequence is only carried by a subset of the An. moucheti samples in the database, whilst for other targets the nearest neighbour sequence is also carried by individuals of other species. For Agam-35 the heatmap shows that, for many targets, the nearest neighbour sequences are not only found in An_gambiae_coluzzii samples, but also in samples from other species in the An. gambiae complex. There are only two matches to samples outside the Pyretophorus series, not shown here. This results in a high-confidence assignment to the Pyretophorus series at the coarse level as well as a high-confidence assignment to the An. gambiae complex at the intermediate level. At the fine level, the largest assignment proportion is to the An_gambiae_coluzzii species-group, but it does not meet the 0.8 classification threshold because of the relatively high assignment proportions to other species-groups within the An. gambiae complex. So at the fine level, this sample cannot be classified with sufficient confidence to a single species. Later, we will present a method to resolve the species identities of samples within the An. gambiae complex.

The species-group assignment has been tested on the reference database itself, by dropping out one sample at a time. The majority of samples could be assigned to the correct species-group at the fine level when using a threshold of 0.8 assignment proportion, see Figure 3, Figure 3—figure supplement 1 and Supplementary files 2–4. To provide context, we have included a phylogenetic tree constructed from pairwise 8-mer distances using FastME (Lefort et al., 2015) and displayed using TreeViewer (Bianchini, 2021). If we ignore the samples that form a species-group on their own, because we do not yet have sufficient representation for those species, 61.8% of samples are assigned correctly at the fine level, and 98.8% and 100% at the intermediate and coarse level respectively. In most cases, the fine level species-groups consist of a single species, although in some cases they comprise multiple species. The jump in assignment success from the fine to the intermediate level is mostly caused by the An. gambiae complex, which is well-represented in the reference database. Most samples within the complex can be assigned to the correct fine level species-group to some extent, but they only meet the assignment threshold at the intermediate level species-group assignment, where all samples in the An. gambiae complex are grouped together. This effect is seen in a few other groups as well, and motivates the VAE part of our assignment procedure (described below).

Figure 3 with 1 supplement see all

Download asset Open asset

For the same group of samples, the average correct assignment proportion per sample at the coarse level is 99.4%. This shows that for samples of species that are well-represented in the reference database, there is a near perfect assignment at the coarse level. And in fact, the average correct assignment proportion of 95.8% at the intermediate level shows that these assignments are also generally with high confidence. At the fine level the confidence starts to break down for some samples, in particular the An. gambiae complex, which is represented by many samples in the reference database. Then the average correct assignment per sample is 81.3%. This shows that an additional classification method for species complexes is desirable.

In NNv1, 21 species-groups at the fine level consist of a single sample. These cannot be assigned to the correct species-group in the drop-out assignment experiment. 12 of these samples are members of a larger species-group at the intermediate level and 11 of them can be assigned at the intermediate level, the other sample is not assigned at the intermediate, nor at the coarse level. The remaining nine samples only become a member of a larger species group at the coarse level. Six of them can be classified at this level, the other three remain unclassified.

The four unclassified samples are An. christyi, An. atroparvus, An. oryzalimnetes, and An. cruzii. All of these are quite diverged from everything else in the reference database and as such do not exhibit a strong matching to any of the coarse level species-groups. In particular, the Kerteszia and Nyssorhynchus subgenus are underrepresented, both in number of species and number of samples, but also the basal species in for instance the Pyretophorus series are not well represented and are too diverged from the other species in this series to exhibit strong similarity to the other samples from this series. In the case of An. cruzii in the Kerteszia subgenus it is actually impossible to assign it to its coarse level species group with the current threshold of 0.8, because of its 26 amplified targets, only 15 are also amplified in the single other sample in the reference database from the Kerteszia subgenus, which results in a theoretical maximum assignment to the correct species-group of 0.58. Note that this is the only species-group at the coarse level for which the theoretical maximum correct assignment for a sample in the reference database is smaller than the threshold. But it again underlines the necessity to extend the reference database, both in number of species and number of samples per species.

When a sample is of a species not represented in the reference database, three things can happen. If its species is much more closely related to a single species in the database than to all the others, it will likely be assigned to the species it is related to. Alternatively, if the database contains multiple closely related species, it will be assigned at a higher level to the group that contains all these closely related species. If the sample is highly diverged from all species represented in the database, its nearest neighbours will essentially be chosen at random, and the assignment threshold will not be met. This emphasises the importance of extending the reference database, both by increasing the number of species represented and by increasing the number of samples per species, with a particular focus on capturing the within-species diversity (e.g. representing the geographic species range, representing all karyotypes when polymorphic chromosomal inversions are present).

The species-groups with a single representative can be used to explore the three possible scenarios for when the test sample belongs to a species not present in the reference database: assignment to a closely related species, assignment to a group of species at the intermediate or coarse level, or no assignment.

The first scenario is represented by An.bellator. The only other representative of the Kerteszia subgenus in our reference dataset is a sample representing An. cruzii. At the fine level species-groups, the An. bellator sample is assigned for 0.875 to An. cruzii and hence meets the threshold to be classified at the fine level. So it can happen that a sample is classified to the wrong species, if this species is the only reasonably close species in the database. This scenario is more likely for diverged groups of species with little representation in the database, for example the Nyssorhynchus and Kerteszia subgenera.

An example of the second scenario is the sample Amar-3–1, in the Myzomyia series. At the fine level, it has an assignment proportion of 0.474 to the An_marshallii_complex and 0.378 to the An. theileri species-group. All other species-groups have an assignment proportion of less than 0.06. So this sample is more related to these two species-groups than to anything else in the reference database, but it does not belong to either of them. At the intermediate level, it has an assignment proportion of 0.850 to the An_marshallii_group species-group, which is the An_marshallii_complex and An. theileri species-group combined. Now it meets the threshold and it will be classified as a member of the An_marshallii_group species-group.

An. christyi is an example of a sufficiently diverged sample that does not reach high assignment proportions for a single species-group at any level. At the coarse level, it is assigned to the Pyretophorus_series with 0.506 and the Myzomyia_Neocellia_series with 0.392 and other species-groups have much lower assignment proportions. Hence, it is not possible to classify this sample as a member of any species-group. Adding more samples from this and other underrepresented species to the database, would increase its power to classify samples from those species.

The nearest neighbour approach is not able to confidently distinguish between closely related species that share a lot of genetic variation, but it does identify samples of those species as members of a species complex at the intermediate assignment level. To resolve the species identity inside these species complexes, we use a variational autoencoder approach specifically trained for the complex under consideration. We demonstrate this method for the An. gambiae complex, both because many of our samples fall within this complex and it is medically relevant to be able to distinguish them, and also more practically, because we have access to a large dataset of species-labelled samples. The classifier we present here contains 7 out of 8 formally named species in the complex (Coetzee et al., 2013), as well as two putative cryptic species (Tennessen et al., 2021; Barrón et al., 2019). Compared to NNv1, three additional species are represented in this classifier to present as much of the diversity of the An. gambiae complex as possible. We expect that this method can be applied to other species complexes for which large species-labelled datasets are available.

The An. gambiae complex classifier has been constructed using a training set (GCref v1) and a validation set (GCval v1) of species-labelled samples. We included as many described species in the An. gambiae complex as we could find or generate sequence data for. For species with wide geographic ranges and a large amount of genomic data available, we also aimed to represent the diverse geography where possible.

Both GCref v1 and GCval v1 consist of amplicon sequences and in silico extracted published samples (Fontaine et al., 2015; Neafsey et al., 2015; Tennessen et al., 2021; The Anopheles gambiae 1000 Genomes Consortium, 2021). The samples from Nigeria and Madagascar were not species-labelled, but they were unambiguously classified by an earlier version of this classifier and we included these samples because Nigeria and Madagascar fill geographic gaps in our sampling dataset. The species represented in GCref v1 are An. gambiae (406), An. coluzzii (222), An. arabiensis (94), An. quadriannulatus (11), An. melas (3), An.merus (6), An. bwambae (3), An. tengrela (38), and An. fontenillei (4). These samples are generally high coverage: 97% of samples have at least 55 of 62 targets amplified. The average number of targets tends to be lower for the samples representing species other than An. gambiae, An. coluzzii, or An. arabiensis, which are also those species represented by fewer samples, but the geographic ranges of these other species are also much more restricted so the samples we do have are likely good representatives of the species. The species represented in GCval v1 are An. gambiae (80), An. coluzzii (15), An. arabiensis (30), An. melas (1), An. merus (5), and An. tengrela (12). The average number of amplicons for the species other than An. gambiae, An. coluzzii or An. arabiensis is lower than in GCref v1 set. Given that for those species, there is only a small number of samples available, we decided to use ones with at least 45 targets amplified in GCref v1 and the ones with at least 30 targets in GCval v1. Sample information for these datasets can be found in Supplementary files 5 and 6.

The input for the VAE is one 8-mer count table per sample, summed over all targets. If a test sample is heterozygous at a given target, we translate each of its haplotypes to an 8-mer count table and sum them to get the test sample’s 8-mer count table for the corresponding target. If a test sample is homozygous at a given target, we translate its haplotype to an 8-mer count table and double the counts to obtain its 8-mer count table at the corresponding target. The counts are doubled in order to represent the target sequences as diploid sequences, and not introduce artificial differences in the total number of 8-mers between homozygous and heterozygous target sites. If the test sample has more than two different haplotypes at a given target, two haplotypes are chosen at random and the sample is treated as a heterozygote. It happens on average in less than 1% of the amplified targets that there are more than two different sequences, so we expect that this inexact way of dealing with those cases does not have a major impact on the results. Because we model the 8-mer counts as the observations of a Poisson distribution, the counts have to be integers, hence we cannot average over all observed alleles as in the nearest neighbour method. If a given target did not amplify in the test sample, the associated 8-mer count table will just contain zeroes, equivalent to simply ignoring missing data. To obtain the VAE input, we sum the 8-mer tables over all 62 targets. This results in a single table per sample, with 65,536(=4⁸) integer entries, roughly summing to twice the number of basepairs covered by the amplified targets, so a little under 20,000 for a sample in which all targets amplified.

The within-complex assignment is based around a variational autoencoder. We considered two other non-linear dimension reduction methods to perform the projection step: UMAP (McInnes et al., 2018) and t-SNE (van der Maaten and Hinton, 2008). t-SNE does not have the ability to embed new samples onto an already existing projection. This means that the distribution of the training set samples depends on the test set we are projecting, which is not desirable for a classification problem. UMAP can distinguish between most species in the An. gambiae complex. The projections have a distinctly different quality from the VAE: the clusters are much tighter and well separated, there are no ‘fuzzy’ boundaries as in the VAE. However, these tight and well separated clusters come at the cost of containing some hard misclassifications; that is samples of one species, which lie well inside the cluster of a different species, which makes it hard to attach any measure of uncertainty to the assignments. For a full description of the UMAP projections and their comparisons with the VAE projection, see Appendix 2. We also ran ADMIXTURE on our data, but it did not separate the species in GCref v1 by assigning them to a unique ancestral population, see Appendix 3.

The VAE consists of an encoder, a latent space projection and a decoder. The specific design we used was inspired by popVAE (Battey et al., 2021). The encoder is a fully connected neural network that takes high-dimensional data as input and encodes that as a point in a latent space of much lower dimension. The decoder is also a fully connected neural network, and it takes as input a point in the latent space and outputs ‘simulated’ data of the same dimensions as the input data. The VAE learns a ‘good’ encoding by adjusting the weights in the encoder and decoder to obtain an output similar to the original input. To prevent overfitting, the input of the decoder is not the exact output of the encoder, but a nearby point in latent space. Furthermore, the loss function used to update the encoder and decoder weights contains a regularisation term on the latent space, in addition to the term measuring the similarity of the decoder output and the original input. Because of the introduced sampling noise and the regularisation constraint, the most efficient way to encode the data is to represent samples that are similar in the high-dimensional data by nearby points in the low-dimensional latent space. In summary, we expect that species identity shapes the structure in the 8-mer count tables and that the VAE projects this structure to the low-dimensional latent space, resulting in clustering by species in the latent space.

In our case, the encoder input is the 8-mer count table of the training set, so a table of dimension n x 65536 with non-negative integer entries, where n is the number of samples in the training set. The output of the encoder is a set of 2d parameters for each sample, where d is the dimension of the latent space. For each dimension, one parameter corresponds to the mean position in latent space and one parameter corresponds to the variance of the position in latent space. In our case, we use a three-dimensional latent space (d=3). The input of the decoder is a position in latent space for each sample, sampled from the distribution determined by the encoder output. The decoder output is an n x 65536 dimensional table of strictly positive entries, however, unlike the input table, the entries are not necessarily integers.

The loss function is the sum of two terms: one measuring the difference between the input and output data and one acting as a regulariser on the latent space. The relative weight of these terms can be adjusted. If we model the count tables as observations of independent Poisson variables, the difference between the input and output can be measured as the Kullback-Leibler divergence (KL divergence) of the Poisson distribution with the means given by the output from the Poisson distribution with the means given by the input. The KL divergence is the same up to a constant as the negative of the Poisson loglikelihood with the input as the observed counts and the output as the means. So minimising the KL divergence is equivalent to maximising the loglikelihood with respect to the output. The difference term of the loss function is obtained by summing over all unique 8-mers. The theoretical minimum of the difference term of the loss function is zero, and this is attained if the output is exactly the same as the input. However, this theoretical minimum cannot be attained in practice, because the sparsity of the input implies that there will be entries equalling zero and the activation function of the decoder generating the output results in strictly positive entries.

The regularisation term of the loss function is based on the KL divergence of the normal distribution parameterised by the encoder output from a standard normal distribution, that is $N\left(\mathrm{0,1}\right).$ The regulariser is computed for each latent space dimension separately. Again, the theoretical minimum equals zero and is attained when the mean outputted by the encoder is zero and the variance outputted by the encoder is one. The regularisation term of the loss function is defined as the KL divergence summed over the latent space dimensions. This effectively enforces the distribution specified by the encoder output to look like a multi-dimensional Gaussian distribution with mean zero, variance one and covariance zero. In layman’s terms, it pulls the projected positions of the samples in latent space towards the origin and establishes a natural scale for them. The regulariser prevents overfitting by making it expensive for the encoder to place samples far away from the origin. The loss function used to train the VAE is the weighted sum of the similarity term and the regularisation term described above, with a parameter w that controls the relative strength of the two terms.

We set most of the parameters involved in training the VAE by comparing the latent-space projections for different parameter values, using a subset of GCref v1 containing only samples from An. gambiae, An. coluzzii, and An. arabiensis. The criteria we used to pick parameter values were species classification accuracy of the reference set and the validation set, using assignments based on convex hulls (described below), visible within-species structure, and, as a secondary criterion, useful visualisation. Further detail on the choice of parameter values is provided in Appendix 4. The latent space projection is also affected by the training dataset. We observed that the presence or absence of most countries does not affect the classification, except for the Gambia and Guinea-Bissau. When these countries are removed, the accuracy to distinguish between An. gambiae and An. coluzzii considerably reduces. This is not surprising, since these samples lie on the boundary of the An. gambiae and An. coluzzii clusters and as such are crucial for assigning those species. The complete results can be found in Appendix 5.

We use the trained and tuned VAE to assign species as follows. We input the summed 8-mer table of the test samples into the encoder of the VAE. The encoder outputs a position in latent space for each sample. Importantly, the VAE is agnostic to species labels; the species assignment happens based on the position in latent space of the test samples in relation to the latent space positions of the species-labelled reference dataset GCref v1.

The top two panels of Figure 4 show the latent space projection of GCref v1. While most species form nicely isolated clusters, An. gambiae and An. coluzzii border each other closely. Interestingly, the boundary is formed by mosquitoes from The Gambia and Guinea-Bissau. These mosquitoes are labelled as An. gambiae by conventional molecular barcoding, but they cannot be confidently assigned to either An. gambiae or An. coluzzii using over 500 ancestry informative markers (AIMs) or whole genome PCA (Anopheles gambiae 1000 Genomes Consortium, 2017; Clarkson et al., 2020). The clusters containing An. bwambae and the putative new species An. fontenillei are placed very close to each other, and can also not be reliably distinguished. These species are closely related, but up until now they have only been discovered in Uganda and Gabon, respectively, and so, since they do not seem to have overlapping geographic species ranges (Barrón et al., 2019), the species identity of samples falling into either of these two clusters can be resolved by their geographic origin.

Figure 4

Download asset Open asset

We perform species classification using the convex hulls of species clusters. A convex hull is the mathematical notion of the smallest convex set containing all points of interest. A nice metaphor is to imagine that you are wrapping all points corresponding to samples of a single species together in such a way that requires the minimal amount of wrapping paper. We constructed one convex hull for each species represented in the dataset, using latent space positions of the samples from GCref v1 as well as of 363 additional species-labelled samples that were not used in training the VAE, to account for possible effects caused by projecting the samples to the latent space (sample information to be found in Supplementary file 7). For our classification procedure, it is important that the convex hulls of different species do not overlap. When constructing convex hulls from the full sample set, only the convex hulls corresponding to An. gambiae and An. coluzzii overlap. We trimmed these hulls by iteratively removing samples from the set of points used to construct them until they did not overlap. In total, we removed 17 An. gambiae and 6 An. coluzzii samples. The samples from An. bwambae and An. fontenillei are combined in one convex hull because they are so close together in latent space.

The classification of new samples happens as follows. If the latent space position of the test sample falls inside a convex hull, the sample is classified as that species. If the latent space position of the test sample falls outside all convex hulls, there are two options. If the sample is much closer to one convex hull than to all others, it is classified as the species corresponding to the hull it is closest to. To be precise, this happens if the euclidean distance to the closest convex hull is at least 7 times smaller than the distance to all other convex hulls. This allows for ‘fuzzy’ boundaries of the convex hulls that are proportional to the separation between the different hulls. We fitted the parameter value 7 on the dataset from Gabon to reflect the assumption that An. tengrela is not believed to be found in Gabon. If the latent space position of the test sample falls outside all convex hulls and outside their fuzzy boundaries, the sample is assigned ‘uncertain’ followed by the labels of all the species whose convex hulls are within a radius of 7 times the distance to the closest convex hull, in order of proximity. This assignment reflects the uncertainty in our classification of samples that fall in an area in latent space where no species-labelled samples fall. At the same time, it gives information on our best guess or guesses for the species identity and leaves open the possibility to modify the assignments based on prior knowledge, for example species ranges or habitat restrictions.

There is one exception to the above assignment procedure: if the closest two convex hulls are from An. gambiae and An. coluzzii and the euclidean distance to both of these hulls is less than 14 then the test sample is classified as uncertain_gambiae_coluzzii or uncertain_coluzzii_gambiae, depending on which convex hull is closer, because we cannot reliably distinguish between these species in this part of the latent space.

In addition to structure driven by species, the latent space projection of GCref v1 also exhibits geographical structure. Within the An. gambiae species cluster, there is a distinct subcluster of samples from Madagascar. Similarly, there is a distinct subcluster of samples from Madagascar in the An. arabiensis species cluster. As mentioned before, the boundary between the An. gambiae and An. coluzzii species clusters is formed by samples from the far West of Africa. These samples also stand out as a separate group in a study on whole genome data (Caputo et al., 2021; Clarkson et al., 2020). There also appears to be a cluster of East African An. gambiae samples that are distinct from the main cluster of An. gambiae as well as the Madagascar samples. It is promising to see signatures of geographic structure within species from the amplicon panel data, because this suggests that the panel can also be useful to explore population structure within species.

We applied the species assignment procedure to GCval v1 (Figure 5). 134 out of 142 samples (94.4%) samples are assigned to a single species and 132 of those (98.5%) are assigned to the species concordant with their species label. One sample labelled as An. coluzzii is classified as An. gambiae and one sample labelled as An. gambiae is classified as An. coluzzii. For all eight samples classified as uncertain, the reference species label is among the set of assigned labels. Seven of the samples classified as uncertain had fewer than 45 targets and we know that the proportion of missing targets affects the position in latent space of the projected samples. The other sample classified as uncertain falls in the space between An. gambiae and An. coluzzii. Further information can be found in Supplementary file 6.

Figure 5

Download asset Open asset

The resulting haplotypes are stored in a table with columns recording the sample name, target, haplotype sequence, read count of supporting reads (for amplicon data only), fraction of supporting reads (for amplicon data only). So each row corresponds to a unique haplotype for a sample target combination, hence samples that are heterozygous at a certain target will have two rows for the same target.

Species labels are assessed and species-groups constructed using custom Python scripts implemented in python 3.8 (Van Rossum and Drake, 2009). The VAE is implemented in keras 2.3 (Chollet, 2015) using a custom python script. Convex hull construction and distance computations rely on scipy 1.6 (Virtanen et al., 2020) and pygel3d 0.2 (Baerentzen, 2018). Plots are created with matplotlib 3.3 (Hunter, 2007) and seaborn 0.11 (Waskom, 2021). All scripts and environments are available on GitHub.

Alignments of amplicon target sequences from highly diverged species are often poor and it is difficult to define a ‘fair’ distance metric based on these alignments. Moreover, there is not a straightforward way to account for small indels and structural variants with alignment based distances. K-mer based distances naturally incorporate indels and structural variation and account for highly diverged sequences in an objective way and provide a solution to the problems arising from relying on alignments. Therefore, our species assignment method uses k-mers to support better comparisons between the sequences in the database and the sequences of the test sample.

There is a trade-off in the choice of k. For large k there is little tolerance for errors, while for small k there is a high chance that the same k-mer is found in multiple locations in the sequence. For example, in a 149 bp sequence, 5 evenly spread SNPs result in no 25-mers matching the reference. On the other hand, the chance that all 4-mers are unique in a sequence of the same length is incredibly small (<10^–22). Based on these trade-offs, we selected 8-mers as a reasonable length. The total sequence length of the amplicon panel targets for the current An. gambiae PEST reference genome sequence AgamP4 is 9928 bp, with a mean target length of 160 bp. There are 65536 unique 8-mers, so the chance that all 8-mers within a target are unique is 84% on average. Across the nearly 10 kb of amplified sequence, the chance that all 8-mers are unique is vanishingly small, but the expected number of unique 8-mers is approximately 8533 (sd 46) and the expected number of non-unique 8-mers is approximately 680 (sd 22).

The methods we present here work with k-mer tables created from each haplotype from each target. A k-mer table consists of 4^k columns, each corresponding to a unique k-mer. To translate a sequence to a k-mer table, we record in each column how often the corresponding k-mer occurs in the sequence. This results in a sparse table (the sparsity of course depends on the choice of k) with non-negative integer entries. As an example, the 2-mer table for the sequences AACTACTCT (first row) and AGCTACTT (second row) is shown below. Note that all possible k-mers are represented in the table, even when they do not appear in any sequence.

The k-mer distance between two sequences is defined as follows. Translate both sequences to k-mer tables as described above and call these q₁ and q₂. Then the k-mer distance between them is given by $d_k\left(q_1,q_2\right)=\frac{\Sigma \left|q_1-q_2\right|}{\Sigma \left|q_1+q_2\right|}$ , that is the number of non-matching k-mers divided by the total number of k-mers in both sequences. The normalisation is required to correct for a bias attributing smaller distances to shorter sequences. In the nearest neighbour method we use this k-mer distance to compare one haplotype sequence of a test individual against one haplotype sequence from the reference database. Note that this definition relies on the k-mer counts, not simply on the presence or absence of each k-mer. Using this definition, the 2-mer distance between AACTACTCT and AGCTACTT is 7/15.

All mosquitoes collected as part of this work were done so under the local guidelines and with permission, which vary depending on the country of origin. ABS (Nagoya) compliance for receiving mosquitoes or their DNA in the UK was completed for all mosquitoes used in this research. Raw sequencing data will be made available on ENA (accession to be confirmed). Pipelines and analysis code, together with processed target haplotypes are available on GitHub: https://github.com/mariloubodde/NNoVAE.

The wild-caught mosquitoes in the reference database were labelled either morphologically or confirmed with species diagnostic PCR. Species identities were further determined by Sanger sequencing the molecular barcodes Cytochrome c oxidase I (COI) and Internal transcribed spacer 2 (ITS2). Both markers were compared to the NCBI database to assess species identity, and COI was also compared to the BOLD database. (See Makunin et al., 2022 for more information).

Some partner labels were revised when distances between the samples in the reference dataset suggested mislabelling. Preferably the relabelling was further supported by molecular barcodes. Relabelled samples are:

• Amar-5: An. marshallii → An. gambiae. The average distance to samples in the An. marshallii complex is 0.54 (0.52–0.58), while the average distance to samples in An. gambiae/coluzzii is 0.09 (0.07–0.13). Additionally, both COI and ITS2 have An. gambiae as best hit.

• Amar-42: An. marshallii → An. jebudensis. The average distance to samples in the An. marshallii complex is 0.38 (0.35–0.42), while the distance to An. jebudensis is 0.05. Unfortunately, neither An. marshallii nor An. jebudensis is present in the barcode databases.

• Amar-3–1: An. marshallii → An_marshallii_cp_sp1. It has the lowest distances to samples of An. hancocki, An. brohieri, An. demeilloni, An. theileri, on average 0.26 (0.22–0.29), while the distances of An. hancocki, An. brohieri and An. demeilloni samples to each other is on average 0.03 (0.02–0.04) and the distances of the An. theileri samples is 0.07 (0.05–0.08). The average distance between the group containing An. hancocki, An. brohieri and An. demeilloni and the An. theileri is 0.29 (0.25–0.33). So it seems like Amar-3–1 does not belong to either of these two groups, but it is as related to these groups as they are to each other. Unfortunately it does not have any informative matches on the molecular barcodes. Therefore this sample is labelled as an unnamed species in the An. marshallii complex.

• Adem-15: An. demeilloni → Myzomyia_sp1. Its lowest distance to another sample in the database is 0.39, to an An. funestus. This is a much higher distance than a sample usually has to other samples of the same or of a closely related species. In particular, the distance of An. hancocki, An. brohieri and An. demeilloni samples to each other is 0.03 (0.02–0.04), while the average distance of Adem-15 to these samples is 0.41 (0.39–0.44). As it does not show similarity to any of the samples in the database, but it is closer to almost all samples in the Myzomyia series than to almost all samples outside the Myzomyia series, it is labelled as an unnamed species in the Myzomyia series.

• Acol-645: An. coluzzii → An_gambiae_cp_sp1. This sample is different from other An. coluzzii and An. gambiae, this was found both in the sample-pair differences and in the subsequent VAE method which is tailored at the An. gambiae complex. The molecular barcodes also suggest that this sample is not An. coluzzii or An. gambiae, but that it is a member of the complex. Therefore, it is labelled as an unnamed species in the An. gambiae complex.

• Apal-257: An. paludis → An_coustani_cp_cl3. Its smallest distances are to samples of An. tenebrosus, An. ziemanni, An. coustani and An. paludis, on average 0.29 (0.28–0.32). We found that the samples of these species split into two clades and the distance of Apal-257 to either of these clades was much higher than the distances between samples in the same clade. Additionally, the average distances between the two clades is 0.18 (0.14–0.21), so Apal-257 is also further away from the two clades than they are from each other. Yet, it is closer to these two clades than it is to An. sinensis and An. hyrcanus, which are the next closest species. Regarding the molecular barcodes, Apal-257 matches to An. coustani on COI, like most samples in the An. coustani complex. On ITS2 it matches to An. junlianensis and An. yatsushiroensis, species in the hyrcanus group. We therefore decided to relabel it to a third clade in the An. coustani group.

• An. nili samples. Anils-7: An. nili s.s. → An_nili_gp_sp1. Anil-237 & Anil-239: An. nili → An_nili_gp_sp2. Anil-233, Anil-236 & Anil-238: An. nili → An_nili_gp_sp3. Both the distances and alignment of the molecular barcodes support this split. It is possible that some of these samples represent different member species of the An. nili group, which have been shown from cytogenetic analysis to differ substantially (Sharakhova et al., 2013), but molecular barcodes are not yet available in public databases for all member species.

• An. hyrcanus samples. VBS00085 & VBS00086: An. hyrcanus → An_hyrcanus_gp_sp1. VBS00082 & VBS00083: An. hyrcanus → An_hyrcanus_gp_sp2. The distances within these pairs are 0.08 for both pairs. The distances between the pairs are 0.35 (0.34–0.36). The molecular barcode matches support the split. Even though hyrcanus is present in all databases, there are few matches to it. On COI in both BOLD and NCBI, the first pair matches to nitidus and the second pair to crawfordi; these species are in two distinct subgroups of the hyrcanus group. On ITS2 the first pair matches to hyrcanus and then nitidus, and the second pair matches to sinensis (sinensis is in the main hyrcanus group, not in either of the aforementioned subgroups). The alignments of COI and ITS2 sequences for these four samples also clearly support the split into pairs.

• Samples in the An. coustani complex. This complex contains the species An. tenebrosus, An. ziemanni, An. coustani and An. paludis. The between sample distances indicate a split into two clades (disregarding sample Apal-257, discussed above); however, the split is not correlated with the species labels. The proposed clades are An_coustani_cp_cl1 containing Aten-191, Aten-185, Azie-334, Acou-956, Acou-959, Acou-962 and An_coustani_cp_cl2 containing Aten-333, Aten-79, Aten-954, Azie-1032, Azie-1055, Azie-70, Azie-77, Acou-71, Acou-80, Apal-81. The average distances between members of the same clade is 0.07 (0.05–0.09) and 0.07 (0.05–0.10) respectively, while the distance between members of different clades is 0.18 (0.14–0.21). The barcode matches are to An. coustani for the vast majority of the samples, even though COI sequences for An. ziemanni and An. tenebrosus are present in both BOLD and NCBI database. Alignment of the COI sequences shows extremely little variation. Alignment of ITS2 does support a split into the two proposed clades.

All labels and groupings are displayed in Supplementary file 1. In the majority of cases, the fine species-groups are supported by at least one molecular marker and not contradicted by different species labels. Some species are not present in one or both databases, so for these a match to a closely related species is allowed. The exceptions are:

• An. jebudensis: one sample was labelled as An. jebudensis, the other one originally as An. marshallii, but showed sufficient evidence for relabelling. Neither An. jebudensis nor An. marshallii is present in either database. On ITS2 they both match to An. moucheti, which is the closest species in the dataset, neither has a match in BOLD, and on COI in the NCBI database the best matches are to An. lindesayi and An. sawyeri, which are in the Anopheles and Nyssorhynchus subgenus respectively. However, it has to be noted that the NCBI database does not contain an An. moucheti COI sequence and the best hits are based on 81% and 91% identity respectively.

• An. rhodesiensis: All samples match to uninformative ‘Anopheles_sp’; even though a COI sequence is available in BOLD and NCBI. The best hit to a named species in the dataset is to An. funestus for COI and An. aconitus for ITS2.

• An. jamesii: Is predicted for BOLD and ITS2, but not for COI in the NCBI database, even though an An. jamesii COI sequence is present there.

• An. maculatus A: has a few matches to An. sawadwongporni, which is also a member of the maculatus group, even though sequence of An. maculatus A is present in all databases.

• An. rampae: the four different samples match to An. rampae(2 x), An. sawadwongporni(1 x) and An. maculatus(1 x) in the BOLD database; all these are in the same species group. In the NCBI database, all COI hits are for An. maculatus (An. rampae COI not present in database), and all ITS2 hits are An. rampae and the uninformative Anopheles_sp.

• An. balabacensis: matches to An. introlatus on COI in both databases, even though An. balabacensis sequence is present. The two species are in the same complex.

• An. carnevalei: all samples match to An. carnevalei on BOLD and ITS2. An. carnevalei sequence is not present in NCBI database for COI; four out of five samples match to An. nili, one to An. darlingi, which is in a different subgenus.

• An. dureni: sequence not present in either database; on ITS2 one sample matches to An. minimus and one to An. leesoni, which are both in different groups; but placement of An. dureni is not very clear on the species tree.

• An. vinckei: sequence not present in either database; on COI there are distant matches to An. gambiae and An. maculatus, which are both in different groups; but placement of An. vinckei is not very clear on the species tree (seems to be close to An. dureni).

• An. coustani group: we have samples of four species, An. coustani (present in all databases), An. paludis (present in no databases), An. ziemanni (COI present in both databases), and An. tenebrosus (COI present in both databases). In BOLD, all samples match to An. coustani. In the NCBI database, on COI most matches are to An. coustani, there is one match to An. ziemanni, and a few matches to Anopheles_cf. On ITS2 there are many matches to Anopheles_sp. and Anopheles_cf., and the remaining matches are to An. coustani and one sample to An. junlianensis and An. yatsushiroensis in the hyrcanus group (this is the An. paludis sample, that is further removed from all other samples in this complex). These samples also show a ‘checkerboard’ pattern, which splits them into two clades, which are not correlated with species labels.

• An. barbirostris: there are matches to An. barbirostris and to the closely related An. dissidens. However, the samples do not form two separate groups based on the distances.

• An. oryzalimnetes: no Sanger sequencing done.

• An. cruzii: no Sanger sequencing done.

• An. bellator: no Sanger sequencing done.

Lastly, there are samples which are further than 0.10 distance away from other samples in the reference database representing the same species. Above we have discussed some samples where there was good evidence to adjust the species label, but there are also some for which there is good reason to retain the partner label. Those are:

• VBS00001: An. annularis. It is a bit more diverged from the other three An. annularis samples, which causes it to fall just above the threshold. However, the distances are not very different from the distance between the other An. annularis, the next closest sample in the reference set are much further away and the molecular barcodes support the partner label.

• Agam-37: An. gambiae. It is a bit more diverged from the other An. gambiae that were sequenced by the panel, which causes it to fall just above the threshold. But is it closer to An. gambiae than to other samples in the reference dataset. The molecular barcodes also support the partner label.

• VBS00149 and VBS00150: An. tessellatus. They are 0.11 distance from each other, but for both the molecular barcodes match to An. tessellatus and the next closest sample in the database is at 0.44 distance.

• anopheles-sinensis-chinascaffoldsasinc2 & anopheles-sinensis-sinensisscaffoldsasinc2: An. sinensis. They are 0.11 distance from each other and that is still clearly closer than the next best match in their sister species An. hyrcanus (namely An_hyrcanus_gp_sp2).

We projected the An. gambiae complex training set (GCref v1) with UMAP, using the same 8-mer count tables as for the VAE as input. Using the default UMAP parameters does a very poor job at distinguishing species (Appendix 2—figure 1). By experimenting with different parameter values and different metrics, the three best-performing projections (judging by eye after colouring samples by their species label, see Appendix 2—figure 1) are:

• metric = hamming, n_neighbors = 3, min_dist = 0, n_components = 2

• metric = canberra, n_neighbors = 3, min_dist = 0.1, n_components = 3

• metric = braycurtis, n_neighbors = 3, min_dist = 0.99, n_components = 3

Appendix 2—figure 1

Download asset Open asset

We have two species that are only represented by three samples and by setting the parameter n_neighbours higher than 3, these species will not split out as separate clusters. Although most of the species form distinct clusters in these three projections, there is still significant overlap between An. gambiae, An. coluzzii and An. tengrela.

It is noticeable that UMAP tends to generate tighter clusters than the VAE (see Figure 4 in the main text); there are less fuzzy boundaries, but it comes at the cost of some samples being placed in the wrong cluster. The convex hull assignment we used with the VAE projection is specifically aimed to deal with the fuzzy boundaries and it doesn’t cope well with outlier samples. Therefore, we use a different classification method, namely support vector machines, to perform species assignment using the UMAP projections. We trained an SVM with the default parameters on the reference dataset GCref v1. The accuracy results and comparison are shown in the table below.

We suspected we could improve the results by treating separate clusters of the same species as their own assignment category, for example relabel the samples in the separate An. arabiensis cluster as An_arabiensis-2. However, this did not improve the accuracy, neither for the reference set GCref v1 nor for the validation set GCval v1.

Overall the VAE has a lower misclassification score, although this comes at the cost of a larger percentage of samples assigned to multiple species. We have consciously decided for a more conservative classifier and we intend to resolve the samples which cannot be assigned to a single species on an ad hoc basis; either by more in depth analysis of specific amplicons (to be developed), by considering carefully documented species ranges where available (e.g. to distinguish between An. bwambae and An. fontenillei, both of which have very specific geographical areas where they occur) or by performing whole genome shotgun sequencing if the identity of the samples is of greater interest.

The UMAP projections can probably be tuned better to achieve similar assignment performance to the VAE, but UMAP does not outperform the VAE.

To assess the amount of population differentiation between the different species in the An. gambiae complex, we ran ADMIXTURE on the reference dataset GCref v1. We ran ADMIXTURE on the 8-mer count tables; first we filtered for variable 8-mers which take count values in {0,1,2} for all samples; there are 25,301 such 8-mers. We treated these 8-mers essentially as unlinked SNPs. We are aware that these input data do not exactly follow the assumptions behind the probabilistic model used by ADMIXTURE, but using only unlinked SNPs would allow for at most one SNP per amplicon, so 62 in total and it would be quite challenging to select one without bias to a certain species or group of species. Despite not following the model assumptions, we believe that ADMIXTURE used on k-mers should allow us to judge whether ADMIXTURE constitutes a useful approach to differentiate between species in the An. gambiae complex.

The inferred ADMIXTURE proportions for different values of K can be found in Appendix 3—figure 1. At K=3 and K=4, ADMIXTURE does a good job at distinguishing the well-represented species: An. gambiae, An. coluzzii and An. arabiensis. However, it struggles to differentiate between the species represented by a smaller number of samples, even though many of these species are much more diverged from each other than An. gambiae and An. coluzzii are (Fontaine et al., 2015).

Appendix 3—figure 1

Download asset Open asset

For the higher values of K, ADMIXTURE reveals substructure and mixed ancestry in the An. gambiae and An. coluzzii samples. The first glimpse of substructure we see already at K=3, where the An. gambiae samples from Guinea-Bissau and the Gambia are modelled as a mixture of An. gambiae and An. coluzzii. This substructure is also found in the VAE projection.

While the An. arabiensis from Madagascar split out as a separate cluster in the VAE and in the three UMAP projections discussed in the next section, ADMIXTURE only identifies them as a separate group at K=11. So similarly as for the species represented by few samples, it seems like ADMIXTURE is hesitant to assign a unique ancestry to a small group of samples.

ADMIXTURE identifies some structure in the An. gambiae complex when treating variable 8-mers as SNPs. It is particularly good at differentiating between groups represented by a larger number of samples and additionally finds some substructure within these groups. It is pronouncedly worse than the VAE at differentiating between species for which we have fewer samples. This is a shortcoming given that a major aim of applying ANOSPP widely to samples from across the globe is to identify novel species that may or may not be contributing to transmission.

To assess which values of K are most informative, we recorded the cross-validation error for these runs, see Appendix 3—figure 2. There are 9 species represented in the dataset GCref v1, so ideally we would find a minimal cross-validation error around K=9. The cross validation plot indicates that we need a minimum of about 9 populations, but it does not unambiguously indicate a suitable value of K. To properly assess the cross validation error, we’d have to run multiple replicates.

Appendix 3—figure 2

Download asset Open asset

Although ADMIXTURE does find structure in the An. gambiae complex, it is principally a method for assessing population differentiation, not for species assignment. It is possible to build an assignment method around it, but given that it performs worse than the VAE at differentiating between all species in the An. gambiae complex, we prefer using the VAE based method for assignments.

Our VAE architecture is inspired by popVAE (Battey et al., 2021). We used their default parameter choices for learning rate, training iterations and validation proportion as well as the ‘elu’ propagation between the layers of the neural network and the ‘linear’ activation function to transform the output from the encoder to the latent space variables. Because we model k-mer counts as independent Poisson variables, we use the ‘softplus’ activation function to transform the output from the last layer of the decoder.

To set the width and depth of the encoder and decoder, we ran a grid search over a combination of widths and depths. We kept the width and depth the same for the decoder and the encoder and we trained the VAE with three different seeds for each width-depth combination. For high width values (more than 500 nodes) the run time required to train the VAE was very long and the process would often crash when conducted on a desktop. For very small networks (width 32, depth 4) the resolution was poor, but in the middle range the visible structure was not severely affected by the choice of width and depth, so we went with width 128 and depth 6, which gave a good resolution and reasonable run time.

As mentioned in the main text, the loss function contains a parameter, w, that controls the relative importance of the data driven term and the regularisation term. We ran a grid search ranging over four orders of magnitude, training with three different seeds for each search point. For small values of w, which gives higher weight to regularisation term, the resulting latent space dimension shows a strong correlation between the latent space dimensions, suggesting that only highly differentiated samples can overpower the effect from the regularisation term. For high values of w, the training process of the VAE becomes more unstable: it crashes sometimes and the resulting latent space projections look less similar to those obtained with the same w value, but a different seed. We chose w equal to 1000, which resulted in reasonably stable projections, with good visible structure.

To assess the effect of using 8-mers over k-mers of other length, we ran the VAE with 6-mers and with 10-mers as input, see Appendix 4—figure 1. The projections using 6-mers, in both two and three dimensions, show less tight and more overlapping clusters by species compared to the 8-mer projection. There are only 4,096 unique 6-mers, compared to >65,000 unique 8-mers, and with almost 10 kB of sequence, there will be much fewer unique 6-mers than unique 8-mers. Although the VAE loss function is based on counts rather than presence/absence of k-mers, it probably is more sensitive to the difference between 0 and 1 than to the difference between 2 and 3. The VAE based on 10-mers shows clusters quite similar to those found in the 8-mer projection, although An. quadriannulatus is much closer to An. bwambae and An. fontenillei. It took much longer to train the 10-mer projection (5.5 hr compared to 20 min for the 8-mers), which is not a major drawback, because the training has to happen only once, but it does make it harder to tweak the parameters for the 10-mer projection.

Appendix 4—figure 1

Download asset Open asset

To determine the number of latent space dimensions, we ran a grid search over 2, 3 and 4 latent space dimensions. In fact, the grid search for w and the number of latent space dimensions was run simultaneously, but there seemed to be little interaction between w and the number of latent space dimensions. For the training set containing only An. gambiae, An. coluzzii and An. arabiensis, adding a third or fourth dimension did not add to the visible structure, so in this case we opted for two latent space dimensions for visualisation purposes. However, for the full training set GCref v1, in the two-dimensional projection, several species clusters are very close to each other, making this hard to use for species assignment. Adding a third dimension resolves this problem, even though the second and third dimensions are strongly correlated (see Figure 4). The benefit of the strong correlation is that we can visualise the structure in two dimensions.

Here we test the robustness of our VAE projection by removing all samples collected in one geographic location from the training set and then projecting those samples and the validation set and comparing these results to those obtained using the VAE trained on the full training set. In particular, we will pay attention to the separation of the species clusters, the classification accuracy of the validation set and the visible structure within the species clusters.

For this analysis, we use a subset of the reference dataset GCref v1 and the validation set GCval v1, containing only the species An. coluzzii, An. gambiae and An. arabiensis, because for the other species in the An. gambiae complex we have very few samples, the samples in this dataset are collected at only one or two locations and the species have more limited geographic ranges than An. coluzzii, An. gambiae and An. arabiensis. With only these species, two latent space dimensions are sufficient to exhibit the relevant structure, so that is what we will use here. We have performed these drop-outs for a subset of the geographic locations, selecting for those we thought would be most likely to affect the structure of the projection.

Below, we discuss the drop-out experiments in more detail, but to summarise: overall, the separation of the species clusters remained intact and the ability to classify species of the samples in the training set and in the validation set was not significantly affected. Samples from The Gambia and Guinea-Bissau could be less accurately classified when no samples from this geographic region were included in training the VAE, but for all other locations we tested, the ability to classify species remained the same. The degree to which geographic structure is visible within the species clusters was somewhat affected for certain collection countries. However, the visible structure never completely disappeared, nor did we see a considerably different geographic structure, it was merely that the structure became a bit more blurred or more condensed than in the original projection.

We removed 10 An. coluzzii samples from the training set, see Appendix 5—figure 1. The separability of species clusters is similar compared to the projection of the VAE trained on the full dataset. The geographic structure within the An. gambiae and An. arabiensis clusters remains intact, but the An. coluzzii cluster loses some of its visible structure. Angola is on the edge of the geographical range of the An. coluzzii samples in this dataset and these samples are also on the very top of the VAE projection, suggesting that they are responsible for a considerable amount of the structure within An. coluzzii.

Appendix 5—figure 1

Download asset Open asset

We removed 10 An. coluzzii and 66 An. gambiae from the training set, see Appendix 5—figure 2. The separability of the species clusters is similar compared to the projection of the VAE trained on the full dataset. The geographic structure within the An. coluzzii and An. arabiensis clusters remains intact, but the An. gambiae cluster loses some visible structure in its main subcluster, while the structure in the other four subclusters (formed by samples from The Gambia and Guinea-Bissau; by samples from Nigeria and Gabon; by samples from Madagascar and by samples from various East African countries respectively) remains largely intact. This is probably due to the fact that the samples from Cameroon form a large proportion of the samples in the main An. gambiae cluster.

Appendix 5—figure 2

Download asset Open asset

We removed 33 An. coluzzii and 19 An. gambiae from The Gambia and 24 An. gambia from Guinea-Bissau from the training set, see Appendix 5—figure 3. Considering the training samples, the separability of the species clusters is very clear. However, the classification accuracy of the projected samples from The Gambia and Guinea-Bissau is only 70%. This is probably due to the fact that these samples form the border between the An. coluzzii and An. gambiae clusters in the projection of the VAE trained on the complete dataset. So by leaving those samples out, we miss the important label information necessary to classify the edge cases. The subclusters of the An. gambiae cluster are arguably more pronounced than in the original projection, but the structure within the subclusters of the An. gambiae cluster and in the entire An. coluzii cluster reduces.

Appendix 5—figure 3

Download asset Open asset

We removed 20 An. arabiensis and 18 An. gambiae from the training set, see Appendix 5—figure 4. The separability of the species clusters is similar compared to the projection of the VAE trained on the full dataset and the geographic structure within the three species clusters remains largely intact. When the Madagascar samples are included in training the VAE, they form a tight subcluster, both for An. gambiae and An. arabiensis. When they are not included in training, these clusters stand out less, due to the fact that the VAE does not recognise the features that make them stand out.

Appendix 5—figure 4

Download asset Open asset

We removed 57 An. coluzzii and 43 An. gambiae from the training set, see Appendix 5—figure 5. The separability of the species clusters is similar compared to the projection of the VAE trained on the full dataset. The visible geographic structure within the three species clusters reduces a bit. This could be because we removed such a large number of samples.

Appendix 5—figure 5

Download asset Open asset

We removed 6 An. coluzzii and 20 An. gambiae from the training set, see Appendix 5—figure 6. The separability of the species clusters is similar compared to the projection of the VAE trained on the full dataset and the geographic structure within the three species clusters remains largely intact. When the Nigerian samples are included in training the VAE, the An. gambiae individuals are more separated from the main An. gambiae cluster than when they are projected using the VAE excluding them from training. This is because the VAE has not been trained to recognise the features that distinguish these Nigerian An. gambiae individuals from the other samples.

Appendix 5—figure 6

Download asset Open asset

We removed 37 An. arabiensis and 20 An. gambiae from the training set, see Appendix 5—figure 7. The separability of the species clusters is similar compared to the projection of the VAE trained on the full dataset. The subclusters of An. gambiae and An. arabiensis remain well separated, but the visible geographic structure within the (sub)clusters reduces a bit.

Appendix 5—figure 7

Download asset Open asset

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

We thank Vickie Brooks for general logistic support for ANOSPP. We thank Petra Korlević for support in wet lab work, valuable discussions and conducting the mitochondrial analysis. For field sampling in Burkina Faso, we would like to thank Franck Adama Yao, Patric Stéphane Epopa Ngomé, Nouhoun Traoré, and the Target Malaria Team. We thank Samuel Anifowose, Akinkunle Adeniyi, Joshua Oduwu, Idris Otun, Juwon Adeniji, Wale Enisemo, Cynthia Umunnakwe and Adedotun Bayegun for sample collection and processing in Nigeria. In Gabon, we thank the ESV team at the CIRMF. We also thank the Institut Pasteur de Madagascar for sample collection in this country. We thank Sanger’s Scientific Operation Teams for carrying out all PCRs, library generation, and sequencing on data presented here and Catherine McCarthy for her support in ensuring all samples are compliant with Access and Benefit Sharing of sequence data as laid out by the Nagoya Protocol. The Wellcome Sanger Institute is funded by the Wellcome Trust (206194/Z/17/Z), which supports MKNL and the work contained here. MB was supported by Wellcome 4 year PhD studentship (RG92770), and RD by Wellcome award (WT207492). Gabon collections were supported by an ANR grant to DA (ANR-18-CE35-0002-01 – WILDING). LB was supported by an IRD student fellowship (Bourse ARTS/IRD).

© 2022, Boddé et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.