Eastern North America is home to the most species-rich non-tropical freshwater fish fauna on Earth (Lundberg et al., 2000). Starting with the work of Rafinesque (1820), the discovery and delimitation of North American freshwater fish species relied entirely on external morphology such as numbers of fin elements and scale rows, proportional measurements, and patterns of pigmentation and coloration. Molecular phylogenies have resolved the relationships of the lineages that comprise the rich North American freshwater fish fauna (Near et al., 2004, 2011; Wright et al., 2012; Ghedotti and Davis, 2017; Bagley et al., 2018; Schönhuth et al., 2018; Near and Kim, 2021), but have also led to the discovery of new species and elevation of species long considered synonyms (Wood and Raley, 2000; Wood et al., 2002; Near, 2008; Near et al., 2017). North American freshwater species with large geographic ranges are prime candidates for undiscovered species diversity and the elevation of synonymized species (Piller et al., 2008; Piller and Bart, 2017). Molecular phylogenetic analyses provide a framework of lineage and species diversification, which aids in partitioning morphological trait variation potentially revealing important differences among lineages that are not apparent without a phylogenetic perspective (Berendzen et al., 2009).

The Stripetail Darter, Etheostoma kennicotti (Putnam), is a widely distributed species forming a circle around the Nashville Basin (Fig. 1). The species occurs throughout the Tennessee River, in the Cumberland River system above Cumberland Falls and below the falls in the Laurel River and portions of the Big South Fork, the Green River system, and tributaries of the lower Ohio River in southern Illinois and northwestern Kentucky (Page and Smith, 1976; Smith, 1979; Burr and Warren, 1986; Etnier and Starnes, 1993; Tiemann and Sherwood, 2011). Etheostoma kennicotti was described from specimens collected by Robert Kennicott from a tributary of the Ohio River in Union County, Illinois (Putnam, 1863; Collette and Knapp, 1966). Later, Jordan and Swain (1883) described Etheostoma cumberlandicum from Wolf Creek in the Clear Fork system, a tributary of the Cumberland River above Cumberland Falls. Curiously, Jordan and Swain (1883) considered E. cumberlandicum closely related to E. flabellare with no reference to Putnam's (1863) E. kennicotti, which was known at the time only from its type locality in southern Illinois. Later Jordan and Evermann (1898: 1098) placed E. kennicotti into the synonymy of E. flabellare, but treated E. cumberlandicum as a subspecies of E. flabellare. Forbes and Richardson (1909: 311–313) identified populations of E. kennicotti in Illinois as Etheostoma obeyense Kirsch; however, the species was resurrected from synonymy with E. flabellare in a list of the fishes of Illinois and a contribution on the type specimens of the species (Smith, 1965; Smith and Page, 1975). Page and Smith (1976) interpreted morphological variation in E. kennicotti and E. cumberlandicum as forming an east–west cline, reflecting hypothesized dispersal westward from an ancestral area centered in the Cumberland River drainage. As a result of the analysis of meristic trait variation, Page and Smith (1976) placed Etheostoma cumberlandicum into the synonymy of E. kennicotti.

Fig. 1

Geographic distribution of Etheostoma kennicotti sensu lato and Etheostoma cumberlandicum and sampling locations of specimens used in molecular and morphological analyses. The populations in the Green River system, the Laurel River system, and the Tennessee River system (partitioned between lower and upper Tennessee River) are each considered undescribed species currently listed as E. kennicotti.

In this study, we reexamine the systematics of Etheostoma kennicotti sensu lato (s.l.) using phylogenetic analysis of the mitochondrial cytochrome b gene, genetic variation at 25 loci, and disparity in meristic traits. Our analyses include specimens sampled throughout the geographic distribution of E. kennicotti (s.l.) and indicate Etheostoma cumberlandicum is a distinct species. Additionally, we identify three newly discovered and undescribed species masquerading as E. kennicotti. The relationships of species within the Etheostoma kennicotti complex are discussed in the context of the pre-Pleistocene configuration of rivers in the Eastern Highlands and patterns of speciation that involve species endemic to the upper Cumberland River.

MATERIALS AND METHODS

Specimen collection.—Specimens of Etheostoma kennicotti (s.l.) were collected from shallow riverine habitats using a minnow seine from a time period extending from 1999 to 2019. Specimens from the main channel of the Green River in Green County, Kentucky and in the mainstem of the Tennessee River were sampled using a boat-assisted benthic trawl (i.e., mini-Missouri trawl; Herzog et al., 2005). Specimens were anesthetized in MS-222 and a portion of the right pectoral fin was dissected at the base and preserved in 95% non-denatured ethanol and cataloged into the Yale Fish Tissue Collection (YFTC). Whole-body specimens were fixed in ∼10% buffered formalin for seven to 14 days, rinsed and soaked in water for one to three days, soaked in 50% ethanol for five to seven days before long term preservation in 70% ethanol. Specimens sampled from the main channel of the Green River were not fixed in formalin. An additional 15 field-collected specimens of E. kennicotti (s.l.) that measured less than 20 mm standard length were preserved whole in 95% non-denatured ethanol, were cataloged into the YFTC, and were available only for molecular analyses. All specimens of E. kennicotti and E. cumberlandicum collected during this study were deposited in the research collections of fishes at Austin Peay State University (APSU), the Illinois Natural History Survey (INHS), the Peabody Museum, Yale University (YPM), Tulane University (TU), and the David A. Etnier Collection of Fishes, University of Tennessee (UT).

