Taxonomic revision of the Pyrgulopsis gilae (Caenogastropoda, Hydrobiidae) species complex, with descriptions of two new species from the Gila River basin, New Mexico

Abstract We describe two new species of springsnails (genus Pyrgulopsis) for populations from the middle Fork and upper East Fork of the Gila River Basin (New Mexico) that had been previously identified as P. gilae. We also restrict P. gilae to its originally circumscribed geographic range which consists of a short reach of the East Fork Gila River and a single spring along the Gila River (below the East Fork confluence). These three species form genetically distinct lineages that differ from each other by 3.9–6.3% for mtCOI and 3.7–8.7% for mtNDI (the latter data were newly obtained for this study), and are diagnosable by shell and penial characters. Collectively the three species form a strongly supported clade that is distinguished from other congeners by the unique presence of two glandular strips on the dorsal surface of the penial filament. These findings suggest that the conservation status of P. gilae, which was recently removed from the list of candidates for listing as endangered or threatened by the United States Fish and Wildlife Service, should be revisited and that the two new species may also merit protective measures given their narrow geographic ranges.


Introduction
Pyrgulopsis is a large genus (137 species ;Hershler et al. 2013) of freshwater gastropods that is distributed in North America west of the Mississippi River basin. The tiny species in this genus live in spring-fed habitats and usually have very small geographic ranges. Pyrgulopsis is a current focus of conservation efforts owing to threats posed by groundwater extraction, livestock grazing and other anthropogenic activities ). Recent molecular studies have shown that several congeners are composites of genetically divergent lineages and are in need of taxonomic revision (e.g., Liu et al. 2003, Hershler and Liu 2010. This is the second in an anticipated series of papers that clarifies the taxonomy of these species . Pyrgulopsis gilae (Taylor, 1987) was described for specimens from single springs along the lower East Fork (type locality) and main stem Gila River in Grant County, New Mexico. Field surveys in the 1990's and 2000's resulted in the discovery of new populations in two other reaches of the upper Gila River watershed (Middle Fork, upper East Fork) that are currently being treated as P. gilae (NMDGF 2012). Hurt (2004) delineated substantial divergence in mtCOI sequences (6.8% average) between specimens from the upper East Fork reach (Wall Spring) and three localities within the originally circumscribed range of P. gilae. In a more comprehensive survey of COI variation within P. gilae, populations from the upper East Fork, Middle Fork, and lower East Fork (and main stem Gila River) reaches were resolved as three divergent (3.9-6.3%) sub-clades which were postulated to be distinct species ). Here we document a congruent pattern of variation in a second mitochondrial DNA marker (NDI) and delineate morphological differences supporting recognition of the upper East Fork and Middle Fork Gila River populations as new species, which are described herein.

