Population genetics and dispersal of the flatworm, Polycelis coronata: a test of the habitat stability hypothesis

ABSTRACT The habitat stability hypothesis states that species in spring-like habitats have little reason to disperse compared to species in temporary habitats. Planarians commonly inhabit springs around the world and they have long been considered poor dispersers. Recently, however, genetic analyses have shown contradictory results on the dispersal of planarians. Asexual planarians that can establish a new population by colonization of a single individual showed little genetic differentiation between sites separated by hundreds of kilometers, whereas species inhabiting springs showed deep differentiation between sites separated by hundreds of meters. The latter results are consistent with the habitat stability hypothesis. We used the cytochrome oxidase subunit I gene from 468 individuals of Polycelis coronata, an asexual species, collected from 50 sites, nested in 26 tributaries, in 4 catchments of the Wasatch Mountains of Utah, USA, to explore the dispersal capabilities of P. coronata. The longest distance between sites was 66 km. Despite this small spatial extent, we found that 77% of the 130 haplotypes were collected from a single site and 89% from a single catchment. FST values between local populations in the same tributary (0.221, 0.266, 0.389) were similar to the average FST values in different catchments for other headwater taxa. Also, variation among individuals accounted for the majority of genetic structuring with little differentiation beyond the scale of a single site. Dispersal is very slow in this species which is consistent with the habitat stability hypothesis. However, we suggest that other explanations also warrant consideration. We also identified two potential cryptic species suggesting a high degree of hidden variation at the level of species in this genus.