Molecular analyses.—Sampling locations of specimens used for genetic analyses are listed in Table 1 and shown in Figure 1. DNA was extracted from 95% ethanol-preserved tissues using a standard DNeasy Qiagen Blood and Tissue Kit (Qiagen, Valencia, CA, USA). To minimize downstream enzymatic inhibition, we purified DNA extractions with an ethanol precipitation: 3 M sodium acetate (pH = 5.2) was added equal to 10% of the total volume of the DNA extraction followed with 100% ethanol equal to 2.5 times the total volume of DNA. After mixing, extractions were incubated for ten minutes at –80°C. Samples were centrifuged for 30 minutes at 8,000 RCF, the supernatant was carefully poured off, and the DNA pellet was washed with 250 µL of cold 70% ethanol. Samples were centrifuged again for five minutes at 8,000 RCF, supernatant was poured off, the pellet was allowed to air dry for ∼15 minutes, and the DNA pellet was resuspended with the desired amount of DNAse-free water.

Table 1

Specimens of Etheostoma kennicotti and E. cumberlandicum used in genetic analyses. Institutional abbreviations follow Sabaj (2020).

The phylogenetic relationships of populations of Etheostoma kennicotti (s.l.) and E. cumberlandicum were assessed with phylogenetic analyses of the mitochondrial encoded cytochrome b (cytb) gene. The molecular phylogenetic analysis included 167 specimens of E. kennicotti (s.l.) and E. cumberlandicum. Two specimens of E. flabellare were sampled as outgroup taxa (Table 1). The cytb gene was amplified using previously published PCR primers and cycling conditions (Near et al., 2000). Amplification products were prepared for DNA sequencing using a polyethylene glycol precipitation. Contiguous sequences were assembled from individual DNA sequencing reactions using the computer program Geneious v.7.2 (Kearse et al., 2012). New cytb sequences were aligned by eye to those previously generated in early studies of darter phylogeny (Porterfield et al., 1999; Near et al., 2011). The optimal data partitioning scheme, among the three codon positions of the cytb gene, and molecular evolutionary models were determined using the Bayesian information criterion in the computer program PartitionFinder v. 2.1 (Lanfear et al., 2017). The mitochondrial gene tree was inferred from the aligned cytb sequences using the optimal molecular evolutionary models and partitioning scheme using the computer program MrBayes v. 3.2 (Ronquist et al., 2012), where posterior probabilities for the phylogeny and parameter values were estimated using Metropolis-couple Markov chain Monte Carlo (MC3; Larget and Simon, 1999; Huelsenbeck et al., 2001). The MrBayes analysis was run for 10⁷ generations with two simultaneous runs each with four chains. Convergence of the MC3 algorithm and stationarity of the chains was assessed by monitoring the average standard deviation of the split frequencies between the two runs, which was less than 0.005 after 3x10⁶ generations. In addition, the likelihood score and all model parameter estimates were plotted against the generation number to determine when there was no increase relative to the generation number in the computer program Tracer v. 1.5 (Drummond and Rambaut, 2007). The first 50% of the sampled generations were discarded as burn-in, and the posterior phylogeny was summarized as a 50% majority-rule consensus tree. All cytb gene sequences generated for this study are available at GenBank MZ078304–MZ078463. The DNA sequence alignments, nexus files used for MrBayes 3.2 analysis, the MrBayes 3.2 posterior tree file, and the summarized consensus tree are available at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg).

A relaxed molecular clock analysis of cytb DNA sequences was performed to estimate the divergence times between the lineages resolved within Etheostoma kennicotti. Because the speciation branching model employed in the analysis assumes that only a single specimen is sampled for each species, the cytb sequence of a single specimen was subsampled from six lineages resolved in the MrBayes inferred cytb gene tree: Green River, upper Tennessee River, lower Tennessee River, Ohio River–Clarks River, the upper Cumberland–Big South Fork, and Laurel River. A single E. flabellare sampled from the Middle Fork Vermilion River was included as an outgroup taxon to root the time calibrated phylogeny (Table 1). The uncorrelated lognormal (UCLN) model of molecular evolutionary rate heterogeneity was implemented using the computer program BEAST v. 1.8.3 (Drummond et al., 2006, 2012). The optimal partitioning scheme and molecular evolutionary models were those used in the MrBayes analysis. A Yule pure birth diversification prior was used for the branching rates in the phylogeny. The relaxed molecular clock was calibrated with a prior on the rate of evolution of the cytb gene of 8.99x10^–3 (95% credible interval: 7.46x10^–3, 1.06x10^–2) substitutions per million years estimated in a previous divergence time analysis of darters (Near et al., 2011: table 2). In the xml file, the molecular evolutionary rate was set with a normal distributed prior with a mean of 0.0089 and a standard deviation of 5.0x10^–4. The BEAST analysis was run three times, with each run consisting of 10⁷ generations. Convergence of model parameter values and estimated node ages to their optimal posterior distributions was assessed by plotting the marginal posterior probabilities using the computer program Tracer v. 1.6 (Drummond and Rambaut, 2007). The resulting trees and log files from each run were combined using the computer program LogCombiner v. 1.83. The posterior probability density of the combined tree and log files was summarized using TreeAnnotator v. 1.8.3. The mean and 95% HPD estimates of divergence times were visualized on the chronogram using the computer program FigTree v. 1.4. The input xml file, a sampling of the posterior trees, and the summarized consensus time tree are available at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg).

