Baylisascaris laevis, Leidy, 1856

Barrera, McIntyre A., Janes, Jasmine K. & Gorrell, Jamieson C., 2022, Molecular phylogenetics and systematics of two enteric helminth parasites (Baylisascaris laevis and Diandrya vancouverensis) in the Vancouver Island marmot (Marmota vancouverensis), International Journal for Parasitology: Parasites and Wildlife 19, pp. 301-310 : 305-307

4.1. Systematics of Baylisascaris laevis

Camp et al. (2018) proposed two clades in Baylisascaris based on reconstructions across eight genes for seven of the 11 recognized species in the genus. Our reconstructions with additional samples and species supported these two clades; however, our analyses generally placed B. tasmaniensis as sister to all other Baylisascaris species, except in the cox 1 and cox 2 trees (Figs. S6, S7, and S9), the BI 28S tree (similar to the BI 28S tree from Camp et al. (2018; Fig. S11)), and the MP and BI ard 1 trees (Figs. S12–S13). The proportion of unambiguous sites selected by Gblocks in our Baylisascaris alignments (Supplementary Table S3) were comparable to the proportion of characters retained by Camp et al. (2018; Supplementary Data) in datasets filtered based on a 60% posterior probability threshold with ProAlign (L¨oytynoja and Milinkovitch, 2003).

The topology of ingroup species (excluding B. laevis ) in our BI mt-tree matched Camp et al. (2018), although B. procyonis from California and Connecticut were placed in a polytomy that was sister to B. columnaris ( Fig. 2 ). Similarly, our BI nr-tree showed an alternate topology for Baylisascaris spp. in the clade with ursid and ailurid hosts and placed B. tasmaniensis as sister to all other Baylisascaris ( Fig. 3 ). Our BI combined tree followed our mt-tree topology for the ursid and ailurid clade, but otherwise reflected our nr-tree arrangement of Baylisascaris spp. ( Fig. 4 ). These differences may reflect the fact that our analyses excluded two nuclear genes included in the analysis by Camp et al. (2018).

Our analyses placed B. laevis in the clade with (but also sister to) species that parasitize skunks, raccoons, and gulonine mustelids ( B. columnaris , B. procyonis , and B. devosi , respectively). This relationship confirms the morphological grouping assigned by Sprent (1968)

305

and constitutes the first molecular phylogenetic placement of B. laevis . It has previously been suggested that the direct life cycle of B. laevis is the result of a capture event in its evolution, during which an ancestral form with a carnivoran host was able to develop to adulthood within its rodent intermediate host ( Anderson, 2000). The position of B. laevis as sister to one clade, but not the other, suggests that this change likely occurred after the common ancestors of each clade diverged.

The pairwise sequence divergence between the VI and mainland samples of B. laevis ranged from 0.11 to 2.86% ( Table 4), with the greatest divergence observed between VI and Idaho at ard 1. However, variation between VI and Idaho was surpassed by variation between Idaho and Alaska at the same locus (3.39%). Locus ard 1 is an EPIC marker selected to capture variability between closely related species ( Camp et al., 2018); it is highly variable by design, and thus likely not representative of mean divergence for the purpose of species delimitation. ITS and cox 1, on the other hand, have been used to estimate interspecific divergence and prospect for cryptic species within nematode genera ( Powers et al., 1997; Blouin, 2002).

P´erez Mata et al. (2016) used ITS to classify the novel species Baylisascaris venezuelensis and reported interspecific pairwise divergences of 2.8–10.0% among B. venezuelensis , B. schroederi , B. transfuga , and B. procyonis . Clearly, the extent of interspecific divergence within Baylisascaris encompasses a range of values depending on both the genes and species pairs examined. Baylisascaris columnaris and B. procyonis are considered distinct species with shallow genetic divergence ( Franssen et al., 2013; Camp et al., 2018). In our analyses, these species formed a sister clade to B. laevis , and are thus a reasonable comparison for interspecific divergence values in closely related species.

The mean pairwise divergences between VI and mainland B. laevis at ITS and cox 1 (2.08 and 0.72%, respectively) were comparable to (or greater than) the divergences we calculated between B. columnaris and B. procyonis for the same genes in our alignment (1.12 and 0.88%, respectively). However, we note that Camp et al. (2018, Supplementary Data) reported mean divergences of 0.32 and 0.96% at ITS and cox 1, respectively, between B. columnaris and B. procyonis in their alignments. These discrepancies may be due to differences in either the length of alignments (e.g., ITS sequences from Camp et al. (2018) were 889–975 bp long, whereas our alignment was truncated to 756 bp), or the

306

algorithms used to construct them. Furthermore, Carlson et al. (2021) found B. procyonis had a maximum divergence of 1.6% at cox 1 when comparing populations on opposite sides of the continental divide.

The divergence of B. laevis between VI and the mainland may be within the range of interspecific divergences for Baylisascaris nematodes. However, our current data are insufficient to fully test our hypothesis of speciation for B. laevis on VI. More sophisticated analyses of species delimitation may be able to resolve this uncertainty. For example, Camp et al. (2018) included ITS and cox 1 in a species delimitation analysis using BP&P v3.3 (Bayesian Phylogenetics and Phylogeography; Yang, 2015). Further analyses must be informed by more extensive population sampling of B. laevis from both M. vancouverensis and mainland marmot populations (e.g., M. caligata ), which we were unable to obtain here, in order to resolve the phylogeography and population genetic structure of B. laevis across its entire geographic and host ranges (P´erez-Ponce de Leon´and Nadler, 2010). Further molecular analyses should ideally be complemented with detailed morphological comparisons of B. laevis from an equally broad range of hosts and localities in order to characterize its intraspecific variation within and among populations.

Figures

Fig. 2. Bayesian consensus phylogram based on Baylisascaris and outgroup alignments of concatenated mitochondrial sequences (12S, cox1, and cox2) from GenBank and this study (B. laevis). Branch labels represent Bayesian posterior probabilities. Branch lengths are scaled to expected number of substitutions per site. Abbreviations refer to sampling sites (AK = Alaska, ID = Idaho, CT = Connecticut; IL = Illinois; CA = California; WV = West Virginia; ALB = Alberta). See also Table 1.

Fig. 3. Bayesian consensus phylogram based on Baylisascaris and outgroup alignments of concatenated nuclear sequences (28S, ITS, and ard1) from GenBank and this study (B. laevis). Branch labels represent Bayesian posterior probabilities. Branch lengths are scaled to expected number of substitutions per site. Abbreviations refer to sampling sites (AK = Alaska, ID = Idaho, CT = Connecticut; IL = Illinois; CA = California; WV = West Virginia; ALB = Alberta). See also Table 1.

Fig. 4. Bayesian consensus phylogram based on Baylisascaris and outgroup alignments of concatenated mitochondrial (12S, cox1, and cox2) and nuclear sequences (28S, ITS, and ard1) from GenBank and this study (B. laevis). Branch labels represent Bayesian posterior probabilities. Branch lengths are scaled to expected number of substitutions per site. Abbreviations refer to sampling sites (AK = Alaska, ID = Idaho, CT = Connecticut; IL = Illinois; CA = California; WV = West Virginia; ALB = Alberta). See also Table 1.