Determining the genetic structure and diversity of Plasmodium malariae populations in sub‐Saharan Africa from microsatellite variation and single nucleotide polymorphisms associated in homologous antimalarial drug resistance genes (Part 2)​​

Targeted amplicon sequencing of antimalarial drug resistance genes

Eleven genes of P. malariae homologous to antimalarial drug resistance genes P. falciparum were amplified and sequenced (Table 2). Specific primers were designed for each gene and the optimal conditions for amplification using the Q5 polymerase (New England Biolabs) are shown in Table 3. Amplification products were checked on a 1% agarose gel and all amplicons from each sample were pooled to prepare a deep sequencing library using the TruSeq HT library prep kit. The library was cleaned using the Agencourt AMPure XP PCR purification kit (Beckman Coulter, Brea, CA, USA) according to the manufacturer’s instructions. Amplicon concentrations and sizes were measured using a Qubit fluorometer (In vitrogen, Carlsbad, CA, Mỹ) and system Tapestation (Agilent).

Paired‐end sequencing was performed on the MiSeq system (Illumina, San Diego, CA, USA) with 10 pM concentration of pooled amplicons, at the MRCG genomics platform, using the Illumina v2 reagent kit to generate 250 bp reads per end, according to the manufacturer’s instructions.

Table 2. List of homologous antimalarial drug resistance genes P. malariae has been sequenced

Table 3. Optimal conditions for generating drug‐resistant amplicons of P. malariae

Q5 Polymerase Enzyme Buffers Available 2.0 mM MgCl₂.

DATA ANALYSIS

Microsatellite or SSR (Simple Sequence Repeat) analysis

The clustered microsatellite alleles were imported as genind objects into the statistical software R (Version 4.1.13) and used for population genetic analysis. Initially, the populations were defined as the countries from which the isolates were collected. Genetic diversity at the population level was assessed based on expected heterozygosity and the number of alleles per locus (i.e., allelic richness). Expected heterozygosity values ​​ranged from 0 to 1 (0 indicating no diversity and 1 indicating all alleles are different). Diversity parameters such as observed and expected multilocus genotypes, Shannon‐Wiener multilocus genotype diversity index, Stoddart and Taylor diversity index of multilocus genotypes, Simpson index, Nei’s unbiased genetic diversity, and evenness were calculated using the “poppr” function in R.

Bruvo pairwise genetic distances were calculated using microsatellite markers for all isolates using the “bruvo.dist” command in R and displayed as a hierarchical heatmap. To determine population structure, the optimal number of initial genetic clusters was determined by running successive K‐means with increasing k values ​​and comparing different clustering solutions using the Bayesian Information Criterion (BIC) in R.

Additionally, the ‘find.clusters’ function in the adegenet package of R (version 2.0) was used to attach each individual to genetic clusters. Discriminant analysis of principal components (DAPC) was applied to describe clusters of genetically related individuals and visualized by scatterplots. DAPC transformed the data using principal component analysis (PCA), and then performed discriminant analysis on the retained principal components using cross‐validation. Ancestry was determined using an admixture model using STRUCTURE version 2.3.4, with individuals attached to K populations based on their multi‐locus genotypes.

The expected K value was set from 1 to 9 and STRUCTURE was run with 100,000 MCMC iterations for 10,000 Burn‐in cycles. The best fitting K value was determined by ΔK which was calculated and implemented on the Structure Harvester Web v0.6.94 tool. The ancestor proportion distribution was then displayed as a bar graph.

Sequence analysis of candidate drug resistance genes and SNPs

The quality of the obtained Fastq files was checked using FASTQC software, then trimmed and filtered to remove poor quality fragments and unwanted indices. The reads were processed as illustrated in Figure 1 ‐ the reads were either mapped directly to the reference set of paired target genes from PLASMODB, or filtered to remove missing data before mapping to the reference set, or processed through a denoising procedure (dada2 version 1.16.0) before performing reference filtering and mapping. The missing data filtering process was performed in three steps:

(i) Samples with more than 80% missing data were excluded from the original dataset;

(ii) Locus with more than 70% missing data from the dataset obtained after the first filtering step (filter 1) is removed;
(iii) Samples with more than 40% missing data from the dataset obtained after the second filtering step (filter 2) were also discarded.;
Details are summarized in Supplementary Table 1.