We genotyped 164 specimens of Etheostoma kennicotti for 25 microsatellite loci using primers developed for this species by Hereditec (Lansing, New York; Table 2). Genotyping was performed by a single-reaction nested PCR method (Schuelke, 2000). We used a forward primer with a universal M13 (5′–TCCCAGTCACGACGT–3′) tail at its 5′ end, a complementary M13 forward primer labeled with one of three fluorescent dyes (6FAM, VIC, or NED), and an unlabeled reverse primer. Microsatellite PCR reactions were performed in a total volume of 12.5 µl, containing 1 µL of isolated DNA, 1X colorless GoTaq Flexi Buffer (Promega Corporation), 2.0 mM MgCl₂, 0.2 mM of each dNTP, 0.04 mM of both the reverse primer and the M13 tagged forward primer, 0.64 mM of the fluorescent-labeled M13 primer, and 0.06 mM of GoTaq DNA polymerase (Promega Corporation). We used an Eppendorf Mastercycler thermo cycler under the following PCR protocol: initial denaturation at 94°C (5 min), followed by 35 cycles of denaturing (94°C, 40 s), annealing, and extension (72°C, 40 s), with the annealing temperature maintained at 57°C (30 s) for the first ten cycles and subsequently lowered to 53°C (30s) for the final 25 cycles. All microsatellite loci were genotyped on 3730xl DNA Analyzer (Applied Biosystems, Inc.) against a LIZ-500 dye size standard (Applied Biosystems, Inc.). We used GENEMAPPER v3.7 (Applied Biosystems, Inc.) to retrieve raw allele sizes and scored them by using the automatic binning function in TANDEM (Matschiner and Salzburger, 2009).

Table 2

Microsatellite primer sequences used in genotyping specimens of Etheostoma kennicotti and E. cumberlandicum.

We tested for departures from Hardy-Weinberg equilibrium (HWE) using the “hw.test” function in the R package Pegas v.0.14 ( http://ape-package.ird.fr/pegas.html). The proportion of polymorphic loci (P) per population was calculated with “locus_table” function in the R package poppr v. 2.8.6 ( https://github.com/grunwaldlab/poppr/). The average number of alleles per locus in each population (A) and private allelic richness (pAR) were calculated using HP-RARE v. June-6-2006 (Kalinowski, 2005), using 12 loci for rarefaction. The number of private alleles per population was calculated with the “private_alleles” function in poppr v. 2.8.6. The mean population expected (H_e) and observed heterozygosity (H_o) was calculated with the “divBasic” function in the R package diveRsity v. 1.9.0 ( https://rdrr.io/cran/diveRsity/). Pairwise population Weir and Cockerham (1984) FST indices were calculated using the “genet.dist” function in the R package hierfstat v. 0.5-7 ( https://github.com/jgx65/hierfstat).

We assessed population structure of Etheostoma kennicotti using the Bayesian clustering algorithm STRUCTURE v2.3.2 (Pritchard et al., 2000; Falush et al., 2003). We used a hierarchical ΔK method to infer the number of genetic clusters by repeating STRUCTURE analyses on each of the K groups inferred in the previous step (Coulon et al., 2008). Hierarchical clustering continued until the log-likelihood for one cluster (K = 1) was larger than the log-likelihoods for all other values of K, or the majority of individuals were not strongly assigned to any cluster (assignment probability < 0.6). In each round of analysis, STRUCTURE was run from K= 1 to K = 10, with ten independent iterations performed for each value of K. The burn-in period was 100,000 replicates followed by 1,000,000 Markov chain Monte Carlo simulations run under a model that assumed admixture and correlated allele frequencies, without prior information on population identity. The optimal K was calculated using the ΔK method implemented in Structure Harvester v0.6.94 (Earl and vonHoldt, 2012). The scored microsatellite alleles used for the analyses are available at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg).

Morphological analyses.—Meristic data were collected from 703 specimens of Etheostoma kennicotti sampled in our fieldwork as well as specimens housed at UT. Catalog numbers of specimen lots used for morphological analyses are given in Material Examined. Detailed sampling locations are shown on Figure 1 and georeferenced data for all sampling locations are available at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg). The numbers of scale rows and fin elements were determined from each specimen as outlined in Hubbs and Lagler (1958) and Page (1981), with the exception of the number of transverse scale rows, which was counted as described by Page (1983: 16, fig. 2). The characterization and counting of caudal bands follows Page and Smith (1976). The meristic data for all 703 specimens are available at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg).

