Integrative taxonomy resolves taxonomic uncertainty for freshwater mussels being considered for protection under the U.S. Endangered Species Act

Objectively delimiting species boundaries remains an important challenge in systematics and becomes urgent when unresolved taxonomy complicates conservation and recovery efforts. We examined species boundaries in the imperiled freshwater mussel genus Cyclonaias (Bivalvia: Unionidae) using morphometrics, molecular phylogenetics, and multispecies coalescent models to help guide pending conservation assessments and legislative decisions. Congruence across multiple lines of evidence indicated that current taxonomy overestimates diversity in the C. pustulosa species complex. The only genetically and morphologically diagnosable species in the C. pustulosa species complex were C. pustulosa and C. succissa and we consider C. aurea, C. houstonensis, C. mortoni, and C. refulgens to be synonyms of C. pustulosa. In contrast, all three species in the C. nodulata complex (C. necki, C. nodulata, and C. petrina) were genetically, geographically, and morphologically diagnosable. Our findings have important conservation and management implications, as three nominal species (C. aurea, C. houstonensis, and C. petrina) are being considered for protection under the Endangered Species Act.

In our study, we implement an integrative taxonomic approach that utilized multi-locus sequence data, morphometric analyses, and geographic distributions to investigate species boundaries in both the C. nodulata and C. pustulosa species complexes. Our findings support that four geographically isolated taxa (C. aurea, C. houstonensis, C. mortoni, and C. refulgens) are not valid species and are considered here as synonyms of C. pustulosa. Within the C. nodulata complex, we provide extensive molecular, morphological, and biogeographical evidence for C. necki, a recently described species 43 . In order to clarify the distribution, phylogenetic position, and morphological variation, we redescribe C. necki based on our more robust integrative dataset. Our findings have important conservation and management implications and are likely to impact ESA listing decisions for at least three nominal species (C. aurea, C. houstonensis, and C. petrina s.s.).

