Genes linked to schistosome resistance identified in a genome-wide association study of African snail vectors

Abstract

Schistosomiasis, a neglected tropical disease, is transmitted by freshwater snails. Interruption of transmission will require novel vector-focused interventions. We performed a genome-wide association study of African snails, Biomphalaria sudanica, exposed to Schistosoma mansoni in an endemic area of high transmission in Kenya. Two snail genomic regions, SudRes1 and SudRes2, were significantly associated with snail resistance to schistosomes. SudRes1 includes receptor-like protein tyrosine phosphatases while SudRes2 includes a class of leucine-rich repeat-containing G-protein coupled receptors, both comprising diverse extracellular binding domains suggestive of host-pathogen interaction. Resistant and susceptible haplotypes show numerous coding differences including presence/absence of entire genes. No loci previously tied to schistosome resistance in a neotropical snail species showed any association with compatibility suggesting that loci involved in the resistance of African vectors are distinct. Snail ancestry was also strongly correlated with parasite compatibility. These results will inform future efforts to predict and manipulate immunity of a major schistosome vector.

Fig. 1. Results of the pooled genome-wide association study (pooled-GWAS) identifying variants associated with resistance to Schistosoma mansoni in the Biomphalaria sudanica genome.

Fisher’s exact test p values for all pooled-GWAS variants are arranged horizontally based on contig orthology to 18 chromosomes (x-axis labels) of the B. glabrata genome (xgBioGlab47.1, NCBI RefSeq: GCF_947242115.1) and linkage map analysis (Dataset S2, Figs. S8, S11, S16). Pooled-GWAS dual-variants, defined as two or more proximate (<50 kb and >1.5 kb apart) variants strongly associated with S. mansoni resistance (uncorrected two-sided p ≤ 2.5e-9) are red. All others are blue. All dual-variants and singleton-variants with uncorrected two-sided p ≤ 1e-17 are labeled red and blue, respectively, with their contig and contig position. Unassigned contigs could not be unambiguously mapped and are mostly small and/or repetitive.

Fig. 2. Amplicon panel validation of pooled-GWAS results.

a Proportion of ancestry from Biomphalaria sudanica Population 1 (red) or Population 2 (blue) as predicted by ADMIXTURE framed on a subset of unrelated individuals (represented by gray dots) is correlated with Schistosoma mansoni infection [Logistic regression, effect on odds of infection = 1.79 + /− 0.28, uncorrected two-sided p = 8e-11; n = 496 biological replicates (positive n = 260, negative n = 236)], such that >90% Population 1 snails (Group A) are 73% positive while >90% Population 2 snails (Group B) are 26% positive. Ancestry analysis sample size was n = 503 including 496 genotyped-validation and genotyped-pooled-GWAS B. sudanica and seven false negative B. sudanica. b Ordered Fisher’s exact test p values per variant of genotyped-validation samples within ancestry groups. SudRes1 variants c582:65,596 and c6844:72,816 are significant validated-variants (uncorrected two-sided p = 5e-05 and 9e-05, respectively; Bonferroni adjusted significance threshold shown by dotted line) within ancestry Group B (n = 34 biological replicates, blue dots) and all other top-outliers (1e-04 <p ≤ 1e-03) are also in SudRes1. c Ordered regression p values per variant of genotyped-validation samples, after accounting for ancestry (n = 220 biological replicates). SudRes2 variant sc94:2,166,296 is a significant validated-variant (logistic regression, uncorrected two-sided p = 2e-06, Bonferroni adjusted significance threshold shown by dotted line) and most other top-outliers (1e-04<p ≤ 1e-03) are in SudRes1 or SudRes2. d Allele and genotype counts for the most significant validated-variants from SudRes1 (c6844:72,816) and SudRes2 (sc94:2,166,296), which were both dominant markers, and the representative codominant marker variants from SudRes1 (c582:65,596) and SudRes2 (sc94:2,174,117), for S. mansoni infection positive (+) and negative (-) snails. Values are read counts for each allele in the pooled-GWAS (“GWAS”), and genotype counts in genotyped-validation snails (“Validation”). e Means and standard errors for each variable in the best-fitting multiple regression model (shown in box). Ancestry (A) is the proportion of Population 2 ancestry, and loci SudRes1 (SR1) and SudRes2 (SR2) are dominant effects at the representative codominant markers c582:65,596 and sc94:2,174,117. **p < 0.01, ***p < 0.001. Logistic regression, uncorrected two-sided p = 3e-07, 5e-04, and 2e-03, respectively; n = 380 biological replicates (positive n = 191, negative n = 189).