A principal component analysis (PCA) of the meristic traits was executed using the “prcomp” function R version 3.2.0 (R Core Team, 2015). The ability of the meristic data to assign individual specimens to one of six groups based on relationships inferred in phylogenetic analysis of mtDNA gene sequences or genomic clusters estimated using microsatellite loci was assessed using a cross-validation linear discriminate analysis (LDA) as applied in the MASS package for R (Venables and Ripley, 2002). Bayesian posterior assignment probabilities for each group were calculated, with group assignment determined by the highest posterior probability. The meristic traits analyzed with PCA and LDA included: number of lateral line scales, number of transverse scale rows, number of scales around the caudal peduncle, number of first dorsal-fin spines, number of second dorsal-fin rays, number of anal-fin rays, and number of pectoral-fin rays. We calculated pairwise group mean Mahalanobis distances in full principal component space using the pairwise.mahalanobis function from the HDMD package (McFerrin, 2013).

RESULTS

Molecular analyses.—The PartitionFinder v. 2.1 analysis of the cytb alignment containing 167 specimens of Etheostoma kennicotti and two E. flabellare outgroups identifies an optimal scheme where each codon position is treated as a separate data partition. The optimal molecular evolutionary models are K80+I+G for the first codon position, HKY+I for the second codon position, and GTR+G for the third codon positions. The Bayesian inferred cytb gene tree, which is a 50% majority-rule consensus tree calculated from the set of 5,000 posterior trees, is shown in Figure 2A. Seven inclusive monophyletic groups within Etheostoma kennicotti (s.l.) are resolved in the phylogeny and are supported with high (>0.94) Bayesian posterior probabilities. A clade comprising specimens sampled from the Green River system is the sister lineage of all other populations of E. kennicotti (s.l.) and E. cumberlandicum (Fig. 2). Within the Green River clade there is a split between the geographically disjunct Pond River and upper Green River (Figs. 1, 2A).

Fig. 2

Phylogenies of Etheostoma kennicotti sensu lato and Etheostoma cumberlandicum inferred from the mitochondrial cytb gene. (A) Bayesian inferred cytb gene tree. Bayesian posterior probabilities are indicated with filled circles. (B) Maximum clade credibility chronogram of Etheostoma kennicotti sensu lato and Etheostoma cumberlandicum. Bars indicate 95% posterior density of age estimates at nodes in the phylogeny.

There are three lineages from the Tennessee River system resolved in the cytb gene tree: a clade comprising specimens sampled from the upper Tennessee River system above Sale Creek in Rhea County, Tennessee and two lineages that comprise populations from the lower Tennessee River system (Fig. 2A). One of the lower Tennessee River clades comprises populations sampled from the Paint Rock, Flint River, and Sequatchie River systems and the Robinson Creek and Chambers Creek systems (Fig. 2A). Individual specimens collected from the same locations in Horse Creek, Hardin Co., Tennessee and Aldridge Creek, Madison Co., Alabama contain mtDNA haplotypes that resolve in both of the lower Tennessee River clades (Fig. 2A).

Populations of Etheostoma kennicotti (s.s.) from the Clarks River in Tennessee and Kentucky and the lower Ohio River system in Illinois and Kentucky resolve as reciprocally monophyletic sister groups in the cytb gene tree (Fig. 2A). The Ohio–Clarks lineage, which includes the type locality of Etheostoma kennicotti, is the sister lineage of a clade that encompasses all sampled populations from the Cumberland River system (Fig. 2A). The Cumberland lineage contains two clades: populations sampled from the Laurel River system and populations sampled from the Big South Fork and tributaries of the Cumberland River above Cumberland Falls that correspond to Etheostoma cumberlandicum (Jordan and Swain, 1883; Page and Smith, 1976; Etnier and Starnes, 1993: 500).

The relaxed molecular clock analysis in BEAST using the Yule pure birth branching model resulted in a posterior phylogeny of Etheostoma kennicotti (s.l.) and E. cumberlandicum that is congruent with phylogeny inferred using MrBayes (Fig. 2B). The time tree of E. kennicotti (s.l.) depicts a history of diversification that initiates in the Late Miocene and extends through the Pleistocene (Fig. 2B). The estimated age of the most recent common ancestor (MRCA) of all the lineages currently delimited as Etheostoma kennicotti and E. cumberlandicum is 6.05 million years ago (mya; 95% HPD: 4.53, 7.86 mya), which is the split between E. cf. kennicotti from the Green River and all other lineages of E. kennicotti and E. cumberlandicum (Fig. 2B). The age of the MRCA of the upper Tennessee and lower Tennessee River clades is 4.38 mya (95% HPD: 3.01, 5.83 mya). The age of the MRCA of E. kennicotti (s.s.) in the lower Ohio and Clarks River system and the two lineages in the Cumberland is 2.73 mya (95% HPD: 1.84, 3.76 mya). The timing of the split between E. cf. kennicotti in the Laurel River system and E. cumberlandicum is 1.27 mya (95% HPD: 0.69, 1.89 mya).