Results
Taxon Sampling and Molecular Analyses. We generated and analyzed 391 CO1, 391 ND1, and 217 ITS1 DNA sequences for this investigation. Collection details, museum catalog numbers, and GenBank accession numbers are presented in Supplementary Table S1 (also available at https://doi.org/10.5066/P9SRSHV2). Our three-locus molecular matrix consisted of 217 individuals representing 8 genera and 21 recognized species (Table 1). Each taxon was represented by CO1 (avg. ≈ 642 nucleotides [nt]), ND1 (avg. ≈ 797 nt), and ITS1 (951 nt with avg. ≈ 49.13% gaps) and the concatenated three-locus alignment consisted of 2397 nt. Protein coding mtDNA genes did not contain any gaps or stop codons. The large proportion of gaps in the ITS1 alignment was a consequence of partial duplication in the gene region (294-298 nt) found in Cyclonaias tuberculata, which was previously reported 33 . Five partitions and nucleotide substitution models were selected by Partitionfinder for implementation in both IQ-TREE and BEAST: CO1 and ND1 1st position-TrNef+I+G, CO1 and ND1 2nd position-HKY+I+G, CO1 3rd position-HKY+G, ND1 3rd position-TrN+G, and ITS1-K80+I+G. To reduce redundancy of IQ-TREE and BEAST analyses, invariant sites were not modeled in instances when a gamma distribution for rate heterogeneity was estimated 44,45 . Convergence of BEAST runs was supported by ESS > 200 for all parameters. We present the ML phylogenetic reconstruction of the concatenated 3-gene matrix containing ML and BI nodal support values ( Fig. 2; Supplementary Figs S1 and S2; also available at https://doi.org/10.5066/ P9SRSHV2).
Phylogenetic analyses recovered C. nodulata nested between two monophyletic and geographically isolated clades representing C. petrina s.s. and C. necki (Fig. 3). In contrast, five of the six recognized species in the C. pustulosa species complex were not monophyletic in the optimal topology (Fig. 4). Specifically, C. succissa was sister to a clade containing C. aurea, C. houstonensis, C. mortoni, C. pustulosa, and C. refulgens. For the C. nodulata complex, totals of 80 and 55 individuals were included in the mtDNA and ITS1 haplotype networks, respectively (Fig. 3). Three geographically isolated groups were recovered in both networks: C. petrina from the Colorado River basin, C. necki from the Guadalupe River basin, and C. nodulata from the Neches, Ouachita, Red, and Salt river basins. For the C. pustulosa species complex, 263 and 114 individuals were included in the mtDNA and ITS1 haplotype networks, respectively (Fig. 4). Cyclonaias succissa was molecularly diagnosable from other taxa and clearly divergent in both the mtDNA and ITS1 haplotype networks. All other species shared ITS1 haplotypes and showed weak phylogeographic structuring among mtDNA haplotypes.
We observed no overlap between intraspecific variation and interspecific divergence in genetic distance among members of the C. nodulata complex (Fig. 5). Additionally, all three clades contained diagnostic nucleotides: C. petrina (CO1/ND1/ITS1 = 4/16/2), C. necki (CO1/ND1/ITS1 = 10/12/3), and C. nodulata (CO1/ND1/ ITS1 = 6/5/3). However, uncorrected p-distances show a high degree of overlap between intraspecific variation and interspecific divergence among members of the C. pustulosa complex, with the exception of C. succissa (Fig. 5), which also exhibited diagnostic nucleotides (CO1/ND1/ITS1 = 3/4/6). None of the other taxa were molecularly diagnosable. The AMOVA results parallel the levels of genetic distances observed in each species complex ( Table 2). The AMOVA for members of the C. pustulosa complex indicated that genetic variation within species was roughly equal to variation between species, with 52.42% and 51.32% of the variation within species, and 47.58% and 48.68% between all species species for CO1 and ND1, respectively. In contrast, AMOVA between members of the C. nodulata complex revealed high levels of genetic structuring, with 87.45% and 88.98% of the variation between the three species groups and 12.55% and 11.02% within species groups for CO1 and ND1, respectively.
Revised Diagnosis: Specimens of C. necki are distinguished from C. petrina by having a shell that is more elongate and more compressed with less fluting along the posterior slope and a periostracum that is typically more yellow with subdued broken green rays (Fig. 1b). It also lacks the two rows of nodules present on the shell disk of C. nodulata. Cyclonaias necki also has 10 diagnostic nucleotides at CO1 Redescription: Maximum shell length to 69 mm, height to 50 mm, and width to 30 mm. Shell subquadrate to subovate or nearly ovate in outline, often with slight corrugations or parallel ridges on the posterior slope that are sometimes obscured by the accumulation of precipitates on the posterior portion of the shell. Posterior ridge rounded, shells thick, solid, and moderately compressed laterally. Periostracum often with a cloth-like texture, yellow to tan or brown to black in color, sometimes with broken green rays or blotches on the disc or along the posterior slope. Nacre color white, usually iridescent posteriorly. Umbo high and usually extends well above hinge line. Umbo sculpture with 2-4 rows of nodules following a 45° angle relative to hinge line with cross-hatching in younger specimens. Left valve with two thick lateral teeth, straight to slightly curved; two Distribution: Cyclonaias necki is endemic to the Guadalupe River drainage in Central Texas (Fig. 6). The historical distribution of C. necki in the Guadalupe River drainage is known from observations in the Guadalupe   In the San Antonio River drainage, historical records for C. necki are limited to a single specimen (UMMZ 77200) purportedly collected from Salado Creek, a tributary of the San Antonio River. Based on shell morphology, this specimen looks more like C. petrina from the Colorado River than C. necki from the Guadalupe, suggesting the locality information associated with this single specimen may be inaccurate. Furthermore, Strecker 46 did not report C. necki from the San Antonio River basin. In the 1970s, Joseph Bergmann reported shell material, of unknown condition, of C. necki from the Medina River near Von Ormy and Macdona, TX, and Salado Creek near Fort Sam Houston 48 , but these collections have since been lost and so the identifications cannot be verified. Around the same time P. Barker, an amateur naturalist, reported collecting a single specimen of C. necki from Medina Reservoir 48 . This species is not known to occur in lakes or reservoirs so the true collection locality for this specimen remains in question. More recently, investigators reported finding shell fragments thought to be C. necki from the San Antonio River at the San Antonio River Walk 49 , but the weathered condition of these fragments precludes confident identification to any species. Recent fieldwork in central Texas has led to the discovery of live individuals or very recently dead specimens of C. necki in the following rivers within the Guadalupe River drainage: San Marcos River in Caldwell, Guadalupe, Gonzales, and Hays counties; and Guadalupe River in Comal, Gonzales, Kerr, Kendall, and Victoria counties 48,[50][51][52] (Supplementary Table S1).

Common name: Texas Pimpleback
Type Material Holotype: MCZ 169291 (Fig. 1f), length 38 mm, Texas, Llanos River, collected by T. H. Webb. Material Examined: All material examined available as Table S1. Diagnosis: Cyclonaias petrina is distinguished from C. necki by having a shell that is thicker, less elongate, and less compressed with more fluting along the posterior slope and a periostracum that is typically more tan to brown, occasionally with vague green rays or blotches. Unlike C. nodulata, it lacks the two radiating rows of nodules located on the shell disk (Fig. 1).   Table 2. Analysis of molecular variance (AMOVA) among members of the Cyclonaias nodulata and Cyclonaias pustulosa species complexes. Samples were grouped according to current taxonomy. All values were significant (P < 0.0001).

Type Material
Holotype: Original figured type specimen not found, which was designated as a lectotype by Johnson (1974). The length of figured shell in original description reported as about 53 mm (Fig. 1g). Syntypes, ANSP 43058, are from Ohio.
Material Examined: All material examined available as Table S1. Diagnosis: Cyclonaias pustulosa is difficult to distinguish from other species of the genus due to high variability in shell morphology and nacre color (Fig. 1). Cyclonaias pustulosa can be distinguished from C. succissa based on three diagnostic nucleotides at CO1 Redescription: Maximum shell length to 83 mm, height to 73 mm, and width to 52 mm. Shell subquadrate to subovate but usually ovate in outline, posterior margin rounded to truncated, posterior slope broad and flat with or without pustules or corrugations. Posterior ridge rounded, shells thick, solid, and moderately to greatly inflated. Shell disk smooth to heavily pustulated. Periostracum smooth to cloth-like, yellow to tan or brown to black in color, occasionally with a broad, broken green ray or blotches extending from the umbo to the posterior ventral margin. May also have green blotches on the posterior slope. Nacre color typically white, although can be light to dark purple. Umbo high, inflated, oriented anteriorly, and usually extends well above hinge line. Umbo sculpture with two to four course ridges, somewhat more nodulous along posterior ridge. Left valve with two thick, straight to slightly curved lateral teeth; two large, robust, erect, triangular pseudocardinal teeth. Right valve with single, lateral tooth, two pseudocardinal teeth, also triangular. The interdentum is moderately wide. Umbo cavity deep, somewhat compressed, and extending well under the interdentum. Distribution: Cyclonaias pustulosa is known from the eastern reaches of the Great Lakes, Lake St. Clair and Lake Erie, and widespread throughout the Mississippi Basin from southern Minnesota south to Louisiana, and from western New York west to South Dakota 39,56 . Based on our findings, we have expanded the range of C. pustulosa to include Gulf Coast rivers from the Pascagoula Basin in Mississippi west to the Nueces Basin in southwest Texas (Fig. 6).

Discussion
Our primary goal was to implement an integrative taxonomic approach 14,18,57 that utilized multi-locus sequence data, morphometric analyses, and geographic distributions to investigate species boundaries in the genus Cyclonaias. For our assessment, we allowed current taxonomy 35,37,38,40,58 to represent our null hypotheses regarding species-level boundaries. Based on this approach, we identified nine well-supported species-level clades, including two species complexes containing taxa of immediate conservation concern (Figs 3-5).
For the C. nodulata complex, both BI and ML analyses resolved C. nodulata nested between two monophyletic and geographically isolated clades representing C. petrina s.s. and C. necki (Figs. 2 and 3). A clear gap between intraspecific variation and interspecific divergence among the three geographically isolated clades (Fig. 5) was exhibited by mtDNA sequences, indicative of species-level divergence and similar to values reported for several other freshwater mussel species 8,9,22,24,32,[59][60][61] . Sequence divergence at ITS1 was lower relative to both mtDNA genes but consistent with patterns observed in previous studies utilizing these loci 8,9,24 . Morphometric analyses revealed little overlap of C. nodulata, C. petrina, and C. necki. when compared to the C. pustulosa complex (Figs 3  and 4). We further tested species boundaries in the C. nodulata complex by implementing the coalescent-based species delimitation model STACEY, which aligned with other molecular and morphological assessments, recognizing C. necki and C. petrina as distinct evolutionary units without a priori designation.
Previous researchers questioned the validity of taxa in the C. pustulosa complex due to difficulties distinguishing between morphological forms, geographic variants, and distinct species [35][36][37][38][39]46,56,58,[62][63][64] . For example, several species in the C. pustulosa complex in western Gulf of Mexico drainages were considered to be sympatric 36,37 . Most recently, the distributions of nominal species in the C. pustulosa complex in western Gulf of Mexico drainages have been revised with all species considered allopatric 38 . If each of these "species" are allopatric and restricted to distinct drainages, then phylogenetic analysis should have recoveedr individuals from each drainage as monophyletic. All data and analyses provided a congruent signal and credible evidence that current taxonomy overestimates species-level diversity in the C. pustulosa species complex.
Our molecular and morphometric data indicate that current taxonomy overestimates species-level diversity in the C. pustulosa complex. In fact, our data show greater genetic divergence and morphological distinctiveness between C. petrina and C. necki than between all C. aurea, C. houstonensis, C. mortoni, C. pustulosa, and C. refulgens sampled. All five taxa previously recognized as species or subspecies in the C. pustulosa species complex exhibited extensive paraphyly (Fig. 4), with no clear distinction between intraspecific variation and interspecific divergence at mtDNA loci (Fig. 5) or clear signals for diagnosis using morphological characters (Fig. 4). With the exception of C. succissa, relationships among mtDNA haplotypes show weak associations with currently recognized taxonomy and several nominal taxa share ITS1 haplotypes (Fig. 4). Additionally, morphometric analyses illustrated limited ability to distinguish between members of the C. pustulosa complex using shell measurements. Specifically, C. houstonensis, C. mortoni, C. pustulosa, and C. refulgens were indistinguishable. The PCA indicated that both C. aurea and C. succissa were more compressed than other members of the complex (Fig. 4), yet only 74% of individuals identified morphologically as C. aurea were binned correctly, with 25% assigned to C. succissa. Our molecular-based analyses, however, do not support the recognition of C. aurea as a distinct species.
Species within the C. pustulosa complex were not molecularly nor morphologically diagnosable, indicating that current taxonomy is vastly overestimating species diversity in this group. Our STACEY analyses were unable to resolve all currently recognized species in the C. pustulosa complex as monophyletic entities, which was consistent with our BI and ML analyses. This is likely the result of excessive haplotype sharing and limited sequence divergence at both mtDNA and nDNA loci. Despite MSCs being a powerful approach to account for incomplete lineage sorting in multi-locus data 3,17,65 , our models in STACEY do not support the recognition of more than two species in the C. pustulosa species complex.

Implications for Taxonomy and Conservation.
In this study, we used an integrative approach that considered molecular, distribution, and morphology data to evaluate species diversity within Cyclonaias. At the species level, congruence across all lines of evidence indicates that current taxonomy overestimates diversity in the C. pustulosa species complex. Considering the lack of diagnosis across multiple independent lines of evidence, we consider C. aurea, C. houstonensis, C. mortoni, and C. refulgens to be synonyms of C. pustulosa. This expands the distribution of C. pustulosa from the Pascagoula River drainage west to the Nueces River drainage in south Texas (Fig. 6). We do see weak evidence for population structure coinciding to drainage of origin in the C. pustulosa complex; however, revising taxonomy such that three or four species are recognized in the C. pustulosa complex would result in taxa with extremely disjunct distributions and sympatric species that cannot be distinguished morphologically or genetically. This information could be useful if future management actions are considered for populations of C. pustulosa. Our findings may impact ESA listing decisions by resource management agencies considering that two species (C. aurea and C. houstonensis) are synonyms of C. pustulosa, and another species (C. petrina s.s.) contains a cryptic lineage, C. necki, which is redescribed herein. Our evaluation of the description of C. necki reveals several deficiencies in the work of Burlakova et al. 43 . First, their findings were based on only ten specimens from three localities, none of which were collected from the mainstem Colorado or Guadalupe rivers. Thus, their description does not represent the full range of morphological or molecular variation, an important consideration for barcoding 66 . Second, Burlakova et al. 43 report that both C. necki and C. petrina are closely related to C. nodulata; however, no evidence is presented to support this claim and no data or material from C. nodulata were examined. In fact, no phylogeny was presented by Burlakova et al. 43 , which is necessary to support their claim that C. necki has been split from C. petrina. Our findings show that the two species (C. necki and C. petrina) are not sister taxa and reveal that C. necki is distinct from both C. nodulata and C. petrina. The exclusion of a phylogeny fails to capture much about the important evolutionary history of the group and our rigorous phylogenetic approach provides a framework for future evolutionary, ecological, and conservation research. Third, neither the results of the Barcode Index Number designations nor halotype networks were presented despite stating these two methods were used to implement "molecular-based species delimitation, " albeit only from a single mtDNA gene. There has long been concern in the scientific community about taxonomy based on single mtDNA gene [67][68][69] . Empirical rigor is necessary to ensure confidence in species descriptions and taxonomic stability 18 . Our example illustrates the importance of detailed integrative taxonomy when important management decisions require taxonomic clarity.

Methods
Taxon Sampling and Molecular Data. Our taxon sampling concentrated on the following recognized taxa: Cyclonaias asperata (Lea, 1861), C. aurea, C. houstonensis, Cyclonaias infucata (Conrad, 1834), Cyclonaias kleiniana (Lea, 1852), C. mortoni, C. nodulata, C. petrina, C. pustulosa, C. refulgens, and C. succissa. Efforts were made to sample throughout the range of each species with an emphasis on each type locality. Outgroup taxa were selected based on relationships resolved in previous phylogenetic studies 32,33,70 . All specimens involved with DNA analyses were sacrificed for vouchering in museum collections except four individuals of C. houstonensis that were non-lethally swabbed 71 .
Fingings of an undescribed species in the C. nodulata complex (voucher UF440979), initially reported by the lead author, were broadly shared with the scientific and regulatory community 51,55 . After submission of our current manuscript for review, the species was named Cyclonaias necki by Burlakova et al. 43 . We have modified our paper using the name, C. necki, and provide a more thorough description based on quantitative morphometrics, robust phylogenetic analyses, and broad taxonomic and geographic sampling.
We utilized two protein-coding mitochondrial (mtDNA) genes and one nuclear (nDNA) locus for phylogenetic reconstruction: cytochrome c oxidase subunit 1 (CO1), NADH dehydrogenase subunit 1 (ND1), and internal transcribed spacer 1 (ITS1). Tissue samples and DNA swabs were preserved in 95% ethanol and DNA was extracted using a modified plate extraction protocol 72 . Primers used for polymerase chain reaction (PCR) and sequencing were as follows: CO1 dgLCO-1490-GGTCAACAAATCATAAAGAYATYGG and CO1 dgHCO-2198-TAAACTTCA GGGTGACCAAARAAYCA 73 ; ND1 Leu-uurF-TGGCAGAAAAGTGCATCAGATTAAAGC and LoGlyR-CCTGCTTGGAAGGCAAGTGTACT 32 ; ITS1-18S-AAAAAGCTTCCGTAGGTGAACCTGCG and ITS1-5.8S-A GCTTGCTGCGTTCTTCATCG 28 . Thermal cycling profiles for COI were as follows: an initial denaturation at 95 °C for 3 min followed by 5 cycles of 95 32,73 . The PCR protocol for plate amplifications was conducted in a 12.5 µl mixture: distilled deionized water (4.25 µl), MyTaq TM Red Mix (6.25 µl) (Bioline), primers (0.5 µl) and DNA template (20 ng). Bidirectional sequencing was performed at the Interdisciplinary Center for Biotechnology Research at the University of Florida on an ABI 3730 (Life Technologies). Geneious v 9.1.5 74 was used to edit chromatograms and assemble consensus sequences. The mtDNA genes were aligned in Mesquite v 3.2.0 75 using the L-INS-i method in MAFFT v 7.299 76 and translated into amino acids to ensure absence of stop codons and gaps. The ITS1 alignment was performed using the E-INS-i method in MAFFT to better account for the presence of indels following developers' recommendations.
Phylogenetic and Phylogeographic Analyses. We estimated phylogenetic relationships using a concatenated three-locus dataset (i.e. CO1, ND1, ITS1) for members of Quadrulini using maximum likelihood (ML) searches in IQ-TREE v 1.6.1 77,78 and Bayesian inference (BI) in BEAST v 2.4.8 79 . We used PartitionFinder v2.1.1 80 to identify partitions and substitution models for IQ-TREE using BIC and to test all models of nucleotide evolution available for BEAST. ML analyses included an initial tree search before implementing 10000 ultrafast bootstrap (UFBoot) replicates to estimate nodal support and nodes are considered supported if they have a UFBoot >95% 81 . BI analyses executed two runs of 3.75 × 10 8 for a total of 7.5 × 10 8 generations sampling trees every 10000 generations with an initial 25% burn-in. Trace logs and species trees for the two runs were combined using LogCombiner v 2.4.8. Tracer v 1.6 82 was used to calculate the standard deviation of log rate on branches and the coefficient of variance was >0.1 for all partitions 45 ; therefore, a relaxed log-normal molecular clock was used on all partitions. The relaxed log-normal molecular clock for the first partition was fixed at 1.0 and remaining partitions were estimated by BEAST. Yule process was used as the species tree prior and trees were linked across all partitions. To ensure adequate sampling and proper burn-in, effective sampling (ESS > 200) of all parameters was ensured in Tracer. We used SumTrees v 4.4.0 in DendroPy v 4.4.0 83 to estimate a consensus tree with an initial 25% burn-in.
TCS haplotype networks 84 were generated from mtDNA and nDNA independently for each group using PopART 1.7 85 to visualize the geographic distribution of genetic diversity within and between the members of two species complexes (Fig. 1) C. refulgens, and C. succissa) and the C. nodulata species complex (C. necki, C. nodulata, and C. petrina). We included samples lacking ITS1 sequences in the mtDNA haplotype networks to increase sample sizes and expand geographic coverage.
To investigate DNA sequence divergence between and within members of both species complexes, we calculated uncorrected pairwise genetic distances in MEGA7 86 for CO1, ND1, and ITS1 independently. Uncorrected p-distances were chosen over model-based distances because the later have been shown to inflate OTU assignments [87][88][89] . Specimens were identified and sequences were grouped according to drainage of collection following recent taxonomic and distributional assessments 38,40,63 : C. aurea (Guadalupe and Nueces), C. necki (Guadalupe), C. houstonensis (Colorado and Brazos), C. mortoni (Trinity, Neches, and Sabine), C. nodulata (Neches, Red, Sabine, and Mississippi), C. petrina (Colorado), C. pustulosa (Neosho, Ohio, Osage, Ouachita, Red, St. Croix, and St. Francis), C. refulgens (Pascagoula and Pearl), and C. succissa (Escambia, Yellow, and Choctawhatchee) (Fig. 6). Gaps and missing data were treated by pairwise deletion between taxa and each taxon was evaluated for diagnostic nucleotides at each mtDNA locus. Additionally, we conducted an analysis of molecular variance (AMOVA) 90 following 1000 permutations to evaluate inter-and intra-population diversity among members of both the C. pustulosa and C. nodulata species complexes using ARLEQUIN v 3.5 91 . These groupings align with the null hypothesis based on current taxonomy 35,37,38 and were not based on distinct genetic groups or phylogeographic results.
Morphometric Analyses. We collected morphometric data for members of the C. pustulosa and C. nodulata species complexes by measuring external shell dimensions on all specimens used in genetic analyses and individuals encountered during field surveys. Three morphological measurements were made to the nearest 0.01 mm using digital calipers: maximum length, height, and width. Measurement values were log e -transformed to produce a scale-invariant matrix while preserving information about allometry 92,93 . Log e -transformed variables were converted into three ratios: height/length, width/length, and width/height. We examined morphological variation using principal components analyses (PCA) in the ggbiplot package 94 and canonical variates analyses (CVA) in the package Morpho 95 using R v 3.3.1. The PCA analyses were performed to visualize whether morphological groupings were apparent without a priori assignment to a specific group. Canonical variate scores were used for cross-validated discriminant analyses (DA) to measure how reliably morphometric data could assign individuals to each a priori defined species in the C. nodulata or C. pustulosa complex.