Introduction
Planariidae possess several traits that may limit their ability to disperse within stream ecosystems. They are obligate aquatic species without a terrestrial phase in their life cycle, which eliminates active, overland dispersal. Also, they may rarely encounter passive vectors (e.g. birds) because the adults are nocturnal and cryptic, and if they produce resting eggs, their cases are firmly attached to the undersides of rocks (Walter 1907;Kolasa et al. 2009). Plus, macro-turbellarians, like Polycelis coronata, do not swim (Kolasa et al. 2009) and they are rarely found in the drift (Minshall &Winger 1968). Consequently, biologists have long thought that planarians have a poor dispersal capability (e.g. Ball 1969Ball , 1974. However, we still have much to learn about dispersal in flatworms because CONTACT Russell B. Rader russell_rader@byu.edu

Sampling design and study system
We sampled 26 headwater segments nested in 4 adjacent catchments (Provo, American Fork, Little Cottonwood, and Weber) to determine the distribution of haplotypes of P. coronata in the central portion of the Wasatch Mountains ( Figure 1). All segments were permanent, and included 3 small third-order streams, 11 second-order tributaries, and 12 first-order tributaries. All segments were sampled at an upstream and downstream site typically separated by »500-1000 m except for two segments that contained a single site. Thus, we sampled a total of 50 sites in this study (Appendix 1). A site consisted of a stream reach »50 m in length with both riffles and pools.
We defined the local population as the planarians collected at a site. The Weber (6070 km 2 ) and Provo (1761 km 2 ) are two of the largest catchments in the central Wasatch Mountains, whereas Little Cottonwood (119 km 2 ) and American Fork (107 km 2 ) are more typical of smaller drainages (Shiozawa & Rader 2005). Each catchment eventually drains into The Great Salt Lake, and all catchments are separated by steep mountainous terrain. In addition to sampling tributaries, we conducted a thorough search for P. coronata in mainstem segments in each catchment, except for the Weber River where private property forced limited access. Also, P. coronata and several aquatic insect species were absent in the central portion of the Provo drainage because of frequent dewatering by agriculture ( Figure 1).
During July and August of 2008, we collected 10 P. coronata by hand from the undersides of submerged rocks at each site. We also collected 10 individuals of P. coronata from five distant locations to determine how many haplotypes, if any, had a broad distribution within the region (Logan River, UT; Salt River, WY; Colorado River, CO; Snake River, ID; Rio Grande River, NM). We could show the ability of P. coronata to disperse over long distances, if some haplotypes are distributed between the central Wasatch Mountains and these distant locations.  (Provo River, Weber River, American Fork River, and Little Cottonwood Canyon). Closed circles represent sites where Polycelis coronata occurred; open circles show sites where this species was absent. Dashed lines indicate the boundaries of each catchment. In the lower right corner is the state of Utah showing Salt Lake City (star), Utah Lake, and the extent of the Wasatch Mountains (shaded area).
We collected a timed (1 minute) benthic sample (D-frame sweep net with a 250-mm mesh) at each site to estimate the average relative abundance of P. coronata in the central Wasatch Mountains. We also made visual estimates of substrate stability (presence of bryophytes and embeddedness) and collected chemical data (pH, conductivity, and water temperature) as spot measurements during mid-day at each site. Stable sites would be spring-like and fed by a high proportion of groundwater. Thus, stable sites would show high conductivity, low pH, cool water temperatures, embedded substrate, and an abundance of bryophytes. Flow records were not available for any of these small tributaries. However, bryophytes are slow-growing plants indicative of rocks that have seldom been rolled by high flows and are thus a reliable indicator of stable substrate conditions (e.g. Bowden 1999;Suren & Duncan 1999;Milner et al. 2000).
Final concentrations for the polymerase chain reaction (PCR) per 25 mL were 25 ng template DNA, 0.25 mM of each primer, 0.625 units of Taq DNA polymerase, 0.1 mM of each dNTP, 2.5 mL of 10X reaction buffer, and 2.5 mM MgCl 2 . Amplification proceeded at 94 C for 2 min followed by 35 cycles of 94 C for 30 s, 48 C for 30 s, and 72 C for 75 s. Following these cycles, final elongation was held at 72 C for 7 min. We examined PCR products on a 1% agarose gel using SYBR safe DNA gel stain (Invitrogen, Eugene, OR, USA). We purified PCR products using a Montage PCR 96 plate (Millipore Billerica, MA, USA). Sequences were obtained via cycle sequencing with Big Dye 3.0 dye terminator ready reaction kits using 1/16th reaction size (Applied Biosystems, Foster City, CA, USA). Sequencing reactions were run with an annealing temperature of 52 C following the ABI manufacturer's protocol. We purified sequenced products using sephadex columns, and sequences were obtained using an Applied Biosystems 3730 XL automated sequencer at the Brigham Young University DNA Sequencing Center. All haplotypes were deposited in GenBank, accession numbers KY057039-KY057186.
Analysis of mtDNA sequence data DNA sequences were edited using Chromas Lite 2.0 (Technelysium, Tewantin, Queensland, Australia) and imported into BioEdit 7.0.5.2 (Hall 1999), then aligned by eye. An alignment algorithm was not necessary because there were no gaps in the sequences. Sequences were checked for unexpected frame shift errors or stop codons in Mega 6.06 (Tamura et al. 2013). Individual sequences were then reduced to 148 haplotypes using DnaSP 5.10.01 (Librado & Rozas 2009). We used RAxML 8.0.24 (Stamatakis 2014) for our maximum-likelihood (ML) analysis to examine the evolutionary relationships of our sequence data on the CIPRES cluster at the San Diego Supercomputer Center (Miller et al. 2010). We used the GTRGAMMA model for 1000 bootstrap replicates and for the final tree calculation. Due to a lack of more appropriate outgroups, trees were rooted with the European species, P. felina and P. tenius (GenBank accession numbers DQ666049 and AF178321, respectively). Mean p-distances for DNA were calculated using Mega based on the individual haplotypes.
All individuals from outside of our four catchments and those identified as potential cryptic species were eliminated from our analyses of P. coronata in the Wasatch Mountains. We used an analysis of molecular variance (AMOVA) implemented in the Vegan package of 'R' 2.8.1 (R Development Core Team 2010) with significance assessed using 10,000 random permutations to examine spatial patterns of genetic variation (Excoffier & Lischer 2010). We partitioned genetic variation into the following hierarchically nested groups: 'F IL ' was the variation among individuals (I) within a local population (L), while F ST was the variation among local populations. F ST was partitioned into F LS , F SC , and F CT . 'F LS ' was the variation among local populations (L) nested within segments (S), 'F SC ' was the variation among segments (S) within catchments (C), and 'F CT ' was the variation among catchments (CT). We do not expect significant differentiation among local populations at any level of the spatial hierarchy if this species has been recently dispersed by humans.
We used a Mantel test in 'R' (R Development Core Team 2010) to test the hypothesis that pairwise genetic distances were either correlated with straight-line geographic distances between local populations, or distances following the stream network measured by ArcGIS 10 (Environmental Systems Research Institute, Redlands, CA). A positive correlation between F ST and geographic distances would suggest the potential for passive overland dispersal by terrestrial vectors (e.g. humans) with a greater tendency to disperse between near sites than between far sites. However, a positive correlation between F ST and stream distances would indicate the potential for dispersal through the stream network.

Environmental conditions
Polycelis coronata occurred at 86% of the headwater sites within our study system. When present, it had an average relative abundance of 21% ( §11%) of the total invertebrates collected at a site. Along with chironomid midges, Baetidae, Simuliidae, and Neothremma alicia (Trichoptera, Uenoidae), it was consistently one of the five most abundant taxa in this area. We did not find P. coronata from mainstem segments where the average water temperature, pH, and conductivity in July and August was 12.8 C (SE § 0.68 C), 8.1 (SE § 0.08), and 185 mS/cm (SE § 11.4), respectively. This is compared to 10.4 C (SE § 0.68 C), 7.6 (SE § 0.06) and 311.6 mS/cm (SE § 17.6) at our headwater sites. Embeddedness of cobbles and boulders was about 30% at headwater sites and about 10%, if at all, in mainstem segments, except for the Provo River. Embeddedness was greater in the Provo River than in other mainstem segments because of the effects of river regulation (armoring). Also, patches of bryophytes were common and extensive at headwater sites and nearly absent in mainstem segments. These data indicate that our headwater sites were cool, stable, spring-like habitats with a high influx of groundwater.

Cryptic species
Editing resulted in a 763-bp fragment of COI from 541 individuals, which collapsed into 148 haplotypes. ML recovered one tree with a ¡ln score of ¡5046.65 (Figure 2). We identified two potential cryptic species shown in the upper left of Figure 2. The first branching in-group lineage, herein called Polycelis sp. A, and the sister lineage to P. coronata, herein referred to as Polycelis sp. B Figure 2. Maximum likelihood phylogram of all 148 COI haplotypes generated in this study. In the upper left corner are two European species used to root the tree and two potential cryptic species. The remaining tree is a blow-up showing relatedness of 130 P. coronata haplotypes. The scales for the branching lengths differ between the cryptic species (0.5) and the blow-up for P. coronata (0.01). This tree is split for ease of display as indicated by the dotted line. Bootstrap values were obtained from 1000 replicates, only values above 50 are shown. Each haplotype is followed by the number of individuals from the catchment(s) it occurred in: Provo (P), Weber (W), American Fork (AF), and Little Cottonwood (LC) or, from the state it was collected from for haplotypes of distant populations: Wyoming (WY), Colorado (CO), Idaho (ID), New Mexico (NW), and Logan (LO), Utah.
( Figure 2), were both strongly divergent from P. coronata. That is, the average p-distance between them and P. coronata was 14.4% as compared to average within species divergences for Polycelis sp. A, Polycelis sp. B, and P. coronata equal to 2.6% ( Table 1). The difference between P. coronata, Polycelis sp. A, and Polycelis sp. B was nearly as great as the difference between P. coronata and the two European species of Polycelis used to root the tree (Table 1). Thirty-five of the 36 individuals of sp. A were collected from within our study system (Figure 2). Of the 22 haplotypes collected from outside our study system, 5 were sp. B, whereas the remaining 17 haplotypes were P. coronata.

Haplotype distribution and relatedness in P. coronata
Of the 763 base pairs, 622 characters were constant, 121 characters were parsimony informative, and 20 were parsimony uninformative. Eight of the 15 haplotypes of P. coronata from distant sites were not collected from any of the four catchments (Haps 137,138,142,143,144,146,147,148) showing that there is considerable genetic diversity within this species outside our study system. We deleted haplotypes from distant sites and from Polycelis sp. A and sp. B, and used the remaining 130 haplotypes from 468 individuals of P. coronata to test our first hypothesis of no genetic differentiation across our study system. Relationships within P. coronata were poorly resolved with little bootstrap support ( Figure 2). Most haplotypes were quite closely related being separated from others by only one or two base pairs. Although most haplotypes were confined to a single catchment, several were closely related to other haplotypes in different basins, especially between the Provo and Weber catchments (lower right corner) and the American Fork and Little Cottonwood catchments (upper right corner).
Given the short distances between sites in this study (66 km), a well-connected species would show little differentiation even between catchments but especially between sites within segments. Contrary to the first hypothesis, the AMOVA showed the opposite result. The difference between local populations within segments (F LS ) was significant suggesting that there was little gene flow between populations in the same stream segment (Table 2). Also, differences among individuals within a local population explained the greatest amount of the total genetic variation (F IL in Table 2). Of the 10 individuals collected at a site, the number of haplotypes ranged from 1 to 9 with an average of 5. That is, 50% of the individuals at a site had different haplotypes. Note: n/a is not applicable. Table 2. Hierarchical analysis of molecular variance of P. coronata. 'F IL ' was the variation among individuals (I) within a local population (L) and F ST , the variation among local populations was partitioned into F LS , F SC , and F CT : 'F LS ' was the variation among local populations (L) nested within segments (S); 'F SC ' was the variation among segments (S) within catchments (C); 'F CT ' was the variation among catchments (CT There was a positive relationship between haplotype abundance and haplotype distribution. The majority of P. coronata haplotypes were rare in our study system. Seventy-one percent (92 haplotypes) were represented by a single individual (Figure 3), whereas 12 haplotypes were considered abundant (10 individuals). The abundant haplotypes were distributed across the greatest number of sites ( Figure 3). However, most of the abundant haplotypes were confined to one or rarely two catchments, and only one was collected from all four drainages (Hap 26). Thus, the majority of haplotypes were very rare and had a very restricted distribution. That is, 77% of the haplotypes of P. coronata in our study system were collected from a single site and 89% occurred within a single catchment (Figure 2), which is also inconsistent with our first hypothesis that all haplotypes could disperse to all sites.
Of the 15 P. coronata haplotypes from distant sites (Figure 2), 7 from Wyoming, Idaho, and Colorado were also collected from our study system (Haps 7,22,41,55,116,117,140). This confirms our second hypothesis and shows the potential for long-range dispersal in this species. However, we found no relationship between genetic distances (F ST values) and distances following the stream channel (R 2 D 0.16; p D 0.24) suggesting that P. coronata typically does not disperse through the stream network. Instead, we found a weak positive correlation between genetic distances and the straight-line distance between sites (R 2 D 0.43; p D 0.001) indicating a slight degree of spatial autocorrelation. That is, dispersal appears slightly more common among near sites than distant sites.

Discussion
Polycelis coronata in the central Wasatch Mountains is a poorly connected species. That is, 77% of the haplotypes were confined to a single site and 89% were restricted to a single catchment. Also, average F ST values for P. coronata between local populations in the same tributary (0.221, 0.266, 0.389) were similar to the average F ST values in different catchments (0.17, 0.318, 0.38, 0.39, 0.57) for other headwater taxa (Blephariceridae, Simuliidae, shrimp) that also poorly disperse (Hughes et al. 1995(Hughes et al. , 1996Wishart & Hughes 2001;Finn & Adler 2006, Finn & Adler 2006. These data show that dispersal is rare in P. coronata even between sites in the same tributary. Rare dispersal is consistent with the habitat stability hypothesis, which states that species in permanent habitats are not likely to disperse compared to species from temporary habitats (e.g . Roff 1990;L azaro et al. 2009). All of the sites in this study were permanent but stability can also be defined by flow characteristics which create sheer stress and dislodged substrate (e.g. M erigoux & Dol edec 2004, Biggs et al. 2005). Our sites undoubtedly showed a gradient from very stable seasonal flows to quasi-stable conditions where smaller substrate particles may be dislodged during spring run-off. However, these sites were very different from mainstem segments and the presence of moss is indicative that most particles were rarely disturbed. Plus, most studies examining the habitat stability hypothesis show a similar gradient in flow conditions within stable sites as in this study (e.g. Baggiano et al. 2011).
Although aquatic species in temporary habitats may frequently disperse to avoid desiccation (e.g. Bilton et al. 2001), it is risky to assume that species rarely disperse from stable habitats. A variety of mechanisms might promote dispersal in stable habitats (e.g. Gandon & Michalakis 1999). For example, high population densities may result in competition causing individuals to disperse to find more abundant resources (e.g. Baguette et al. 2013). As an alternative explanation to the stable habitat hypothesis, low rates of dispersal in P. coronata may be related to habitat specialization. Species should have a low tendency to disperse if there is a low probability of finding a particular type of habitat (e.g. Wishart & Hughes 2001). Narrow thermal tolerance, as in P. coronata, may select against dispersal if there is a high probability of landing in inhospitable habitats (Dobzhansky 1950;Vannote & Sweeney 1980) such as, warm downstream reaches. Species that are habitat specialists and transported by vectors may frequently land in unsuitable habitats, which may drive their evolution away from a dispersive phenotype (e.g. Jacob et al. 2015).
Vectors may be the only means of dispersal in P. coronata because it is an obligate aquatic species, it does not swim, and it is rarely found in the drift. Plus, it has a narrow temperature tolerance which probably prevents it from dispersing though warmer downstream reaches. Although it may have dispersed through the stream network during glacial periods, we found no correlation between genetic distances and the distance between sites following the stream channel. Instead, we found a weak correlation between genetic distance and the straight-line, overland distance between sites. Dispersal by terrestrial vectors (e.g. birds and mammals) is most consistent with this pattern. If vectors are the primary means of dispersal, then we should not be surprised that there was a weak relationship between genetic distances and geographic distances because individuals of P. coronata would be transported various distances to a variety of habitats by different vectors thus, introducing a lot of variation.
Several studies have shown that organisms move non-randomly through the environment to settle in habitats that may enhance their individual performance (e.g. Edelaar et al. 2008, Dreiss et al. 2012. This is called the habitat matching hypothesis (Edelaar et al. 2008). Passive dispersal on vectors should be more stochastic than habitat matching because it depends on the behavior and habitat selection of other species. Planarians transported by vectors cannot control where they land (e.g. Hebert & Payne 1985). Thus, contrary to the habitat matching hypothesis, the distribution of haplotypes and patterns of gene flow should be more random where passive dispersal on vectors is the primary means of dispersal.
Why do individuals account for the majority of genetic structuring with little differentiation at larger scales? Stochastic dispersal by vectors combined with little connectivity between sites would reduce genetic structuring at larger scales. That is, the initial stochastic distribution of P. coronata, possibly attributed to dispersal by vectors would be maintained if there was little dispersal among segments and catchments to homogenize the distribution of haplotypes at larger scales. For example, if haplotypes could disperse within a catchment but not between catchments, then we would expect significant genetic structuring at the scale of catchments. But, if haplotypes rarely disperse beyond a site, then we would expect individuals within sites to account for the majority of genetic structuring, which is what we found in this study. Thus, perhaps we should not expect to find genetic structuring at larger scales when species rarely disperse by vectors.
We found a positive relationship between the distribution and abundance of P. coronata haplotypes, plus some haplotypes had a very broad distribution occurring in our study system and at least as far away as Idaho, Wyoming, and Colorado. Why are a few haplotypes widely distributed when the vast majority have a restricted distribution? There are at least two related possibilities. First, the most widely distributed haplotypes may be the most abundant haplotypes. That is, abundant haplotypes should have the greatest chance of being transported by vectors. Any vector must come in contact with individuals of P. coronata, and that possibility increases for the most abundant haplotypes. Second, the most widely distributed haplotypes of P. coronata may have expanded their distribution during early post-glacial climates. Cooler temperatures throughout the stream network may have facilitated the dispersal of P. coronata before warming restricted its distribution to upper reaches in tributary segments. Catchments in this study were free of extensive ice sheets starting about 12,000 years ago (e.g. Laabs et al. 2006) whereas, warming and desertification became evident at about 9000 years ago (Rhode 2000). This leaves approximately 3000 years of cooler temperatures for some haplotypes of P. coronata (possibly the most abundant) to expand their distribution.
Species that rarely disperse should show low haplotype diversity at local scales because novel haplotypes would not be able to readily disperse across our study system (Connor & Hartl 2004). Free dispersal would increase the number of haplotypes at any given site as all newly derived haplotypes within the region could colonize all sites. We found, on average, 5 different haplotypes out of 10 individuals per site with a few sites showing as many as 9 out of 10. Since dispersal appears to be rare in this species, how did several haplotypes of P. coronata accumulate at local scales? Over time, isolated populations can accumulate haplotypes by mutation and slow lineage sorting. Lineage sorting can be a slow process in large populations and P. coronata was one of the five most abundant species in our study. As population size increases, the rate at which rare alleles are lost via drift decreases (Wright 1931), whereas the rate at which novel haplotypes are generated via mutation increases, which would result in increased levels of mtDNA polymorphisms (Templeton 2006). A slow rate of lineage sorting also explained a large number of haplotypes at local scales in Nematodes (Blouin et al. 1992) and Blackflies (Finn & Adler 2006), which are also species that attain high population densities and show a low degree of connectivity.
Also, asexual reproduction can increase the effective population size for mtDNA. Only females pass their mtDNA to the next generation in species with obligate sexual reproduction. In populations that reproduce by fission, all of the individuals pass their mtDNA to the next generation, thus preserving a greater proportion of variant mtDNA fragments (Bessho et al. 1992). Both factors (large populations and asexual reproduction) can cause incomplete lineage sorting and prolong the persistence of haplotypes at local scales, which may account for the high number of haplotypes at local scales in P. coronata.
The genus Polycelis has a very uniform morphology potentially hiding several new species (e.g. Kenk 1973). Identifying new species in asexual flatworms can be especially difficult because the morphology of the penis, which can be obscured in asexual species, is a critical diagnostic character (e.g. Kenk 1973). We discovered what appears to be two cryptic species based on mtDNA. Polycelis sp. A was widespread across the Provo Basin, being found at seven sites, although usually in low abundance. However, it was common at one site in the Provo catchment (Marjorie Creek) and at one site in the Weber catchment (Shingle Creek). Otherwise, we collected a single individual of Polycelis sp. A from the American Fork catchment and one individual from Wyoming (Fort Bridger). Only five individuals of Polycelis sp. B were found from the Logan, (UT) Snake (ID), and Upper Colorado (WY) basins all outside our study system.
Two other subspecies of P. coronata are currently recognized, P. c. brevipenis (from western Colorado) and P. c. monticola (from California, near Lake Tahoe); while a third form, P. borealis (described from Alaska), has been synonomized with P. coronata (Kenk 1989). We were not able to include samples from their respective type localities to determine if any of these existing names apply to the cryptic species in our study. However, our data suggest the validity of previous observations (Braithwaite 1962), that continued sampling across the entire range of this species might reveal a large complex of unidentified species currently grouped with P. coronata. Future research should expand the spatial extent of this research to include the entire range of this species. Ignoring cryptic diversity distorts our perception of general diversity patterns and can obscure our understanding of a species niche requirements and functional roles (e.g. Smith et al. 2006).
Polycelis coronata is an obligate aquatic species that inhabits cool, headwater tributaries at the tips of the stream network in catchments of the central Wasatch Mountains. Our intensive, hierarchical sampling design over a very small spatial extent leaves little doubt that this is a poorly connected species. Thus, we found no evidence to suggest that asexual reproduction and recent transport by humans have played a role in increasing the dispersal of this species in the central Wasatch Mountains. However, low connectivity may not be attributed to habitat stability, and multiple explanations should be explored. For example, habitat specialization could also select against a dispersive phenotype. Finally, the effect of dispersal by vectors on the distribution of haplotypes warrants further investigation. Do vectors produce a strong stochastic component to dispersal and how does that effect the distribution of haplotypes? Do individuals at a local site become the greatest source of genetic variation with little structuring at larger scales if vectors are the primary means of dispersal?