Immune targets for schistosomiasis control identified by 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 immunity 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 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.

Figure 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, Fig. S2-S4). Pooled-GWAS dual-variants, defined as two or more proximate (<50 kb and >1.5 kb apart) variants strongly associated with S. mansoni resistance (p ≤ 2.5e-9) are red. All others are blue. All dual-variants and singleton-variants with 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.

Figure 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 (31) is correlated with Schistosoma mansoni infection (p < 1e-14), such that majority-Population 1 snails (Group A) are 68% positive while the majority-Population 2 snails (Group B) are 41% positive. Total number of samples included in analysis was 503 including 496 genotyped-validation and genotyped-pooled-GWAS B. sudanica from this study, 5 inbred B. sudanica previously sequenced (27) and 2 outbred B. choanomphala. (B) Ordered Fisher’s exact test p values per variant of genotyped-validation samples within ancestry groups. SudRes1variants c6:3,490, c6:31,057, and c6844:72,816 are significant validated-variants (p < 1.3e-04, Bonferroni adjusted significance threshold shown by dotted line) within ancestry Group B (blue dots) and most other top-outliers (1.3e-04< p < 1.0e-03) are in SudRes1 or SudRes2. (C) Ordered additive regression p values per variant of genotyped-validation samples, after accounting for ancestry. SudRes2variant sc94:2,166,296 is a significant validated-variant (p < 2.3e-04, Bonferroni adjusted significance threshold shown by dotted line) and most other top-outliers (1.3e-04< p < 1.0e-03) are in SudRes1or SudRes2. (D) Allele and genotype counts for the most significant validated-variants from SudRes1 (c6:3,490) 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. Ancestry is the proportion of Population 2 ancestry, and loci SudRes1 and SudRes2are additive effects per each copy of the minor allele at the representative codominant markers c582:65,596 and sc94:2,174,117. ***p < 0.001.

Figure 3

Characterization of Biomphalaria sudanica SudRes1 genomic region. (A) Pooled-GWAS p values in SudRes1 regions (dashed boxes), which contain pooled-GWAS dual-variants (red) and other variants (blue), defined as in Figure 1. Contigs (grey 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, Fig. S2-S4). 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. S8A). 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 is shown at the top (in 100 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 is shown as in C, for 25 aa sligning windows in aligned proteins.

Figure 4

Characterization of Biomphalaria sudanica SudRes2 genomic region. (A) Pooled-GWAS p values in SudRes2 regions (dashed boxes), which contain pooled-GWAS dual-variants (red) and other variants (blue), defined as in Figure 1. Contigs (grey 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, Fig. S2-S4). 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 (69) comparing synteny of the SudRes2 region in Biomphalaria sudanica Bs111 genome (27) 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.25704and its paralog.