PTC2 region genotypes counteract Biomphalaria glabrata population differences between M-line and BS90 in resistance to infection by Schistosoma mansoni

Snail strains

The original M-line was reportedly created in the 1950s by crossing a pigmented Puerto Rican strain that was susceptible to Puerto Rican schistosomes, with an albino, resistant Brazilian strain, and then inbreeding with selection for albinism and susceptibility (Newton, 1955; Richards & Merritt, 1972). BS90 was reportedly isolated in Salvador, Bahia, Brazil, in the 1960s, and has been in labs in the US since the 1990s (Paraense & Correa, 1963, cited in Knight et al., 1999; C. Bayne, 2006, personal communication). Both have been maintained by the Biomedical Research Institute (BRI) schistosomiasis resource center (https://www.afbr-bri.org/schistosomiasis/) and in dozens of labs around the world for decades. We obtained BS90 in 2015 from Timothy Anderson’s lab at the Texas Biomedical Research Institute (originally from BRI), and M-line snails in 2018 from the BRI (their “T0” population). For our crosses we used the outbred M-line population and an inbred line of BS90 that we generated via two generations of selfing (line FSS5). Snails were housed in 7.5 liter plastic tanks containing artificial spring water (Ulmer, 1970).

Schistosome challenges

We challenged 6 to 9 mm diameter snails individually in 24-well plates as follows: 15 miracidia of PR1-strain S. mansoni were added to ~1 ml water in each well of a 24-well plate. Then a single snail was put into each well. Snails were left in the plates with miracidia for two hours, and were then moved to darkened tubs of ~24 snails per tub. Snails were kept in the darkened tubs for four weeks, and then were checked every week or two for shedding for six more weeks. Snails that did not shed within the 10 week window were classified as non-infected. Individuals that died before shedding were not included in the analysis. BS90 is always 100% resistant to PR1 (i.e., no BS90 snails ever got infected), while about 60% of our M-line snails get infected under these challenge conditions. In all experimental challenges for this study we always included at least 12 pure BS90 and 12 pure M-line snails as controls for schistosome infectivity. The Oregon State University Institutional Animal Care and Use Committee provided approval for this research under Animal Care and Use Protocols 2021-0213 and 5115.

PacBio sequencing and characterization of the PTC2 region in BS90 and M-line

We used an inbred line of BS90 (two generations of selfing) and an outbred population of M-line for our crosses. We obtained PacBio genomic assemblies from a snail from the inbred BS90 line and from an M-line snail. The BS90 and M-line samples were run on a Sequel II machine and assembled using Flye (Kolmogorov et al., 2019). We generated assemblies with high contiguity for both BS90 (2,534 fragments, length = 843.8 Mb, N50 = 2.3 Mb, largest scaffold = 15.5 Mb) and M-line (1,834 fragments, length = 852.5 Mb, N50 = 15.1 Mb, largest scaffold = 38.5 Mb). In both genomes, PTC2 occurs on assembled fragments over 4 Mb long, more than 1 Mb from either end, allowing for detailed assessment of the surrounding genomic region. We used alignments of the BS90 and M-line sequences to find regions of synteny around the PTC2 region in which we could place PCR markers for the marker-assisted backcrossing (Fig. 1). We used GENSCAN (Burge & Karlin, 1997), HMMER (Finn, Clements & Eddy, 2011), DeepTMHMM (Hallgren et al., 2022), Vectorbase (Giraldo-Calderón et al., 2015), and BLASTP to annotate coding regions in these assemblies, supplemented by our previous genomic and transcriptomic data from the 1316R population of B. glabrata (Tennessen et al., 2020). Assemblies have been made available through NCBI Genome, BioProject Accession PRJNA847915.

Backcross experiment

We used PCR to score five markers in and around the PTC2 region (Fig. 1). For ease of scoring, we designed primers to target gel-scorable indels that differed between our M-line and BS90 parental strains (Table S1). To drive the BS90 PTC2 region into an M-line genetic background we created F1s, and then backcrossed them to pure M-line (“B into M” backcrossing). We genotyped backcross offspring until we found recombinants inside one of the outer two markers, and then used those offspring in repeated rounds of backcrossing and genotyping until we found individuals that were homozygous M-line at the two outer markers and heterozygous at the inner three markers (i.e., genotype MHHHM at the five loci, where M = homozygous for M-line allele and H = heterozygous; Fig. 1). This took three or four rounds of backcrossing per tank. Individual snails with the desired double recombination event (MHHHM) were backcrossed one final time to further remove background BS90 genome and to produce a large number of siblings which all had the MHHHM genotype. We then let tanks of MHHHM siblings breed randomly and expand to create populations containing MMMMM, MHHHM and MBBBM individuals on a randomized genetic background for challenge testing. Each tank was started from an independent double recombinant event. We created and tested four independent “B into M” populations (i.e., started by MHHHM individuals), named M8.2, M8-2.2, M13.2 and M17. We challenged and genotyped N = 159, 89, 85 and 68 surviving snails from each of these populations, respectively.

We used the same genotyping and backcross approach to drive the M-line PTC2 region into a BS90 genetic background (“M into B” backcross population). We were able to find only a single BHHHB individual (where B = homozygous for the BS90 allele). We used that snail to create a single “M into B” population (population B7.2), following the same protocol as with the “B into M” populations. B7.2 produced only 11 infected individuals out of 116 challenged survivors (9.2%).

F2 Intercross population to test regional effect

We also crossed BS90 × M-line to create an F2 cross population that we allowed to breed randomly for almost two years while we conducted the backcrossing experiments (this population was created as a backup mapping population in case we rejected the PTC2 hypothesis). We pooled F1s from two individual BS90 × M-line individual snail crosses, and then used their offspring to found the F2 breeding population. Our lab was abruptly shut down during the covid pandemic, during which time we were able to maintain the F2 breeding population. However, we lost the tubs that contained the original founding parents, so we don’t have their DNA. The BS90 parents were all from an inbred line, so are likely to be similar to the source inbred line, which we retain and for which we have a PacBio assembly. As with the backcross populations, M-line parents were from our outbred population. So there could be genetic variation within the M-line contribution to the F2 population.

Seven months after setting up the F2 population, we started challenging ~two dozen snails every few months until we had accumulated >300 phenotyped individuals. We challenged a total of 424 individuals, of which 333 survived and 62 (19%) became infected. We successfully genotyped 329 of these 333 individuals at the five loci that span the backcross region. We chose 55 infected and 55 non-infected individuals as cases and controls, and additionally genotyped them via PCR at OPM-04 (primers in Table S1) in order to compare the effects of allelic variation at OPM-04 vs PTC2. The OPM-04 amplicons were scored via Sanger sequencing. In the full sample of 333 individuals we observed evidence of null alleles at marker loci 0, dn1 and dn2, while markers up1 and up2 showed Hardy-Weinberg Equilibrium. So we analyzed genotypes at marker up1 as a proxy for PTC2 genotype (only ~1 Mb away, so should be in high linkage disequilibrium with marker 0). We noticed no evidence for null alleles in the backcross populations, but that is expected because each generation we used individuals that were successfully typed at all five loci.

Statistical analysis

The various populations created, and samples used, for this study are summarized for clarity in Fig. S1. The effect of segregation for BS90 vs M-line alleles on infection status in the backcross and F2 populations was analyzed via binary logistic regression in Systat 13.2. For the four “B into M” backcross populations, we tested main effects of genotype and population, and the genotype × population interaction. We tested additive and completely dominant modes of inheritance. We observed too few infected snails from the “M into B” population to have any power for a statistical test, so we simply show the data from B7.2 and note that the general pattern is consistent with what we observed in the “B into M” populations.

We tested the effect of PTC2 genotypes using marker 0 in the backcross populations, and using marker up1 in the F2 population (to avoid the null allele at marker 0). For locus OPM-04, we observed two sequences from the M-line population segregating in the F2s, both of which were distinct from the BS90 sequence (Fig. S2). So we analyzed those data as bi-allelic, B vs M (where M = either of the two M-line alleles).

Effect of PTC2 in mostly M-line genetic backgrounds

The number of infected and not-infected individuals per PTC2 genotype in each “B into M” backcross population is shown in Fig. 2. The additive model showed only a slightly better fit than an M-allele-dominant model (AIC = 458.8 vs 459.6). There were strong main effects of both genotype and population on infection status (Tables S2, S3). The genotype-by-population interactions were not significant, despite the visual appearance of a different genotypic effect in population 17 in Fig. 2. Under an additive model, each additional BS90 allele roughly doubles the odds of infection (odds ratio ± standard error = 2.2 ± 0.78; P = 0.025). The estimated effect is similar under an M-allele dominant model (odds ratio = 3.8 ± 1.9 for the BB genotype relative to M_ = MM or BM; P = 0.009).

One unexpected result is that we see a significant heterozygote excess at the PTC2 in population M8.2, and a trend of heterozygote excess in the other three “B into M” populations (F_IS at locus 0 = −0.18, −0.12, −0.06 and −0.21 in populations M8.2, M8-2.2, M13.2 and M17, respectively). Combining all four populations together, observed heterozygotes (N = 228) significantly exceed expectation (χ² test, P = 0.003), whether expected heterozygosity is calculated from each population separately (expected N = 198.7) or from all samples pooled (expected N = 199.7). We see similar negative F_IS at all three of the segregating backcross-region loci (up1, 0 and dn 1) in these populations, so the negative F_IS seems unlikely to result from technical artifacts.

Effect of PTC2 in a mostly BS90 genetic background

We see infection in all three genotypes, and the trend is consistent with what we observed in the four “B into M” populations (fractions infected were 3/33 = 0.11 BB, 7/64 = 0.11 BM, and 1/19 = 0.06 MM). Population B7.2 also showed a trend of heterozygote excess at the backcross region loci (F_IS = −0.12 at all three loci), although the deviation from Hardy Weinberg Equilibrium is not significant.

Effect of PTC2 in the F2-cross genetic background, and comparison with OPM-04

We examined effects of the PTC2 and OPM-04 regions in the F2 population in order to test whether the locus linked to OPM-04 in Knight et al. (1999) was still nearby, even if not within the region we introgressed into our backcross lines. As in the backcross populations, the BS90 allele conferred significantly enhanced susceptibility in the full sample from the F2 population (N = 329 successfully genotyped) (Tables S2, S3, Fig. 3). The additive and M-allele-dominant models were equally supported. Odds ratios for increased susceptibility were 1.5 (± 0.32) per B allele under an additive model (P = 0.041) and 1.8 (± 0.54) for BB relative to M_ under an M-allele dominant model (P = 0.044).

The subsample of 55 cases and 55 controls also showed significantly enhanced susceptibility caused by the PTC2 B allele. Odds ratios = 2.0 (± 0.61; P = 0.029) and 3.3 (± 1.49; P = 0.007) for the genotype effect under additive and M-allele dominant models, respectively (Tables S2, S3). There was a similar pattern at OPM-04 (Fig. 4), although the genotype effect was not significant (P = 0.286 and 0.380 for additive and M-allele dominant models, respectively; Tables S2, S3). So unlike our original expectations from Knight et al. (1999), OPM-04 does not appear to be anywhere near a locus at which the BS90 allele confers dominant resistance. Instead, it shows a pattern similar to that of the closely-linked PTC2 region.

As mentioned previously, there were obvious null alleles at loci 0 to dn2. But loci up2 and up1 were in Hardy-Weinberg Equilibrium (F_IS = −0.03), so we do not see the pattern of heterozygote excess that was apparent in the backcross populations.

PacBio assemblies

In the 1316R population there are three PTC2-region haplotypes that are ordered from most to least protective as R > S2 > S1 (Tennessen et al., 2020). Our M-line haplotype is over 99% similar to the previously described S2 haplotype, with a total of only six nonsynonymous differences between M-line and S2 across all 11 genes in the region. Our BS90 haplotype is distinct from other known haplotypes; a portion of it is similar to S2/M-line, while the remainder is closer to S1 and is extremely divergent from S2/M-line (Figs. 5, 6). At the amino acid level, PTC2 genes 1–5 are 95% identical between BS90 and S2/M-line, while the remaining genes are either only present in BS90 (gene 7), only present in M-line (genes 6 and 8), or in the case of gene 9 (BGLB018951) present in both but over 50% divergent in sequence between them. Thus, this divergent portion of PTC2 seems most likely to explain the observed phenotypic difference between BS90 and M-line.

Recall, however, that the ~300 kb PTC2 region is embedded within a larger region (~850 kb in BS90 and ~1,300 kb in M-line) that shows almost no sequence similarity between our M-line and BS90 assemblies (Fig. 6). Thus, it is possible that the causal gene/s lie in the flanking regions, which we also introgressed. Both haplotypes possess unique genes in this zone, although as with many snail genes, little is known about their function and most appear to be transposable elements. Notably, three AIG genes occur in this region in M-line but not BS90 (in possession of an AIG1 domain; homologs of each other and of BGLB016854, BGLB012732, and BGLB004870). AIG genes are suspected to have immune function in Biomphalaria (Lu et al., 2020). We have confirmed with RNA-Seq data from haplotype S2 (Accession SRX8534561) that all three AIG1 genes are expressed. Farther out from PTC2, our BS90 and M-line assemblies are again similar in sequence and so seem less likely to contain the causal genes, though this remains a possibility. For example, a toll-like receptor gene BGLB025483 occurs near locus dn1 in both populations and has 24 amino acid differences between the haplotypes. Pila et al. (2016b) showed that knocking down expression of a different Toll-like receptor makes BS90 more susceptible, so this is an interesting class of proteins.

One interesting observation from our BS90 and M-line PacBio assemblies is that the 1.2 kb OPM-04 amplicon sequence published by Knight et al. (1999) (genbank AF078109) is almost the same in our BS90 assembly, but is interrupted by a ~1.7 kb insertion in our M-line assembly (Fig. S2). So if Knight et al. (1999) M-line and BS90 also differed at this indel, then that might explain why under a low-stringency RAPD-PCR reaction they got amplification from BS90 but not from M-line chromosomes, resulting in a dominant marker for BS90.

Why were we unable to replicate Knight et al. (1999) result?

We were surprised to see no effect of the OPM-04 region in our BS90 × M-line F2 population. Possible explanations for this failure to replicate Knight et al. (1999) result include: (1) the strains of “BS90”, “M-line” or schistosome we used differ from theirs, and (2) we were not measuring the same phenotype.

(1) BS90 and M-line have been in culture in many labs for decades, so there is no guarantee that our lines are genetically the same as the those used by Knight et al. (1999) over 20 years ago. BS90 is phenotypically consistent in that it remains 0% susceptible to PR1 schistosomes. On the other hand, Knight et al. (1999) reported that their M-line was 100% susceptible to PR1, whereas under our experimental conditions M-line was only ~60% susceptible. Thus, our outbred M-line population may contain variation for resistance. Perhaps their M-line was more isogenic and carried a susceptibility allele near OPM-04 that is absent or in low frequency in our M-line. Another possibility is that effects of genetic background change the penetrance of whatever allele was linked to OPM-04 in their M-line. It is also possible that the SmPR1 strain of schistosome we used differs in some way from the strains they used.

(2) Another plausible explanation for the difference in results between the two studies is that we were not measuring the same phenotype. In Knight et al. (1999) study they report phenotype ratios in BS90 × M-line F1s and F2s that are consistent with single-locus Mendelian inheritance, and in which the BS90 allele is dominant and confers complete resistance (i.e., 0% infected F1s, ~25% infected F2s). Older literature on crosses among strains that differed in resistance to PR1 reported data consistent with single-locus inheritance when using mature snails, but data more consistent with polygenic inheritance when using juvenile snails (e.g., Richards, 1984). Thus, it is possible that our experimental conditions (e.g., snail size or maturity, number of miracidia per challenge, or some other environmental factor) reduced the importance of whatever highly penetrant genotype was linked to OPM-04 in Knight et al. (1999) study. For example, although both studies used similar-sized snails, Knight et al. (1999) challenged each snail with 100 miracidia, whereas we challenged with only 15 miracidia (because higher numbers caused excess mortality in our hands).

There is no obvious reason to doubt the data in Knight et al. (1999). Plus, it seems a remarkable coincidence that in a 0.9 Gb genome of 18 chromosomes, such an important region as PTC2 would lie only 5 Mb, and almost exactly the expected number of cM from OPM-04, yet not contain the causal gene/s (Tennessen et al., 2020). Therefore, we suggest that if one could reproduce the highly-penetrant, dominant BS90 phenotype observed by Knight et al. (1999) and previous workers, then it would be worth re-testing whether the OPM-04 region shows an association, and whether the causal gene/s linked to that marker lie in the PTC2 region.

Genetic background effects

The sizable population effects among the four “B into M” populations were stronger than the genotype effect (Fig. 2; Table S3), and presumably reflect their different genetic backgrounds. This shows that other loci have a very large effect on resistance in this system. Each population was started using different M-line individuals from the outbred M-line, and each must retain different residual segments of donor line (BS90) genome that is still segregating. So there is plenty of opportunity for these populations to differ in genetic background, and we don’t know to what extent that results from variation in the base M-line population or from different retained segments of BS90 chromosomes. The population-by-genotype interactions were not significant, although the trend towards a weaker genotypic effect in population M17 was intriguing (Fig. 2). The effect of SOD1-region alleles on resistance appears to depend on genetic background (Bonner et al., 2012), so similar epistatic effects would not be surprising.

Negative F_IS in the PTC2 region

It was intriguing that F_IS in the PTC2 region was negative in all of the backcross populations, and significantly so in one of them, and in the set of all populations combined. The mortality rate after parasite challenge in these populations ranged from 6% to 31%, and if these deaths were primarily BB and/or MM homozygotes, that could entirely explain the observed heterozygote excess. However, observed mortality was not correlated with F_IS (for example, population 13.2 had the highest mortality but the lowest heterozygote excess) and heterozygotes did not show any overdominant effect on odds of infection (Figs. 2, 3). So if there is overdominance for survival, it is not obviously owing to the schistosome challenges. Unfortunately we did not keep the backcross populations at the end of the experiment, so we cannot go back and test for Hardy-Weinberg equilibrium in the unchallenged populations. A slightly negative F_IS can also result when a small effective number of breeders (N_b) contributes to a sample, but to get F_IS in the observed range of −0.15 to 0.20 would require N_b ~2–3 individuals (Waples, 2015). Such a small N_b is not out of the question, but seems low. That there might be overdominance for survival in a region of such obvious immune relevance is a hypothesis worthy of further exploration, as this would be one mechanism that could maintain genetic diversity in the region.

Additional insights into the PTC2 region

Allelic variation in the PTC2 region has the largest known effect on resistance identified to date (Tennessen et al., 2020). The extreme structural and non-synonymous sequence variation in this region, and its enrichment for membrane-spanning protein genes, strongly suggests host-parasite interacting molecules under selection. Of course, schistosomes are an introduced species in the New World, so this region no doubt evolved with other pathogens. It will be important to determine which loci are involved in interactions with schistosomes, and if multiple loci are involved, how they interact. The 2- to 4-fold effect of BB vs MM observed here was weaker than previously reported for PTC2 (Tennessen et al., 2020), but comparable to the effect of other important B. glabrata resistance loci like sod1 or RADres (Tennessen et al., 2015a). Such effects are real and likely to be biologically meaningful, but they are easily missed. Here we benefitted from the use of hundreds of snails. A power analysis shows that if we had used half the sample size in this study, we likely would not have inferred a significant genotype effect. Thus lower-powered studies might miss these more subtle genotype-phenotype correlations.

Although our addition of two new assemblies of the region revealed new structural variants and additional sequence variation, the two new haplotypes closely resemble the other known haplotypes or recombinants of them. Comparing the two new haplotypes to the three previously described haplotypes allowed us to make some a priori predictions about which loci are most likely to be involved. For example, the M-line haplotype closely resembles S2 while the BS90 haplotype is closer to S1, and thus the 7-fold decrease in infection odds conferred by S2S2 relative to S1S1 (Tennessen et al., 2020) could have the same genetic basis as the effect reported here, perhaps due to one of genes 6 through 9 or an AIG gene. We are currently attempting functional knockdowns of these genes to provide definitive tests.