Twenty-five microsatellite loci were genotyped for 164 individuals of Etheostoma kennicotti (s.l.). The resulting data-set contains only 5% missing genotypes. For the six inclusive clades in the cytb gene tree (Fig. 2A), the mean number of private alleles per locus ranges from 3.2 to 5.6 (Table 3). The lower Tennessee River clade contains nearly double the number of private alleles and private allelic richness as any other population (Table 3), and 32% of the microsatellite loci (48/150) are not in HWE in at least one population (Supplemental Table 1; see Data Accessibility). None of the loci are consistently out of HWE in all populations, and one population (the lower Tennessee River) accounts for a disproportionate number (21/48) of departures from HWE. This lower Tennessee population has much lower observed (0.39) than expected (0.75) heterozygosity (Table 3).

Table 3

Microsatellite population summary statistics for Etheostoma kennicotti and E. cumberlandicum: sample size (n), proportion of polymorphic loci (P), mean number of alleles per locus (A), number of private alleles (nPA), private allelic richness (pAR), expected (H_e) and observed (H_o) heterozygosity.

The average pairwise F_ST value across the six inclusive clades of Etheostoma kennicotti (s.l.) and E. cumberlandicum in the cytb gene tree is 0.323 (Table 4). The lower versus upper Tennessee populations exhibit much lower differentiation (F_ST = 0.145) than all other population pairs. In contrast, E. cf. kennicotti from the Laurel River and E. kennicotti (s.s.) are the most differentiated (F_ST = 0.438), though several other lineage contrasts exhibit similarly high F_ST values (Table 4).

Table 4

Pairwise Weir and Cockerham (1984) F_ST estimated from 25 microsatellite loci for populations of Etheostoma kennicotti and E. cumberlandicum.

Hierarchical STRUCTURE analyses identify 18 clusters (K) as optimal to describe the genetic variation within Etheostoma kennicotti (s.l.) and E. cumberlandicum (Fig. 3A). Early stages of the hierarchical analyses identify genetic clusters consistent with the cytb phylogeny. For example, specimens of E. kennicotti (s.s.) from the lower Ohio and Clarks Rivers are identified as a distinct genetic cluster early in the hierarchical analyses (Fig. 3A). In the hierarchical analyses, ten of the 18 genetic clusters are localized in the Tennessee River basin, with small tributaries containing genetically distinct populations (Fig. 3A). This pattern is even stronger when analyzing K = 18 using all samples; 12 genetic clusters are restricted to the Tennessee River. Many genetic clusters in the Tennessee River represent individual sampling localities. The hierarchical analyses identify fine scale population structure within E. kennicotti (s.s.) in tributaries of the lower Ohio River with three genetic clusters corresponding to each sampled locality (Fig. 3A). In contrast, all specimens of E. cf. kennicotti from the Green River system from disjunct localities in the Pond River system and the upper Green River form a single genetic cluster (Figs. 1, 3A). Likewise, specimens of E. cumberlandicum form a single genetic cluster (Fig. 3A). Specimens of E. cf. kennicotti from the Laurel River are separated into two genetic clusters, with one cluster containing specimens sampled from Craig Creek and the other cluster containing all other sampled populations in the Laurel River system (Fig. 3A). The ancestry assignments for the individual specimens in the hierarchical STRUCTURE analysis are presented in Supplementary Figure 1 (see Data Accessibility).

Fig. 3

Population structure inferred from 25 microsatellite loci. (A) Hierarchical STRUCTURE analysis identifies 18 genetic clusters (K). Optimal K values for each round of hierarchical clustering are shown. (B) STRUCTURE results for K = 6. Pie charts represent average ancestry coefficients for each sampling locality. The phylogeny is redrawn from Figure 2A.

To determine whether genetic clustering is concordant with the lineages resolved in the cytb gene tree (Fig. 2A), we also visualized STRUCTURE results for K = 6. Four of the six genetic clusters match clades in the cytb phylogeny: Etheostoma kennicotti (s.s.) from the Clarks and Ohio Rivers, E. cf. kennicotti from the Green River, E. cf. kennicotti from the Laurel River, and E. cumberlandicum (Fig. 3B). There is little admixture among these genetic clusters except for one individual of E. cf. kennicotti from the Laurel River that exhibits appreciable genomic ancestry shared with E. cumberlandicum (Fig. 3B). Two genetic clusters are identified in the Tennessee River system. However, these clusters do not match the cytb gene tree (Figs. 2A, 3B).