Coalescent-based Species Delimitation.
We implemented MSC models using the STACEY v 1.2.2 package 96 in BEAST v 2.4.8 79 on two independent molecular matrices coinciding with each species complex. Additionally, we included Cyclonaias tuberculata (Rafinesque, 1820) and C. succissa as outgroups for the C. nodulata and C. pustulosa complexes, respectively. The refined matrices were realigned to collapse gaps in ITS1 caused by outgroup taxa. Partitions and substitutions models for each matrix were reevaluated using Partitionfinder v2.1.1 80 . MSC models implemented in STACEY improve efficiency of the DISSECT birth-death collapse model 96 and infers species boundaries without a priori species designations. Therefore, we allowed the model to consider all individuals as separate species and assign individuals to species clusters. Each run executed 3 × 10 8 generations and logged every 5,000th tree with an initial 10% burn-in. Adequate sampling and proper burn-in was ensured using Tracer. The number of well-supported clusters was calculated using SpeciesDelimitationAnalyser 96 following an initial 10% burn-in (6,000 trees).

Data Availability
All molecular, morphological, and geographic information generated and analyzed as part of this study are publicly available at https://doi.org/10.5066/P9SRSHV2. We included DNA alignment files, morphological measurements, museum catalog numbers, and collection information for every specimen. All DNA sequences analyzed in this study are novel and available on GenBank (MH361762-MH362664, MH633560-MH633655).