Methods
For the current molecular study we used the same samples that were analyzed in our previous phylogeographic investigation of P. gilae across its entire geographic range ; Fig. 1). Genomic DNA was extracted from single, entire snails using a CTAB protocol (Bucklin 1992); 3-8 specimens were analyzed (separately) from each sample. ND43F and RND592F (Liu et al. 2003) were used to amplify a 530 bp fragment of NADH dehydrogenase subunit I (NDI). This primer pair did not amplify the region for specimens from two localities (G8, G11) and consequently we designed a second set of oligonucleotide primers for these snails, ND30 (5'TCT TAY ATR CAR ATW CGT AAA GG3') and RND490 (5'ATG TTA CAA ATC ATA TAA ATG3'), based on conserved regions of NDI in an alignment from P. gilae and two closely related species (P. deserta [Pilsbry], P. davisi [Taylor]). Degenerate positions are represented by the following ambiguity codes: Y=C/T; R=A/G; W=A/T. Amplification conditions and sequencing of amplified polymerase chain reaction product followed Liu et al. (2003). Sequences were determined for both strands and then edited and aligned using SEQUENCHER™ version 5.0.1.
The 56 newly sequenced specimens of P. gilae were analyzed both separately and together with our previously published COI dataset . We included the same set of outgroup taxa as in our prior study of P. gilae , with Floridobia floridana again being used as the root. The GenBank accession numbers for these sequences (COI, NDI) are given in Appendix 1. Note that we newly sequenced specimens of P. "mimbres" for COI (using the methods of Liu et al. 2003) and NDI as part of this study (GenBank accession numbers: COI, KM205358; NDI, KM205359). The newly obtained haplotypes from each P. gilae sampling locality were deposited in GenBank. Sample information and GenBank accession numbers are given in Table 1. One example of each haplotype detected in a given sample was used in our analyses. The partition homogeneity/incongruence length difference test (Farris et al. 1994) was used to determine whether the COI and NDI datasets were consistent and could be combined for the phylogenetic analysis. This test, which was conducted using parsimony-informative sites only and 1,000 replicates, did not detect significant incongruence (P=0.21) and consequently we combined the two datasets in the phylogenetic analysis. MRMODELTEST 2.3 (Nylander 2004) was used to obtain an appropriate substitution model (using the Akaike Information Criterion) and parameter values for this analysis. Phylogenetic relationships were inferred by Bayesian analysis using MR-BAYES 3.1.2 (Huelsenbeck and Ronquist 2001). Metropolis-coupled Markov chain Monte Carlo simulations were run with four chains (using the model selected through MRMODELTEST) for 2,000,000 generations, and Markov chains were sampled at intervals of 10 generations to obtain 200,000 sample points. We used the default settings for the priors on topologies and the GTR + I + G model parameters selected by MRMODELTEST as the best fit model. The Tracer program was used to analyze runs for Effective Sample Size (ESS, greater than 200) to ensure that sufficient sampling occurred. At the end of the analysis, the average standard deviation of split frequencies was less than 0.01 (0.002) and the Potential Scale Reduction Factor (PSRF) was 1, indicating that the runs had reached convergence. The sampled trees with branch lengths were used to generate a 50% majority rule consensus tree with the first 25% of the samples removed to ensure that the chain sampled a stationary portion. Genetic relatedness within P. gilae was further assessed by a haplotype network that was generated by TCS version 1.21 using the default settings (e.g., 95% connection limit) and fixing the connection limit at 90 steps (Clement et al. 2000). NDI sequence divergences (maximum composite likelihood) within and between lineages were calculated using MEGA6 (Tamura et al. 2013), with standard errors estimated by 1,000 bootstrap replications with pairwise deletion of missing data. Structuring of variation among lineages was evaluated by an AMOVA using Arlequin 3.5 (Excoffier and Lischer 2010).
Types and other voucher material were deposited in the National Museum of Natural History (USNM) collection. Specimens of P. gilae from the Bell Museum of Natural History (BellMNH) were also examined during the course of this study. Series of large adults (n>10) were used for shell measurements. Whorl counts refer to the entire shell. Sexual dimorphism in shells, which is occasionally observed in Pyrgulopsis (Taylor 1987), could not be quantified owing to small sample sizes. The total number of shell whorls was counted (WH) for each specimen; and the height and width of the entire shell (SH, SW), body whorl (HBW, WBW), and aperture (AH, AW) were measured from camera lucida outline drawings using a digitizing pad linked to a personal computer (see Hershler 1989). In addition, three ratios were generated from the raw data (SW/SH, HBW/SH, AH/SH). Descriptive statistics were generated using SYSTAT FOR WINDOWS 11.00.01 (SSI 2004). T-tests (two-tailed) of differences among shell variables were conducted using an on-line calculator (http:// in-silico.net/tools/statistics/ttest); data for type material of P. gilae were from Taylor (1987, table 11). Penial variation was described from series of adult specimens that were relaxed with menthol crystals and fixed in dilute formalin prior to preservation in 70% ethanol. Descriptive penial terminology is from Taylor (1987) and Hershler (1994Hershler ( , 1998. Variation in the number of cusps on the radular teeth (n=5) was assessed using the method of Hershler et al. (2007).
We used a conservative, evolutionary lineage concept in describing new species only for those snails that are morphologically diagnosable as well as phylogenetically independent and substantially divergent genetically ). Inasmuch as the principal goal of our paper was to delimit species, we provide only brief taxonomic descriptions which focus on those aspects of morphology that have proven most useful in previous such studies of Pyrgulopsis (Taylor 1987, Hershler 1994, Hershler 1998.

Results
Sixteen (16) NDI haplotypes of P. gilae were detected, 11 of which were restricted to single populations ( Table 2). The others were shared by pairs of populations along the lower East Fork (haplotype II), upper East Fork (haplotypes VII, IX, XI) and Middle Fork (haplotype XIII) Gila River. Six (6) samples each contained a single haplotype (G3, G4, G5, G7, G13, G14). The TCS analyses (not shown) recovered three well differentiated haplotype groups composed of specimens from along the lower East Fork and main stem Gila River (clade I), Middle Fork Gila River (II), and the upper East Fork Gila River (III). These groups differed from each other by 3.7-8.7% sequence divergence; variation within groups was minor ( Table 3). The AMOVA indicated that most of the detected variation (91.7%) was partitioned among these groups; variation within populations, and among populations within the groups was much smaller (1.35, 6.93%) but nonetheless was significant ( Table 4). The three previously reported clades (I-III; Liu et al. 2013) were similarly recovered in Bayesian analyses of both the NDI dataset, and the combined COI + NDI dataset (Fig. 2). Based on the genetic evidence of distinctiveness and the diagnosable shell and penial characters that are detailed below we recognize two of these lineages as new species which are described herein (clade II as P. marilynae, clade III as P. similis) and restrict P. gilae (clade I) to its originally circumscribed geographic range.   Diagnosis. Distinguished from P. gilae and the species described next (P. similis) by its narrower shell (mean shell width/shell height 0.613 vs. 0.682, t=-9.6588, df=36.2176, P<0.0001, n=30 for P. gilae; 0.613 vs. 0.734, t=-16.3617, df=18.9656, P<0.0001, n=11 for P. similis), more pronounced whorl shoulders, and broad overlap of the ventral surface of the penis by the terminal gland (probably reflecting fusion with a distal ventral gland). Further differs from P. gilae in its smaller size (mean shell height 2.77 vs. 3.47 mm, t=-11.3848, df=21.9544, P<0.0001) and (basal) extension of the outer penial gland to mid-line or left edge of penis. Further differs from P. similis in its larger size (mean shell height 2.77 vs. 2.36 mm, t=7.3691, df=15.3701, P<0.0001), smaller number of dorsal glands on the penis, and larger size of the terminal and ventral glands on the penis.
Operculum (Fig. 3C-D) as for genus; edges of last 0.5 whorl weakly frilled on outer side; portion of muscle attachment margin thickened on inner side. Radula (Fig.  3E-G) as for genus; dorsal edge of central teeth concave, lateral cusps four-five, basal cusp one. Lateral teeth having two-three cusps on both inner and outer sides. Inner marginal teeth with 14-20 cusps, outer marginal teeth with 17-22 cusps. Radula data are from USNM 1135067.
Penial filament and penial lobe about equal in length ( Fig. 4A-B). Filament having two (penial) glands on dorsal surface; inner gland shorter. Outer penial gland curving to mid-line (10/24 specimens) or left edge of penis (14/24 specimens), the latter condition probably represents fusion with a gland on the left edge (Dg2). Terminal gland elongate, horizontal, broadly overlapping ventral surface of penis. Dorsal surface of penis  having gland along right edge of lobe (Dg3) and 2-3 additional glands (22/24 specimens); one specimen did not have any additional glands and one specimen had four additional glands. Ventral gland positioned near centrally. Penial data are from USNM 1135067.

Etymology. The specific epithet is a patronym honoring Marilyn Myers (United States Fish and Wildlife Service, retired) for her dedicated efforts to survey Pyrgulopsis habitats in the upper Gila River basin.
Distribution. A series of seeps and springs along the north side of short reach (ca. 0.25 km) of the Middle Fork Gila River just below Jordan Hot Spring (Fig. 1). The type locality is a seep wall which is the lower-most occurrence of P. marilynae along the Middle Fork Gila River; the water temperature at this site was 25°C on 1 October 2009.
Remarks. Pyrgulopsis marilyane was resolved as sister to P. gilae (100% posterior probability) in the molecular phylogenetic analysis (Fig 2). The apparent fusion of the terminal and distal ventral glands of the penis that characterizes this species (in part) was previously reported for P. sadai (Hershler 1998 , fig. 39I). The sample attributed to Jordan Hot Spring (USNM 883175) may have been collected instead from a closely proximal spring as P. gilae has not been found at the former locality during recent surveys (USFWS 2011b). Diagnosis. Differs from P. gilae in its smaller size (mean shell height 2.36 vs. 3.47 mm, t=--22.7297, df=36.4071, P<0.0001, n=30 for P. gilae), larger number of glands on the dorsal surface of the penis, frequent extension of outer penial gland and/or Dg2 to the mid-line of the penis, and smaller size of the terminal and ventral glands on the penis. Contrasted with P. similis above.
Operculum (Fig. 5C-D) as for genus; edges of last 0.5 whorl frilled on outer side; inner side near smooth. Radula (Fig. 5E-G) as for genus; dorsal edge of central teeth concave, lateral cusps four-six, basal cusp one. Lateral teeth having two-three cusps on inner sides and two-four cusps on outer sides. Inner marginal teeth with Penial filament longer than lobe (Fig. 4C-D). Filament having two (penial) glands on dorsal surface; inner gland shorter. Outer penial gland sometimes extending (basally) to mid-line (4/30 specimens) or left edge (7/30 specimens); Dg2 sometimes curving (basally) to mid-line (11/30 specimens). Terminal gland transverse, rather small. Dorsal surface of penis having gland along right edge of lobe (Dg3) and 3-7 additional glands (30/30 specimens) which form long, slightly oblique strips. Ventral gland small, positioned near centrally; second gland rarely present (4/30 specimens). Penial data are from USNM 1135065. Distribution. Springs along a short reach (ca. 10 km) of the East Fork Gila River from just above Wall Lake to slightly above the mouth of Burnt Corral Canyon (Fig. 1). The type locality is a spring brook (ca. one m wide and 0.25 m deep) that discharges at the base of the canyon wall along the east side of Beaver Creek; the water temperature at this locality was 22.1°C on 21 May 2009. The flow at this locality is augmented by numerous small seeps.
Etymology. The specific epithet is an adjective referring to the close resemblance between this species and both P. gilae and P. marilynae.
Remarks. Pyrgulopsis similis was resolved as sister to the clade composed of P. marilynae and P. gilae (100% posterior probability) in the Bayesian analysis of molecular data (Fig. 2).
Remarks. Examination of the large series of penes that Taylor scored for this species (BellMNH 20898, BellMNH uncat.) indicated that neither the outer penial gland nor Dg2 extends appreciably onto the dorsal surface of the penis (Fig. 4E-F; also see Taylor 1987, fig. 7b-c) in contrast with P. marilynae and P. similis. Specimens from the two new populations (Fig. 6A), which are closely proximal to the type locality, and the disjunct "Alum Spring" population ( Fig. 6B) conformed to P. gilae in all morphological details. Pyrgulopsis gilae co-occurs sympatrically with P. thermalis (Taylor) at several localities (Taylor 1987).

Discussion
The results of this study provide additional evidence that the current taxonomy of Pyrgulopsis in some cases masks cryptic species diversity. Our previous revision of widely ranging P. micrococcus revealed this taxon to be a polyphyletic composite of five species, three of which were undescribed . Here we have shown that, in contrast, P. gilae (in the broad sense) is a monophyletic species complex diagnosed by a unique penial character-the presence of two glandular strips on the dorsal surface penial filament. (Note that P. merriami [Pilsbry and Beecher], which is distributed in an isolated basin in southeastern Nevada, has a somewhat different pattern consisting of two glands on the dorsal and one gland on the ventral surface of the filament; Hershler 1994). These disparate findings underscore the complexity of and taxonomic challenges posed by the Pyrgulopsis radiation, which is characterized by endemism on very fine geographic scales and extensive morphological homoplasy (Hershler 1994, Liu andHershler 2005).
The delineation of cryptic species complexes often has important consequences for conservation (Bickford et al. 2007). Pyrgulopsis gilae was recently removed from the federal candidate list for listing as endangered or threatened in part owing to the discovery of populations from along the East Fork and Middle Fork of the Gila River (USFWS 2011a, USFWS 2011b) that are assigned herein to two new species. The resulting restriction of P. gilae to its originally circumscribed geographic range-several groups of springs in the lower reach of the East Fork Gila River (below the mouth of Black Canyon) and a single spring along the Gila River ca. 2 km below the East Fork confluence-suggests the conservation status of this species should be re-visited by the USFWS. The narrow endemism of the two new species suggests that these may also merit consideration for possible listing by the USFWS.
Our findings also underscore the need for additional field surveys to further delineate the occurrences of Pyrgulopsis in New Mexico and to supplement the recent monograph by Taylor (1987). Large portions of the Gila River and other drainage basins in the state have yet to be carefully searched for these tiny animals.