Morphological analyses.—Scale row and fin elements counts are presented for 263 specimens of Etheostoma cumberlandicum and 440 specimens of E. kennicotti (s.l.) in Tables 5–10. The data for the number of scales above and below the lateral line, the number of scales around the caudal peduncle, and the number of anal-fin rays showed little variation and are not summarized in tables, but the meristic data for all specimens are available at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg). Plotting the first two principal components (PC) axes from the PC analysis shows broad overlap in the morphospace among the lineages resolved in the cytb gene tree and characterized as genetically distinct in the analysis of 25 microsatellite loci (Fig. 4A). The LDA of the meristic traits correctly classified 72% (180 of 250) of E. cumberlandicum and 80% (39 of 49) of E. cf. kennicotti from the Laurel River; however, it did not correctly identify a high proportion of the specimens of E. kennicotti (s.s.; 29%), E. cf. kennicotti from the Green River system (42%), and E. cf. kennicotti from the Tennessee River system (52%). Comparing the Mahalanobis distances of the PC values among the lineages of E. kennicotti (s.l.) shows the greatest disparity in the meristic traits is between the sister lineages of E. cumberlandicum and E. cf. kennicotti from the Laurel River (Fig. 4B), which were the two species with the highest proportions of correct specimen identification in the LDA.

Table 5

Counts of lateral line scales in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation.

Table 6

Counts of pored lateral line scales in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation.

Table 7

Counts of transverse scale rows in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation.

Table 8

Counts of first dorsal-fin spines in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation.

Table 9

Counts of second dorsal-fin rays in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation.

Table 10

Counts of left pectoral-fin rays in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation.

Fig. 4

Morphological disparity in Etheostoma kennicotti sensu lato and Etheostoma cumberlandicum. (A) Plot of first and second principal component scores of meristic traits in Etheostoma kennicotti sensu lato and Etheostoma cumberlandicum. (B) Mahalanobis distances of PC scores for contrasts of species in the Etheostoma kennicotti complex. In each comparison, the red and blue lineages on the cytb phylogeny are contrasted. A lineage that is a dashed branch is not included in the contrast. The comparison between the upper (U.) Tennessee and lower (L.) Tennessee River is considered an intraspecific contrast.

Despite broad overlap in the PC plot, moderate disparity in meristic traits among most of the lineages resolved in the cytb gene tree, and generally poor identification in LDA for three of the five species, there is notable variation in meristic and pigmentation traits that allows for morphological characterization of four of the five lineages of Etheostoma kennicotti (s.l.) and E. cumberlandicum. As noted in Page and Smith (1976), populations of E. kennicotti (s.s.) from tributaries of the lower Ohio River in Illinois and Kentucky and populations from the Clarks River system in Kentucky and Tennessee have fewer lateral line scales, fewer pored lateral line scales, and fewer transverse scales than other lineages of E. kennicotti (s.l.) and E. cumberlandicum (Tables 5–7). Etheostoma cf. kennicotti in the Laurel River has a modal count of eight spines in the first dorsal fin (Table 8). The population in the upper Tennessee River system does not exhibit a modal count of dorsal-fin spines as the same number of specimens have seven and eight spines, but E. cumberlandicum and all other lineages of E. kennicotti (s.l.) have a mode of seven spines in the first dorsal fin (Table 8). Etheostoma cf. kennicotti in the Green River system have modally 13 rays in the second dorsal fin (Table 9). Etheostoma cumberlandicum is distinct from all lineages of E. kennicotti (s.l.) in having modally 13 versus 12 rays in the left pectoral fin and a mode of seven bands of dark pigment in the caudal fin of adult male specimens (Tables 10, 11). The fewer caudal bands in E. cumberlandicum are apparent in comparisons with other lineages of E. kennicotti (s.l.), including its sister species E. cf. kennicotti from the Laurel River system (Fig. 5). There are no apparent meristic traits that allow a morphological diagnosis of E. cf. kennicotti from the Tennessee River system, and the Mahalanobis distance between the upper and lower Tennessee River lineages was the lowest contrast among all lineages of E. kennicotti (s.l.) and E. cumberlandicum (Fig. 4B).

Table 11

Counts of caudal bands in Etheostoma cumberlandicum and E. kennicotti. Abbreviations: n, number of specimens; SD, standard deviation. Only specimens 44.0 mm and greather in standard length included. Data included from Page and Smith (1976: table 7).

Fig. 5

Photographs of live specimens of Etheostoma cumberlandicum, E. cf. kennicotti Laurel River, and E. kennicotti. (A) Etheostoma cumberlandicum, YPM ICH 028205, 62 mm SL, male, Wolf Creek, Whitley Co., Kentucky, USA, 28 April 2015. (B) Etheostoma cumberlandicum, YPM ICH 028204, 46 mm SL, female, Little Wolf Creek, Whitley Co., Kentucky, USA, 25 March 2015. (C) Etheostoma cf. kennicotti, Laurel River species, YPM ICH 028202, 63 mm SL, male, Lick Fork, Laurel Co., Kentucky, USA, 14 May 2015. (D) Etheostoma kennicotti, UF 167303, 61 mm SL, Bay Creek, Pope Co., Illinois, USA, 18 April 2007.

DISCUSSION