Fig. 3. Characterization of Biomphalaria sudanica SudRes1 genomic region.

a Pooled-GWAS uncorrected two-sided p values (Fisher’s exact test) in SudRes1 regions (dashed boxes), which contain pooled-GWAS dual-variants (red) and other variants (blue), defined as in Fig. 1. Contigs (gray rectangles) on B. sudanica chromosome 5 arranged horizontally based on contig orthology to 18 chromosomes of the B. glabrata genome (xgBioGlab47.1, NCBI RefSeq: GCF_947242115.1) and linkage map analysis (Dataset S2, Figs. S8, S11, S16). Gene positions are shown (yellow/brown/green/orange boxes), highlighting particularly prevalent classes of genes: in SudRes1 the multiple epidermal growth factor (MEGF) and galactose-binding like domain (GBD) containing receptor-like tyrosine-specific protein phosphatase (RPTP), other protein coding genes containing MEGF domains. b The predicted protein structure of receptor-like tyrosine-specific protein phosphatase (RPTP) coding gene BSUD.17727 (c6844) present in the SudRes1 region of the B. sudanica genome and containing intronic validated-variants (Fig. S12A). A similar RPTP coding gene is contained within adjacent contigs c582 and c5209 within SudRes1, and contig c1041 neighboring SudRes1 region (Fig. 3A). c Nucleotide alignment of paralogous contigs representing a portion of the SudRes1 region in both Bs111 and Bs2280 genomes, showing location of genes and coding regions relative to GWAS variants used in the amplicon panel. Four G/T polymorphisms are shown, two of which act as non-Mendelian dominant markers due to paralogous amplification, one of which shows two alleles with Mendelian segregation (c582:65,596), and one of which segregates in our population but is invariant in these sequenced genomes (c5209:45,695). Exons are colored by protein domain, and mean pairwise sequence identity calculated in Geneious v2022.0.2 (Biomatters Ltd.) is shown in 100 bp sliding windows across aligned contigs (green = 100% identity, brown = 30% to <100%, red = <30%). d Amino acid alignment of RPTP genes in SudRes1, demonstrating extracellular diversity and loss of functional EGF domains in Bs2280 contig c32471. Protein domains are indicated by color, and non-synonymous polymorphisms in comparison to the majority consensus of paralogs/orthologs are shown in black. Sequence identity calculated in Geneious v2022.0.2 (Biomatters Ltd.) is shown as in (c), for 25 aa sligning windows in aligned proteins.

Fig. 4. Characterization of Biomphalaria sudanica SudRes2 genomic region.

a Pooled-GWAS uncorrected two-sided p values (Fisher’s exact test) in SudRes2 regions (dashed boxes), which contain pooled-GWAS dual-variants (red) and other variants (blue), defined as in Fig. 1. Contigs (gray rectangles) on B. sudanica chromosome 6 are arranged horizontally based on contig orthology to 18 chromosomes of the B. glabrata genome (xgBioGlab47.1, NCBI RefSeq: GCF_947242115.1) and linkage map analysis (Dataset S2, Figs. S8, S11, S16). Gene positions are shown (yellow/red/orange boxes), highlighting particularly prevalent classes of genes in SudRes2 encoding a class of leucine-rich repeat containing G-protein couple receptors (GRL101) with C-type lectin and low-density lipoprotein extracellular domains (partial GRL101 genes included), and a Zinc-RING finger and inhibitor of apoptosis containing protein. b A representative predicted protein structure of a GRL101-like G-protein coupled receptor coding gene, twelve of which were predicted through manual annotation within the SudRes2 region of contig sc94 (1.82–2.26 Mb) in the B. sudanica genome. c Dot plot constructed using D-GENIES comparing synteny of the SudRes2 region in Biomphalaria sudanica Bs111 genome and the Bs2280 resistant snail genome, highlighting regions of divergence and structural rearrangements between the two genomes. GRL101 genes are indicated along with Zn-RING-IAP gene BSUD.25704 and its paralog.