Phylogenetic relationships and divergence times inferred from mitochondrial DNA sequences, population structure at 25 microsatellite loci, and divergence of meristic and pigmentation traits support the conclusion that Etheostoma kennicotti (s.l.) is composed of multiple species and E. cumberlandicum warrants recognition as a distinct species. The mitochondrial cytb gene tree resolves six major lineages which began diversifying approximately six million years ago, with the most recent speciation involving E. cf. kennicotti from the Laurel River system and E. cumberlandicum ∼1.3 million years ago (Fig. 2B). With the exception of samples from the Tennessee River system, population structure inferred from microsatellite loci is consistent with the mtDNA phylogenetic structure (Fig. 3). We suggest there are at least five species in the Etheostoma kennicotti complex: E. kennicotti (s.s.), E. cumberlandicum, E. cf. kennicotti distributed in the Tennessee River system, E. cf. kennicotti distributed in the Green River system, and E. cf. kennicotti distributed in the Laurel River system (Fig. 1).

The cytb gene tree resolves three deep branching lineages in the clade that comprises E. cf. kennicotti from the Tennessee River system (Fig. 2A). The clades resolved in the Tennessee River lineage correspond to a split between populations from tributaries of the upper Tennessee River and those from the lower portions of the Tennessee River (Fig. 2A). It is long accepted that phylogeographic structure can evolve without barriers to gene flow in nonrecombining genetic units such as mtDNA (Irwin, 2002). The splits in the cytb gene tree likely do not reflect geographic isolation among the populations in the lower Tennessee River system as haplotypes sampled from two populations, Horse Creek and Aldridge Creek, resolve in both of the lower Tennessee River system clades. A very similar pattern where there is little associating between geography and distribution of haplotypes in the Tennessee River system was observed in Etheostoma simoterum (Harrington and Near, 2012).

Populations of Etheostoma cf. kennicotti from the Tennessee River system contain the highest diversity of microsatellite alleles (Table 3) and the highest number of genetic clusters (Fig. 3A). However, genetic structure does not strongly match geography, and there are signals of admixture throughout the lower Tennessee River (Fig. 3A). Pairwise F_ST between the lower and upper Tennessee River is the lowest of any comparison among the six lineages resolved in the cytb phylogeny (Table 4). These patterns suggest the Tennessee River contains disjunct populations of E. cf. kennicotti connected by gene flow. The Tennessee River is well suited to produce such mosaic genetic structure, with many small but geologically stable tributary systems separated by larger, shifting river courses that may facilitate a mosaic of localized isolation with periods of gene flow among populations in the Tennessee River system. Despite high levels of genetic variation, we recognize all populations of E. cf. kennicotti in the Tennessee River, sans the Clarks River system, as a single undescribed species.

While the phylogenetic relationships and substantial genetic diversity strongly suggest that Etheostoma kennicotti (s.l.) is composed of five distinct species, meristic and pigmentation traits traditionally utilized to discover and delimit species of ray-finned fishes allow for diagnosis of only four of these lineages. The nominal E. kennicotti differs from all other species in the clade by having a lower average number of lateral line scales (Table 5). The species E. cf. kennicotti from the Green River system is distinguished by a larger number of transverse scale rows and modally 13 versus 12 rays in the second dorsal fin (Table 9). The species E. cf. kennicotti endemic to the Laurel River differs from all other lineages of E. kennicotti and E. cumberlandicum by having modally eight versus seven dorsal-fin spines (Table 8). Etheostoma cumberlandicum is diagnosed with 13 versus 12 pectoral-fin rays and modally seven versus nine or ten caudal bands (Tables 10, 11).

Page and Smith (1976) interpreted the lower scale counts in Etheostoma kennicotti (s.s.), the intermediate scale counts of E. cf. kennicotti from the Tennessee River system, and the higher scale counts in E. cf. kennicotti in the Laurel River and E. cumberlandicum as evidence for a pattern of east–west clinal variation in E. kennicotti (s.l.). This hypothesized cline was used to justify treating E. cumberlandicum as a synonym of E. kennicotti (s.l.). While the morphological variation in the scale counts detected in our study is similar to that presented by Page and Smith (1976), phylogenetic relationships inferred in the cytb gene tree do not support that E. cf. kennicotti from the Tennessee River system is “intermediate.” Instead, the Tennessee River lineage is sister to a clade containing E. kennicotti (s.s.), E. cf. kennicotti from the Laurel River system, and E. cumberlandicum, which are the lineages that exhibit the lowest and highest scale counts (Fig. 2A, Tables 5, 7).

The phylogenetic relationships and timing of diversification of species in the Etheostoma kennicotti complex are congruent with the paleogeographic history of rivers of the Eastern Highlands of North America. The pre-Pleistocene configurations of the present-day Green, Cumberland, and Tennessee River systems are sufficient to explain the origin of the lineages endemic to these systems (Mayden, 1988). The pre-glacial Green River was a tributary of the Old Ohio River that drained to the Old Mississippi River and was independent of the Old Cumberland and Old Tennessee Rivers (Burr and Page, 1986; Burr and Warren, 1986; Starnes and Etnier, 1986; Mayden, 1988). Etheostoma cf. kennicotti in the Green River system diverged from all other lineages of E. kennicotti in the Late Miocene (Fig. 2A, B), consistent with the fractured nature of the pre-glacial Eastern Highland River systems and the long isolation of the Green River system (Mayden, 1985; Mayden and Matson, 1992). In addition to E. cf. kennicotti, there are five other species of freshwater fishes endemic to the Green River system: E. barrenense Burr and Page, E. rafinesquei Burr and Page, E. barbouri Kuehne and Small, Nothonotus bellus (Zorach), and Thoburnia atripinnis (Bailey). The divergence times of these Green River endemics range from ∼15 mya for T. atripinnis and a clade comprising all three species of Hypentelium that is widespread throughout freshwater habitats in eastern North America (Bagley et al., 2018), 12.7 mya for the split between the clade containing the Green–Barren endemics E. rafinesquei and E. barrenense and the three species of the Etheostoma simoterum complex (E. simoterum, E. atripinne, and E. planasaxatile) distributed in the Tennessee, Cumberland, and Duck River systems (Near et al., 2011), 9.8 mya for the MRCA of E. barbouri and the Etheostoma basilare complex that is endemic to the Caney Fork River system (Hollingsworth and Near, 2009; Near et al., 2011), to 1.6 mya for the MRCA of the Green River endemic N. bellus and N. camurus that is widespread through the Ohio, Cumberland, and upper Tennessee River systems (Near et al., 2011). The relationships and divergence time of E. cf. kennicotti contribute to a set of observations that indicate the endemism of fishes in the Green River system is not the result of a single shared paleogeographic event or biogeographic process.

The elevation of Etheostoma cumberlandicum out of synonymy with E. kennicotti highlights patterns of endemism of freshwater fishes in the Cumberland River system above Cumberland Falls, which is a barrier to upstream dispersal of aquatic organisms. The darter E. susanae (Jordan and Swain) is endemic to the upper Cumberland upstream and in the immediate vicinity of the Cumberland Falls (Jordan and Swain, 1883; Starnes and Starnes, 1979; O'Bara, 1991; Strange, 1998). Etheostoma cumberlandicum and E. sagitta (Jordan and Swain) are near endemics to the upper Cumberland above the falls, but both species are distributed in the Roaring Paunch Creek system, a tributary of the Big South Fork that empties into the Cumberland River below the Cumberland Falls (Jordan and Swain, 1883; Kuehne and Bailey, 1961; Burr and Warren, 1986). Headwater stream capture between the upper Kentucky River system and the upper Cumberland is invoked as a mechanism of allopatric speciation leading to the origin of E. susanae and E. sagitta (Kuehne and Bailey, 1961; Starnes and Starnes, 1979; Strange, 1998); however, the origin of Etheostoma cumberlandicum involves allopatric speciation within the Cumberland River system. The Laurel River endemic Etheostoma cf. kennicotti and E. cumberlandicum are sister species with an estimated divergence time of 1.3 mya (Fig. 2A, B). This implies that the Cumberland Falls was an important geographic isolating barrier in the diversification of these two species (Fig. 1). The sister lineages of both E. susanae and E. sagitta are endemic to the upper Kentucky River system (Kuehne and Bailey, 1961; Strange, 1998; Heckman et al., 2009). As determined for Green River endemics, it appears endemism of freshwater fishes in the upper Cumberland River system is not the result of a single event involving vicariance between two specific river basins, but several events that involved allopatric processes between different areas adjacent to the upper Cumberland River system.

Conclusions.—Phylogenetic analysis of mtDNA sequence data and assessment of genetic clusters using microsatellite loci identify five distinct lineages in the Etheostoma kennicotti complex that we treat as allopatrically distributed species. Combining inferences of phylogenetic resolution, genetic clustering, and morphological distinctiveness, we elevate Etheostoma cumberlandicum out of synonymy with E. kennicotti. The delimitation and description of the three new species distributed in the Laurel River system, the Tennessee River system, and the Green River system (Fig. 1) require more data on male nuptial pigmentation and testing the phylogenetic resolution of the cytb gene tree using a method such as ddRADseq to collect nuclear DNA sequences from tens of thousands of loci (e.g., Near et al., 2021). The results of the phylogenetic and morphological study of E. kennicotti are part of a larger story where the applications of molecular and phenotypic data to species of darters with widespread distributions in the Eastern Highlands of North America results in the discovery of additional biodiversity masquerading as a single species. Ichthyologists studying the North American freshwater fish fauna are entering the last stages of describing all species of darters. As reflected by research of the past decade, this last phase of biodiversity discovery in darters will likely combine traditional morphological data with genetic and genomic analyses as has been deployed over the past decade (Keck and Near, 2013; Robison et al., 2014; Near and Thomas, 2015; Powers et al., 2015; Kozal et al., 2017; Near et al., 2017, 2021; Sterling and Warren, 2020).

DATA ACCESSIBILITY

Supplemental material is available at  https://www.ichthyologyandherpetology.org/i2021053 and at the Dryad Digital Repository ( https://doi.org/10.5061/dryad.47d7wm3jg). Unless an alternative copyright or statement noting that a figure is reprinted from a previous source is noted in a figure caption, the published images and illustrations in this article are licensed by the American Society of Ichthyologists and Herpetologists for use if the use includes a citation to the original source (American Society of Ichthyologists and Herpetologists, the DOI of the Ichthyology & Herpetology article, and any individual image credits listed in the figure caption) in accordance with the Creative Commons Attribution CC BY License.