Comparative Phylogeography, Historical Demography, and Population Genetics of Three Common Coastal Fauna in Spartina Marshes of the Northwestern Gulf of Mexico

Abstract

Coastal wetlands worldwide are experiencing high rates of loss and degradation that may lead to a reduction in diversity in faunal populations. Since salt marsh habitats are subject to a multitude of stressors, evaluations of the genetic diversity, connectivity, and potential resilience of faunal communities within salt marsh habitats are relevant. This study characterizes mitochondrial DNA (mtDNA) diversity for three common faunal residents of salt marshes along the northern Gulf of Mexico. Gulf Killifish (Fundulus grandis) samples were characterized for 1077 bp of the concatenated nucleotide sequence corresponding to the Control Region and Nitrogen Dehydrogenase, Subunits 2 and 5. Daggerblade grass shrimp (Palaemon pugio) samples were characterized using 466 bp of 16sRNA sequence, and phloem-feeding planthoppers (Prokelisia marginata) were characterized using 372 bp of Cytochrome c Oxidase Subunit I (COI) sequence. For F. grandis, our data revealed high levels of haplotypic diversity, evidence of isolation by distance (IBD), and regional population structuring associated with the distribution of two distinct phylogroups and distinct historical demography signatures. P. pugio and P. marginata displayed low levels of haplotypic diversity and evidence of population structure, but both appear to contain only snapshots of the total potential diversity for these species in the Gulf of Mexico. Greater resolution of the patterns of historical demography of Gulf Killifish may be obtained in future studies by including localities from Florida and Mexico. For both P. pugio and planthoppers, future studies would benefit from the characterization of genetic markers with a higher degree of polymorphism. We conclude that despite these three species inhabiting the same habitats along the same stretch of coast, each is subject to a different combination of evolutionary forces, and this study was able to reconstruct differences in how the genetic variation in each of these species emerged, and how it is maintained.

1. Introduction

Coastal wetlands function as essential wildlife habitats and provide vital ecosystem services to nearby coastal communities [1]. Despite their importance, coastal wetlands have decreased in area throughout much of the continental United States [2,3]. In a comprehensive review over a five-year period from 2004 to 2009, Dahl and Stedman [4] reported a loss of over 360 thousand acres of coastal wetland area in the conterminous U.S., or roughly 72 thousand acres per year. Such high rates of habitat fragmentation and loss can be expected to have enormous negative impacts on biodiversity [5,6], and the concomitant loss of faunal populations is expected to adversely affect the ecosystem services salt marshes provide [7,8,9].

Habitat fragmentation reduces population connectivity and genetic outcrossing, and populations with small effective population size (N_e) are prone to genetic erosion due to genetic drift. In turn, such losses potentially reduce fitness and adaptive responses to stressful conditions, consequently increasing the risk of extinction [10,11,12]. Due to the hierarchical nature of diversity, genetic diversity can have important ecological consequences at the population, community, and ecosystem levels, and in some cases, the effects are comparable in magnitude to the effects of species diversity [13]. For example, genetic studies on restored seagrass beds found that increased genetic diversity of dominant flora enhanced ecosystem services [14]. Since salt marshes and their resident fauna are subject to a multitude of stressors, including drought and dredging [15], nutrient loading [9,16], and tidal surges and fluctuations [17], it is germane to investigate the genetic diversity, connectivity, and potential resilience of faunal communities within salt marsh habitats, particularly when facing the challenges associated with climate change [15,18].

This study focuses on characterizing genetic variability in Spartina marsh fauna. To ensure a representative coverage of faunal salt marsh communities, we selected two aquatic species, a fish and a crustacean, and an aerial insect, all known to be intimately associated with this habitat [19,20,21,22,23]. Each of the chosen species displays distinct life history characteristics that uniquely illustrate patterns of gene-flow potential relative to their habitat (aquatic versus aerial), and modes of dispersal and reproduction, with the expectation that measurable differences in genetic variability exist both among and within species. Comparative studies of species with similar geographic distributions can help reveal shared evolutionary events over both geography and time [24,25,26]. We chose Gulf Killifish, Fundulus grandis (Cyprinodontiformes), and daggerblade grass shrimp, Palaemon pugio (Decapoda), previously Palaemonetes pugio, for the two aquatic species because they are among the most abundant nekton in salt marsh habitats along the northern Gulf of Mexico (hereafter, Gulf) coast, and represent vital links in coastal marsh food webs [27,28,29]. For the aerial example, we selected the phloem-feeding planthopper, Prokelisia marginata (Hemiptera). Prokelisia spp. are host-specific sap-feeders that, together with leafhoppers and mirid bugs, may reach combined densities that account for more than 90% of herbivore biomass in some habitats [22,30,31].

The three species characterized in this study encompass a variety of life history traits and dispersal potentials. Fundulus grandis have relatively low fecundity and are benthic spawners that attach their eggs to substrate within the marsh [32,33,34,35], giving them low dispersal potential. Mark-recapture studies along the Gulf Coast have shown that F. grandis exhibits high site fidelity, with individuals traveling 100 m or less between connected marsh habitats [36]. Their reproductive characteristics and limited movement as adults combine to make F. grandis an excellent indicator species in terms of genetic and physiological responses to exposure to a variety of industrial toxins, and they have been used extensively as such [37,38,39,40,41,42,43].

P. pugio, by contrast, reach sexual maturity quickly and have a short generation time, high fecundity, and high dispersal potential in their nauplii stages [27,44,45,46]; however, they appear to exhibit high levels of site fidelity as adults [47]. They play a vital role in the food webs of coastal marsh habitats by breaking down organic materials and making those nutrients available at other trophic levels [44,45,48,49], and for these reasons, they are frequently used as indicator species in studies of coastal habitat quality and toxicity response [50,51,52].

Phloem-feeding delphacid planthoppers (hereafter, planthoppers), such as P. marginata, are closely associated with plant cover of native marshes. Previous studies have shown an association between abundance of Prokelisia spp. and the ecosystem services of Spartina-dominated salt marshes [7,8]. Further, their mode of dispersal is relevant to our understanding of population connectivity between marshes. Planthoppers may present in a wingless, more fecund brachypterous form or a winged, less fecund macropterous form with higher dispersal capabilities, depending upon environmental conditions and stressors [22,53]. Along the Gulf coast, the wingless form with limited dispersal capabilities is the most common [54].

Previous genetic studies of the F. grandis in the Gulf indicate limited dispersal, primarily between adjacent estuaries, and significant levels of isolation by distance using both single nucleotide polymorphisms (SNPs) and microsatellites [55,56]. Population genetics studies of P. pugio are few, but include an allozyme study that compared nine collections on and around Galveston Island that reported lower levels of variation in isolated, recently formed ponds compared to those from larger populations occupying older bodies of water open to migration [57,58]. By contrast, high levels of gene flow over a wide geographic range in marshes along the S. Atlantic coast of the U.S. (hereafter S. Atlantic) were invoked to explain the reduced genetic partitioning among distant populations of P. pugio using mitochondrial DNA (mtDNA) 16sRNA single strand conformation polymorphism (SSCP) [59]. In P. marginata, mtDNA Cytochrome Oxidase Subunit I (COI) data revealed a pronounced population structure at very large geographic scales (e.g., Gulf and S. Atlantic versus the N. Atlantic U.S. coast, hereafter N. Atlantic), and the signature of the corresponding clades was used to identify the source of putatively introduced populations [60].

Mitochondrial data for many individuals, particularly across a group of co-occurring species, can provide a baseline from which to generate questions for deeper investigation [25]. This study is part of a larger study aiming to characterize genetic variation in both mtDNA and nuclear DNA, seeking to add to the current knowledge of these three important salt marsh residents. Here, we analyze mtDNA sequence data to investigate the genetic diversity and connectivity of populations from estuaries along the north and west Gulf coast for the three species, which differ in life history patterns and dispersal potential. Dispersal potential was expected to be highest for the insect P. marginata because of its aerial capabilities, and can be hypothesized to translate into levels of gene flow that are higher than that of the two aquatic representatives. The crustacean P. pugio, due to its high fecundity and planktonic larval dispersal, should rank second highest, while the fish F. grandis, due to its low fecundity and low dispersal potential, was expected to rank last and consequently display the highest level of genetic structuring. For this interpretation to be sound, it is necessary to invoke comparative phylogeography and historical demography to examine the extent to which vicariance, population expansion, and variance in reproductive success account for their respective regional phylogenies and local genetic structuring.

2. Materials and Methods

F. grandis, P. pugio, and P. marginata were captured in the spring and summer months of 2014–2017 in Spartina alterniflora salt marshes along the northwestern Gulf of Mexico coastlines of Texas, Louisiana, and Mississippi (Figure 1). All specimens were immediately preserved in the field using 95% ethanol, and were transferred to 70% ethanol for long-term preservation within 24 h of collection. F. grandis were sampled using minnow traps baited with dog or cat food kibbles, placed in shallow (<0.5 m) water in a low marsh habitat, and allowed to soak for up to 12 h. Specimens were visually identified in the field and humanely sacrificed via immersion in MS-222 as per U.S. federal policies on the use of laboratory animals as subjects (AUP# 2014-0111 and 2017-0105). Palaemon spp. were collected using a dip net along the edge of the marsh habitat. Morphological identification to species level was carried out on well-preserved, intact specimens using a dichotomous key for this genus [44,61], and verified via the mtDNA sequences generated in this study. Arthropods, including planthopper specimens, were obtained using professional insect 15” muslin sweep nets (www.gemplers.com (accessed on 1 May 2014)). Preserved samples were examined under a dissecting microscope at 10× magnification, and specimens with the characteristics of Prokelisia spp. were separated from other arthropods. P. marginata and P. dolus are extremely similar morphologically, and distinguishing these two species is very difficult in males and impossible in females [19,62]. Therefore, to avoid inflating estimates of diversity via the inclusion of cryptic species [63,64], Prokelisia spp. were identified to species level by characterizing their mtDNA with a high-resolution melting analysis (HRMA) assay developed specifically to distinguish these two species [65]. Following identification, representative samples of each species from each location were sequenced and analyzed.

DNA from F. grandis and P. pugio was isolated from the axial muscle using a Zymo Quick-DNA Universal Kit, following the manufacturer’s instructions for tissue (Zymo Research, Irvine, CA, USA). Because of their small size (<5 mm), DNA isolation from P. marginata required individually grinding the entire specimen with a sterile disposable pestle. DNA was initially isolated using the Zymo Quick-DNA Universal Kit, but yields were low. However, since higher DNA yields were obtained using a Qiagen Pure-gene extraction kit (Qiagen, Hilden, Germany) with the addition of 1 μL glycogen solution, as recommended by the manufacturer, this kit was used with the majority of P. marginata.

For F. grandis, three sets of primers were designed to target a segment of the mitochondrial Control Region (CR1), and the corresponding segments from Nitrogen Dehydrogenase Subunits 2 and 5 (ND2 and ND5) as their mutation rates are sufficient for investigating population structure [66]. For P. pugio and P. marginata, we targeted segments of 16sRNA (16S) and Cytochrome C Oxidase Subunit I (COI), respectively, as previous studies have used these markers successfully in population studies of these species [59,60]. All primer sets used in this study are either primers from other studies modified to match our target species, or were designed for this study using the Primer 3 software [67,68] within Geneious v.9.1.8 (Biomatters Ltd., Aukland, New Zealand) as the most optimal among the potential primer pairs capable of amplifying an in silico fragments 400–600 bp in length for the targeted regions. The mitogenomes of F. grandis (Accession # FJ445396) and P. pugio (Accession # EU868697) were, respectively used as a reference when designing primers for those species.

Table 1. PCR primer summary for each species and the mtDNA marker sequenced in this study along the northern Gulf of Mexico. Nucleotides that have been changed from the originally published primer sequence are denoted in bold, and nucleotides that have been inserted into the originally published sequence are underlined.

Locus	Primer Names	Primer Sequence	Fragment Size	Annealing Temp. (°C)
F. grandis CR-1 1	L15998-FG CSBD-H-FG	5′ CGC CCC TAG CTC CCA AAG CTA 3′ 5′ AAT AGG AAC CAA ATG CCA G 3′	400 bp	50
F. grandis ND-2 2	L4173ND2-FG H4634ND2-FG	5′ CAT CAT CCC CGA GCC GTT GA 3′ 5′ GGA AGG TTA AGG ATG GGA AG 3′	421 bp	50
F. grandis ND-5 2	L12137ND5-FG H12717ND5-FG	5′ GCA GAA ACG GTA GTG TCC AC 3′ 5′ GTA CTT GAA TGC AGT AGG GC 3′	540 bp	50
P. pugio 16sRNA 3	L-16sRNA-PP H-16sRNA-PP	5′ TGC CCT GTT TAT CAA AAA CAT 3′ 5′ AGA TAG AAA CCC AAC CTG G 3′	470 bp	50
P. marginata COI 4	C1-J-1751 C1-N-2197	5′ GGA TCA CCT GAT ATA GCA TTC CC 3′ 5′ CCC GGT AAA ATT AAA ATA TAA ACT TC 3′	400 bp	57

¹ Modified from Alvarado-Bremer and others [69]; ² Generated for this study; ³ Modified from Crandall and Fitzpatrick [70]; ⁴ Directly from Simon and others [71].

provides a summary of the primers used in this study, and when applicable the source of primer sets that were modified to match the species characterized in this study [69,70,71].

PCR was carried out separately for each locus in 12.5 μL reactions containing 1X EconoTaq Plus Green Master Mix (Lucigen, Middleton, WI, USA), 0.2 μM of each primer, and 10–20 ng of isolated DNA as template. Thermocycling was performed on an Eppendorf Mastercycler (Eppendorf, Hamburg, Germany) with an initial denaturing step at 94 °C for 2 min, followed by 35 cycles of denaturing 94 °C for 25 s, annealing at the corresponding temperature for each primer pair (See Table 1) for 30 s, and extensions at 72 °C for 90 s, before a final extension step at 72 °C for 3 min. Negative controls were included in all reactions. PCR products were then visualized for specificity and yield via electrophoresis on a 2% agarose gel pre-stained with ethidium bromide (EtBR). PCR products that produced a single band were diluted 1:10 in ddH₂O for post-PCR cleanup and sequenced in both directions, with reaction setups and thermocycling profiles as described in Cruscanti et al. [72]. Individual haplotypes derived from the sequences generated herein have been submitted to Genbank, and Accession numbers are provided in Supplementary Table S1.

Multiple sequence alignments were carried out in Geneious Pro v.9.1.8 (Biomatters Ltd., Aukland, New Zealand). Haplotype data files were generated in DnaSP v6.12.03 [73]. Arlequin v.3.5.2.2 [74] was used to estimate genetic diversity within sampling locations, and to calculate sequence diversity indices, pairwise F_ST, coancestry coefficients (Reynold’s Distance), and Slatkin’s Linearized F_ST [75,76,77,78]. Reynold’s distance was calculated in the event that genetic differentiation occurs only via genetic drift without mutations [75]. A Mantel test of isolation by distance (IBD) was calculated using Slatkin’s linearized F_ST correlated against pairwise distances between sampling locations [77,79]. p-values for pairwise comparisons of the six localities were corrected for multiple testing using Benjamini and Hochberg’s method [80], which corrects for significance by controlling the false discovery rate (FDR) and produces fewer false negatives than Bonferroni corrections [81]. Since the distribution of diversity statistics falls on an asymptotic curve rather than a normal curve, the Salicrú χ2 method [82] was used to test for pairwise significant differences in haplotypic diversity between sampling locations. POPart v.1.7 [83,84] was used to build median-joining networks (MJN) for F. grandis markers. The MJN was chosen for F. grandis sequences because of the large number of haplotypes separated by small genetic distances [83]. Due to the small number of haplotypes in both P. pugio and P. marginata, the respective relationship among lineages was reconstructed with minimum spanning networks (MSN) using POPart v.1.7. A representative of F. heteroclitus (Accession # KT869378) was used as the outgroup for the MJN of F. grandis. Representatives of P. vulgaris and P. dolus sequenced in this study were used as the outgroups for the MSNs of P. pugio and P. marginata, respectively. Principle component analyses (PCAs) as implemented in R v3.6.1 [85] were also used to identify structure in the distribution of mtDNA variation of F. grandis along the northwestern Gulf of Mexico coastlines of Texas, Louisiana, and Mississippi.

Spatial analysis of molecular variance (SAMOVA) was conducted to identify population structure [86]. SAMOVA is similar to AMOVA [87] except that it has the potential to identify genetic barriers between sampling groups without a priori constraints on the geographic composition of the groups. Since each species was sampled at six distinct locations, SAMOVAs were tested with two to five groups, and the grouping that produced the highest F_CT value was chosen as the best partitioning scheme.

To investigate patterns of historical demography and estimate female effective population sizes (N_ef), F. grandis samples were pooled into the most optimal hierarchical arrangement scheme based on the SAMOVA that was congruent with the phylogeographic association. For P. pugio and P. marginata, due to overall low levels of haplotypic diversity, it was necessary to pool all the samples within species to obtain a semblance of their demographic histories and to estimate N_ef. Estimates of N_ef were obtained as described in Roman and Palumbi [88]. Briefly, to determine the rate of divergence, each species was compared against the corresponding segments of the mitogenome of the closest relative in DnaSP v6.12.03 to obtain Tamura–Nei gamma-corrected distances (Da) between the species, and the timing of the corresponding speciation events from the literature. Accordingly, F. grandis populations were compared against seven representatives of the Atlantic sister species, F. heteroclitus (Accession # KT869378, FJ445398, FJ445399, FJ445401, FJ445402, FJ445403, NC_012312). P. pugio were compared against the mitogenome of five representatives of their sympatric sister species, P. vulgaris (Accession # JQ042300, JN674358, KP178999, KT959473, KT959519). P. marginata were compared against five representatives, each a different haplotype, of their sympatric sister species, P. dolus, that were sequenced by the authors for a previous study [65]. Generation time was assumed to be one year for F. grandis [33], two months for P. pugio [45], and 1.5 months for P. marginata [53]. For F. grandis, mutation rates were calculated based on the time since divergence from their Atlantic sister species based on Avise [26], who estimated that divergence times between the Gulf and Atlantic Ocean for coastal fish species range from 0.5–4.8 million years (mean = 1.3 million years). For P. pugio and P. marginata, comparisons with sympatric sister species prevent using geological events to estimate mutation rates. Therefore, mutation rates of 0.9%–1.1% per MY for P. pugio [89] and 2.7% per MY for P. marginata [90] were obtained from the literature for closely related taxa. Pairwise mismatch distributions [91,92], the D statistic by Tajima [93], the R2 statistic by Ramos-Onsins and Rozas [94], and the estimated mutational time, tau (τ) since population expansion [95] were generated to estimate patterns of historical demography.

3. Results

3.1. The Gulf Killifish, Fundulus grandis

Sequences for three mtDNA segments, namely CR1 (336 bp), ND2 (344 bp), and ND5 (397 bp), were obtained from 166 F. grandis specimens from six sampling locations. These sequences were concatenated into one single segment 1077 bp long, containing 176 segregating sites (

Table 2. Molecular indices for 1077 bp of the concatenated mtDNA nucleotide sequences of the Control Region, ND2, and ND5 genes of F. grandis (see Table 1) by sample location along the northern Gulf of Mexico. M, no. of haplotypes; h, haplotypic diversity; π, nucleotide diversity; SD, standard deviation; S, no. of segregating (polymorphic) sites; Ts, no. of transitions; Tv, no of transversions; I/D, no. of insertions and/or deletions.

Location	N	M	h (SD)	π (SD)	S	Ts	Tv	I/D
Corpus Christi	26	15	0.945 (0.024)	0.013 (0.007)	64	55	7	2
Matagorda	18	16	0.980 (0.028)	0.019 (0.010)	64	55	9	2
Reitan Marsh	29	24	0.988 (0.012)	0.019 (0.010)	77	71	5	2
Sportsman’s	24	23	0.956 (0.015)	0.020 (0.010)	80	70	10	3
Venice	36	18	0.886 (0.037)	0.004 (0.002)	37	30	8	0
Ocean Springs	33	26	0.983 (0.012)	0.007 (0.003)	54	43	13	1
All Samples	166	109	0.987 (0.003)	0.018 (0.009)	176	142	34	4

) that define 109 distinct haplotypes. For all loci, patterns of genetic variability within and among-localities were estimated. Nearly identical patterns of differentiation and diversity were obtained by analyzing each of these segments separately (Supplemental Tables S2–S8); therefore, only the results for the concatenated segment are reported below.

High levels of haplotypic diversity (h > 0.88) were found in all sampling locations, but differences in levels of genetic variation among some of the samples were observed (Table 2). Values of haplotypic diversity were significantly lower in Venice than all other locations except Corpus Christi. Within Galveston Bay, Sportsman’s Road was significantly less varied than the Reitan Marsh, which contained the highest overall value of haplotypic diversity (

Table 3. Values for pairwise comparisons of northern Gulf of Mexico F. grandis samples estimated from concatenated (1077 bp) sequences. Z-scores from Salicrú Χ² test for pairwise comparisons of haplotypic diversity are above the diagonal. Pairwise F_ST values are below the diagonal. Significant values are in bold, with significance at p < 0.05 denoted by *, and significance at p < 0.01 denoted by **.

	Corpus Christi	Matagorda	Reitan Marsh	Sportsman’s	Venice	Ocean Springs
Corpus Christi		−0.949	−1.603	−0.389	1.338	−1.416
Matagorda	0.18880 **		−0.263	0.756	2.026 *	−0.098
Reitan Marsh	0.13947 **	−0.01239		1.666 *	2.622 **	0.295
Sportsman’s	0.20315 **	−0.02873	−0.00191		1.753 *	−1.406
Venice	0.67054 **	0.42053 **	0.44186 **	0.37772 **		−2.494 **
Ocean Springs	0.63133 **	0.37057 **	0.40474 **	0.34059 **	0.03746 *

). Comparison of pairwise F_ST values for F. grandis identified Corpus Christi as significantly different from all other locations (p < 0.01). Similarly, Venice and Ocean Springs differ significantly from all other locations (p < 0.01) and from each other (p < 0.05; Table 3). Calculations of the Slatkin’s linearized F_ST and Reynold’s distance (

Table 4. Values for pairwise comparisons of F. grandis concatenated mtDNA sequences. The Slatkin’s linearized F_ST between the northern Gulf of Mexico samples is above the diagonal, and the Reynold’s distance is below the diagonal.

	Corpus Christi	Matagorda	Reitan Marsh	Sportsman’s	Venice	Ocean Springs
Corpus Christi		0.03221	0.03032	0.03539	0.08943	0.03664
Matagorda	0.03274		0.00638	0.01289	0.07185	0.01846
Reitan Marsh	0.03079	0.00640		0.00000	0.06613	0.01480
Sportsman’s	0.03602	0.01297	0.00000		0.06741	0.01578
Venice	0.09355	0.07449	0.06836	0.06973		0.03300
Ocean Springs	0.03732	0.01863	0.01491	0.01591	0.03355

) yielded similar relationships. A Mantel test is consistent with IBD (R² = 0.5385; p = 0.0129) among the F. grandis sampling locations characterized in this study (Figure 2).

SAMOVA results for F. grandis produced the highest value of interpopulation variance (F_CT) when samples were allocated into three regional groups consisting of (1) Corpus Christi; (2) Matagorda, Galveston (Reitan Marsh, and Sportsman’s Road); and (3) Venice and Ocean Springs (

Table 5. SAMOVA results for F. grandis concatenated mtDNA sequences. The highest Among Groups (F_CT) value was obtained with three groups, as follows: Population 1 (S.Tex): Corpus Christi; Population 2 (E.Tex): Matagorda, Reitan Marsh, and Sportsman’s; and Population 3 (CN.Gulf): Venice and Ocean Springs.

Source of Variation	d.f.	Sum of Squares	Variance Components	Percentage of Variation
Among Groups	2	485.23	4.5071 Va	39.45
Among Populations Within Groups	3	21.83	0.0135 Vb	0.12
Within Populations	163	1125.51	6.9050 Vc	60.43
Total	168	1632.57	11.42551
Fixation Indices			p-values (≥)
		FSC:	0.45064 ± 0.01396
		FST:	0.00000 ± 0.00000
		FCT:	0.01271 ± 0.00366

), which correspond to the regional sampling carried out in this study. Matagorda samples were originally expected to group with Corpus Christi samples, since the upper Laguna Madre and Matagorda Bay systems are contiguous; however, all analyses performed in this study indicate that Matagorda lineages are more closely related to the Galveston Bay samples, and will therefore be grouped together in ensuing analyses. Accordingly, hereafter, these populations will be referred to as S.Tex, E.Tex, and CN.Gulf, respectively.

An MJN (Figure 3) illustrates the relationships among F. grandis mtDNA haplotypes (n = 109) and identifies two major phylogroups separated from each other by at least 12 mutations, which display distinct phylogeographic associations. The first phylogroup (Phylogroup I) appears to be more closely related to F. heteroclitus. It is important to note that this placement of the root does not coincide with that obtained in a phylogenetic reconstruction of the same mtDNA lineages (Supplemental Figure S1). Phylogroup I includes, near its baseline, lineages primarily (87.1%) from E.Tex, although it includes two haplotypes from S.Tex and three haplotypes from CN.Gulf (two from Ocean Springs and one from Venice), each found only once. About half of the total number of E.Tex haplotypes cluster within this mostly western Phylogroup I subgroup. The second subgroup of Phylogroup I consists of lineages private to the CN.Gulf locations of Venice and Ocean Springs, and includes two haplotypes whose frequency is higher than any other F. grandis haplotype characterized in this study. From these two centroids, very closely related haplotypes, one or two mutational steps apart, respectively, emerge in a pattern concordant with star phylogenies. These two centroids reach the highest frequency in Venice, the location with the lowest value of haplotypic diversity, but are also present in Ocean Springs, which is also significantly less variable than any locations sampled in Texas (Table 2). By contrast, Phylogroup II contains lineages exclusively from S.Tex and E.Tex, with the other half of E.Tex lineages belonging to this second phylogroup. The majority (80.0%) of F. grandis from S.Tex (Corpus Christi) belong to Phylogroup II, and include eight divergent haplotypes located at the terminal branches of the MJN, each separated by >6 mutations from their nearest neighbor. However, two haplotypes from S.Tex were shared with E.Tex (Matagorda and Reitan Marsh), and a third S.Tex haplotype, repeated three times, clusters among other haplotypes from E.Tex. Lastly, intermediate to phylogroups I and II, there are two haplotypes from E.Tex and one from S.Tex, which collectively are separated from these groups by 12 and eight mutations, respectively. It should be noted that this intermediate group lies at the base of the phylogenetic tree, being more closely related to the outgroup F. heteroclitus (Figure S1).

An investigation of the frequency distribution of phylogroups by sampling location (Figure 4) clearly depicts the phylogeographic association shift from west to east, with the majority of S.Tex lineages belonging to Phylogroup II, all CN.Gulf lineages belonging to eastern subgroup of Phylogroup I, and with E.Tex lineages being split nearly evenly between the Phylogroups II and the western component of Phylogroup I. The first two axes of the PCA, which account for about 96.0% of the observed variance, also summarize the patterns of differentiation among F. grandis mtDNA lineages (Figure 5). Congruent with the MJN, PCA loadings along PC1 separate most haplotypes into lineages belonging to the ubiquitous Phylogroup I, characterized by strong positive loading (>1.0), and another cluster of lineages, exclusively from Texas, belonging to Phylogroup II with strong negative loadings (<−2.0). Negative loadings (≈−1.5) along PC1 also identify the three outliers that correspond to the lineages identified as intermediate in the MJN (see Figure 3). The loadings on PC2 largely separate the eastern component of Phylogroup I, with the majority of Venice and Ocean Springs lineages having positive scores (>0.0).

The historical demographic signature of F. grandis differed among the regions of the Gulf of Mexico surveyed. Both D and R2 statistics suggest population expansion (α < 0.05) in in the CN.Gulf population (

Table 6. Historical population demography parameters and N_ef for F. grandis populations along the northern Gulf of Mexico, based on SAMOVA results. Da, Tamura–Nei corrected genetic distance between the population and sister species, F. heteroclitus; τ, estimated mutational time since population expansion; D, Tajima’s D with probability value (p); R2, Ramos-Onsins and Rozas’s R2 with probability value (p); significant tests, in bold; T, time since divergence in millions of years used for mutation rate estimations.

Population	Da	τ	D (p)	R2 (p)	T	Nef
S.Tex	0.076	6.50	−0.483 (0.346)	0.109 (0.368)	0.50	50.48
					1.25	126.21
					4.8	484.65
E.Tex	0.072	8.62	−0.343 (0.424)	0.888 (0.417)	0.50	80.74
					1.25	201.84
					4.8	775.07
CN.Gulf	0.078	2.95	−2.058 (0.003)	0.035 (0.001)	0.50	18.89
					1.25	47.23
					4.8	181.35

). Mismatch distributions of pairwise differences d S.Tex and E.Tex yielded multiple peaks over wide range of pairwise differences (0–40). By contrast, the CN.Gulf distribution is unimodal, with pairwise differences ranging between 0–15 (Figure 6). The corresponding values of τ suggest that the populations E.Tex and S.Tex had long and stable histories, which are two to three times older than the CN.Gulf population. Estimates of N_ef ranged from 19 million to 775 million females, depending on the estimated time since divergence used in the calculation (Table 6).

Since the MJN identifies two phylogroups that have strong phylogeographic associations, the estimates of N_ef may be upward biased, see Alvarado Bremer et al. [24], as they would be the product of two separate evolutionary histories, as would specifically be the case for the E. Texas and S. Texas populations of F. grandis. Accordingly, mismatch distributions and historical demography statistics, and consequently N_ef, were recalculated separately for phylogroups I and II for E.Tex (

Table 7. Historical population demography parameters and estimates of female effective population size for F. grandis phylogroups depicted in the MJN. Da, Tamura–Nei corrected genetic distance between the population and sister species, F. heteroclitus; τ, estimated mutational time since population expansion; D, Tajima’s D with probability value (p); R2, Ramos-Onsins and Rozas’s R2 with probability value (p); T, time since divergence from sister species in millions of years used for mutation rate estimations, and long-term female effective population size (N_ef) in millions.

Grouping	Da	τ	D (p)	R2 (p)	T	Nef
Phylogroup I	0.076	5.63	−2.019 (0.003)	0.034 (0.005)	0.50	17.76
All Lineages					1.25	44.39
					4.8	170.43
Phylogroup I	0.074	6.34	−1.578 (0.043)	0.059 (0.011)	0.50	9.41
E.Tex Linages					1.25	23.51
					4.8	90.29
Phylogroup I	0.078	2.95	−2.058 (0.003)	0.035 (0.001)	0.50	18.89
CN.Gulf Lineages					1.25	47.23
					4.8	181.35
Phylogroup II	0.79	6.19	−1.106 (0.127)	0.069 (0.127)	0.50	26.10
All Lineages					1.25	65.25
					4.8	250.57
Phylogroup II	0.81	5.16	−1.312 (0.076)	0.071 (0.064)	0.50	16.14
E.Tex Lineages					1.25	40.34
					4.8	150.25

), whereas for S.Tex, the historical demography estimates were only recalculated for Phylogroup II, due to the paucity of Phylogroup I lineages. For the CN.Gulf, these estimates remained the same as the original calculations (Table 6), since all the lineages found in Louisiana and Mississippi belong to Phylogroup I (Figure 3 and Figure 4). The recalculated D and R2 tests for Phylogroup I lineages in E.Tex were significant, congruent with the population expansion recorded for this Phylogroup in the CN.Gulf, although the estimates of N_ef for this latter region are twice as large as E.Tex. The separate mismatch distributions of each phylogroup in E.Tex were considerably less variable and unimodal for Phylogroup I. However, these tests were not significant for Phylogroup II, both in E.Tex and S.Tex. Overall, the number of pairwise differences was reduced by half compared to the analysis with all lineages together (Table 7). The number of pairwise differences for S.Tex distribution by phylogroup translated into a reduction from 35 pairwise differences in the original analysis to 21 pairwise differences (Figure 7). Values of τ indicate that Phylogroup II has a slightly longer and stable history than Phylogroup I, and that the two subgroups of Phylogroup II underwent concurrent expansions. By contrast, comparisons of τ values within Phylogroup I indicate that the E.Tex lineages are twice as old as the CN.Gulf lineages. The new estimates of N_ef for E.Tex based on the sum of the two phylogroups are about four times smaller than when all lineages are included. For S.Tex, N_ef dropped by about 15% (see Table 6 and Table 7).

3.2. Daggerblade Grass Shrimp P. pugio

P. pugio specimens (n = 119) from six sampling locations were successfully sequenced for 466 bp of the mtDNA 16sRNA gene containing nine variable sites, which define thirteen haplotypes throughout the region studied. Except for Ocean Springs, the values of h obtained for most sampling locations were extremely low (

Table 8. Molecular indices for 466 bp of 16sRNA sequence of P. pugio by sample location along the northern Gulf of Mexico. M, no. of haplotypes; h, haplotypic diversity; π, nucleotide diversity; SD, standard deviation; S, no. of segregating (polymorphic) sites; Ts, no. of transitions; Tv, no of transversions; I/D, no. of insertions and/or deletions.

Location	N	M	h (SD)	π (SD)	S	Ts	Tv	I/D
Port Mansfield, TX	20	1	0.000 (0.000)	0.000 (0.000)	0	0	0	0
Matagorda, TX	20	4	0.284 (0.128)	0.001 (0.000)	3	1	1	0
Sportsman’s Rd, TX	20	1	0.000 (0.000)	0.000 (0.000)	0	0	0	0
Pt. Arthur, TX	20	2	0.100 (0.088)	0.001 (0.000)	1	0	0	0
Venice, LA	19	4	0.667 (0.086)	0.002 (0.002)	3	1	1	0
Ocean Springs, MS	20	6	0.811 (0.047)	0.003 (0.003)	5	2	2	0
All Samples	119	13	0.779 (0.021)	0.005 (0.003)	9	6	3	0

), with the two E.Tex localities (Pt. Mansfield and Sportsman’s Road) each displaying a single haplotype (h = 0) (Figure 8). The single Sportsman’s Road haplotype was private to that locality, but its two nearest neighbors came from Venice and Ocean Springs, whereas the Pt. Mansfield haplotype was shared with these two CN.Gulf localities. Further, Venice and Ocean Springs collectively encompass all the variability summarized in the MSN, extending from the closest relative to the outgroup, P. vulgaris, to the most divergent lineage from that outgroup. It is also worth noting that in those two CN.Gulf localities, no haplotype exceeds 50% in frequency. The two remaining localities, Matagorda and Port Arthur, both in E.Tex, share a common haplotype also found in Ocean Springs. The relationship of that common haplotype to other lineages in these two E. Texas localities resembles a star phylogeny, with daughter lineages one or two mutational steps apart from the centroid (Figure 8).

The levels of haplotypic diversity varied along the coastal range sampled, with the Salicrú χ2 test indicating significant differences in the levels of variability among P. pugio samples (

Table 9. Values for pairwise comparisons of northern Gulf of Mexico P. pugio samples estimated from 466 bp of 16sRNA. Z-scores from Salicrú Χ² test for pairwise comparisons of haplotypic diversity are above the diagonal. Pairwise F_ST values are below the diagonal. Significant values are in bold, with significance at p < 0.05 denoted by *, and significance at p < 0.01 denoted by **.

	Pt. Mansfield	Matagorda	Sportsman’s	Pt. Arthur	Venice	Ocean Springs
Pt. Mansfield		−2.219 *	0.000	−1.136	−7.756 **	−17.255 **
Matagorda	0.93285 **		2.219 *	1.185	−2.484 **	−3.865 **
Sportsman’s	1.000000 **	0.96026 **		−1.136	−7.756 **	−17.255 **
Pt. Arthur	0.98361 **	0.02105	0.99010 **		−4.608 **	−7.127 **
Venice	0.45823 **	0.81601 **	0.65339 **	0.86007 **		−1.469
Ocean Springs	0.46252 **	0.40897 **	0.73527 **	0.47953 **	0.41360 **

). Notably, the CN.Gulf localities (Ocean Springs and Venice) are more variable (h = 0.667) than any of the four localities in Texas, where Matagorda (h = 0.284) was more variable than the rest, which are devoid (h = 0.000) or nearly devoid (h ≤ 0.100) of variation (Table 9). Despite the overall low levels of genetic variability, pairwise F_ST in P. pugio revealed significant differences (p < 0.01) among the majority of pairwise comparisons, except between Matagorda and Pt. Arthur (Table 9).

SAMOVA results (

Table 10. SAMOVA results for P. pugio based on 466 bp of 16sRNA sequence from the northern Gulf of Mexico. The highest Among Groups (F_CT) value was obtained with five groups, as follows: Population 1: Sportsman’s; Population 2: Pt. Mansfield; Population 3: Matagorda and Pt. Arthur; Population 4: Ocean Springs; and Population 5: Venice.

Source of Variation	d.f.	Sum of Squares	Variance Components	Percentage of Variation
Among Groups	4	104.55	1.1235 Va	80.72
Among Populations Within Groups	1	0.18	−0.0049 Vb	−0.35
Within Populations	113	30.87	0.2732 Vc	19.63
Total	118	135.60
Fixation Indices			p-values (≥)
		FSC:	0.49365 ± 0.01428
		FST:	0.00000 ± 0.00000
		FCT:	0.06940 ± 0.00845

) returned the highest F_CT value (non-significant) for five distinct P. pugio populations, in agreement with the significant F_ST values. However, due to the low haplotypic diversity in the localities from E.Tex, it was necessary to pool those samples in order to obtain meaningful estimates of both historical demography and N_ef. The mismatch distribution of pairwise differences for P. pugio (Figure 9) indicates that about 22% of the shrimp share haplotypes (i.e., zero differences), with an additional 50% of the individuals differing by 1–3 mutations, and the rest by 4–5 pairwise differences. Both D and R2 tests were non-significant for (p > 0.05), as reflected by the multimodal shape of the curve. The estimate of N_ef was between 1.9–2.3 million individuals (

Table 11. Historical population demography parameters and estimates of female effective population size for P. pugio and P. marginata. Da, Tamura–Nei corrected genetic distance between the population and sister species, F. heteroclitus; τ, estimated mutational time since population expansion; μ, mutational rate per million years; D, Tajima’s D with probability value (p); R2, Ramos-Onsins and Rozas’s R2 with probability value (p); N_ef, estimated female effective population size.

Species	Da	τ	μ	D (p)	R2 (p)	Nef
P. pugio	0.0695	1.435	0.009	0.893 (0.824)	0.126 (0.832)	2.359
			0.011			1.930
P. marginata	0.0470	0.258	0.027	−0.927 (0.214)	0.043 (0.168)	0.444

). A Mantel test could not be calculated for P. pugio due to the lack of variability in the two E.Tex sampling locations mentioned above. Similarly, the absence of variation prevented conducting meaningful PCA plots for P. pugio.

3.3. The Planthopper P. marginata

P. marginata specimens (n = 124) from six sampling locations were successfully sequenced for 372 bp of COI. This species displayed less diversity (h < 0.45) and fewer haplotypes (n = 4;

Table 12. Molecular indices for P. pugio 16sRNA sequences by sample location along the northern Gulf of Mexico. M, no. of haplotypes; h, haplotypic diversity; π, nucleotide diversity; SD, standard deviation; S, no. of segregating (polymorphic) sites; Ts, no. of transitions; Tv, no of transversions; I/D, no. of insertions and/or deletions.

Location	N	M	h (SD)	π (SD)	S	Ts	Tv	I/D
South Padre	20	1	0.000 (0.000)	0.000 (0.000)	0	0	0	0
Corpus Christi	19	2	0.409 (0.100)	0.001 (0.001)	1	1	0	0
Matagorda	20	2	0.442 (0.088)	0.001 (0.001)	1	1	0	0
Reitan Marsh	24	1	0.000 (0.000)	0.000 (0.000)	0	0	0	0
Sportsman’s	21	3	0.267 (0.120)	0.001 (0.001)	2	2	0	0
Pt. Arthur	20	3	0.279 (0.012)	0.001 (0.001)	2	2	0	0
All Samples	124	4	0.248 (0.049)	0.001 (0.001)	3	3	0	0

) than both F. grandis and P. pugio, yet the levels of variation and the phylogeographic association were sufficient to reveal significant differences in haplotypic diversity among samples with more than 26% of the variance explained among groups (

Table 13. SAMOVA results for planthopper COI sequences for northern Gulf of Mexico samples. The highest Among Groups (F_CT) value was obtained with two groups, as follows: Population 1: Corpus Christi and Matagorda; Population 2: South Padre, Reitan Marsh, Sportsman’s, and Pt. Arthur.

Source of Variation	d.f.	Sum of Squares	Variance Components	Percentage of Variation
Among Groups	1	2.20	0.399 Va	26.39
Among Populations Within Groups	4	0.27	−0.0023 Vb	−1.50
Within Populations	118	13.40	0.1135 Vc	75.11
Total	123	135.60	1.31703
Fixation Indices			p-values (≥)
		FSC:	0.64027 ± 0.01649
		FST:	0.00000 ± 0.00000
		FCT:	0.07136 ± 0.00884

). Corpus Christi and Matagorda, which showed similar values of h, were more variable than the majority of planthopper localities sampled, whereas South Padre and Reitan Marsh were less variable (h = 0) than any other locality. In addition, these two geographically distant sampling locations shared the same haplotype (Figure 10), and thus were not different from each other (

Table 14. Values for pairwise comparisons of P. marginata samples from the northern Gulf of Mexico estimated from 372 bp of COI. Z-scores from Salicrú Χ² test for pairwise comparisons of haplotypic diversity are above the diagonal. Pairwise F_ST values are below the diagonal. Significant values are in bold, with significance at p < 0.05 denoted by *, and significance at p < 0.01 denoted by **.

	South Padre	Corpus Christi	Matagorda	Reitan Marsh	Sportsman’s	Pt. Arthur
South Padre		−4.090 **	−5.023 **	0.000	−2.225 *	−23.250 **
Corpus Christi	0.93285 **		−0.248	4.090 **	0.909	1.291
Matagorda	1.000000 **	0.96026 **		5.023 **	1.176	1.835 *
Reitan Marsh	0.98361 **	0.02105	0.99010 **		−2.225 *	−23.250 **
Sportsman’s	0.45823 **	0.81601 **	0.65339 **	0.86007 **		−0.099
Pt. Arthur	0.46252 **	0.40867 **	0.73527 **	0.47953 **	0.41360 **

). In each of the six locations sampled, a common haplotype accounted for >85% of the individuals, with three additional haplotypes, one mutational step away from the main haplotype, accounting for the remaining individuals in four localities (Figure 10). The F_ST between Corpus Christi and Matagorda was not significant after corrections for multiple testing, but these two samples differed, respectively, from South Padre, Reitan Marsh, Sportsman’s Road, and Pt. Arthur (Table 14). The Slatkin’s linearized F_ST and Reynold’s distance for P. marginata could not be calculated due to sampling locations with h = 0.

The highest F_CT value for P. marginata in SAMOVA was obtained by placing S. Padre and E.Tex in the same group (Table 13). However, this optimization must be questioned, since the index was not significant, and there is no biological rationale for grouping these two geographically discreet samples, which are nearly equidistant from the intermediate locality of Corpus Christi. Although no temporal samples were obtained to verify the stability in haplotype frequency, the observed patchiness of the genetic signature among locals may result from female variance in reproductive success within demes [96] that may be subject dramatic changes in population size and that may include local extinction and recolonization events. As such, the demographic estimates at the local level may not be very informative given that four out of six localities had less than two haplotypes, and the remaining two had only three haplotypes each. Accordingly, to obtain a regional reconstruction of the historical demography and N_ef of P. marginata, samples were pooled together, and a pairwise mismatch distribution for P. marginata was obtained (Figure 9). Both D and R2 tests were not significant (p > 0.05), suggesting that these populations were not subject to a population bottleneck followed by expansion, and the curve fit was concordant with a stable population at mutation drift equilibrium. The long-term N_ef was estimated at 444,000 individuals (Table 11). A Mantel test could not be calculated for P. marginata due to low levels of variation across sampling locations. Similarly, a PCA is not shown for P. marginata, as there was not enough variation to generate meaningful plots.

4. Discussion

4.1. The Gulf Killifish F. grandis

This study aimed to compare the levels of genetic variation among the populations of three ecologically important residents of the Spartina salt marshes found in estuaries along the northern central and west coast of the Gulf of Mexico. These species were selected because they display contrasting life histories, and therefore substantially different genetic signatures could be expected when analyzing the respective patterns of mtDNA sequence variation. In addition to adding to current knowledge of the patterns of connectivity of these species, this information was expected to shed light on the corresponding historical demographic signatures and the timing of events that may have influenced the populations of these species in the same region.

Our results show that F. grandis are highly variable in all the localities sampled, except in Venice, LA, near the mouth of the Mississippi River. F. grandis showed high levels of genetic structuring, with nearly 40% of the total variance explained by differences among the three regions of the Gulf surveyed (Table 5). F. grandis displays high site fidelity, with individuals typically only moving ~100 m between connected salt marsh sites through their lifetimes [36]. The sedentary nature of F. grandis, in conjunction with self-adhesive and demersal eggs, high predation rates, and a limited (expected) lifespan, limits overall movements, thereby reducing gene flow [97]. Williams et al. [55] concluded from the patterns of isolation (found using distance, spatial autocorrelation, and assignment tests derived from microsatellite data) that dispersal of F. grandis is limited, occurring primarily between neighboring sites. However, our mtDNA data indicate that dispersal over long periods of time is sufficient to overcome the genetic separation among adjacent marshes within regions, such that along the northern Gulf of Mexico, from Mississippi to South Texas, there are at least three populations that are at migration–drift equilibria, separated from each other by barriers to gene flow.

The subdivision of F. grandis is supported by the strong phylogeographic associations of mtDNA lineages and phylogroups, and by the distinct historical demographic signatures that exist within the region. Previous investigations using allozymes and restriction fragment length polymorphisms (RFLP) of total mtDNA indicate that Florida, Louisiana, and Texas F. grandis populations are more similar to each other than to Mobile Bay, Alabama, which stands as an outlier [97]. The distinctiveness of the Mobile Bay population relative to Florida was corroborated in a more recent study using microsatellites [55], although principal component analysis and Bayesian clustering revealed that Mobile Bay has a closer affinity to the samples of Louisiana and Texas. Herein, while we did not sample Mobile Bay, we did characterize Ocean Springs, MS, which lies about 40 miles away. This sample, along with Venice, LA, contained a subset of private lineages belonging to Phylogroup I. Unfortunately, it is not possible to determine the relationship of this phylogroup to those characterized by Gricius [97] for two reasons. First, in that study, total mtDNA was digested with six base cutter restriction enzymes, and the fragments were separated through agarose gels subjected to Southern blot hybridization using an F. heteroclitus probe. As a consequence, the number of fragments characterized in that study was low, resulting in a reduced number of haplotypes. For instance, in Corpus Christi, Gricius [97] only identified 2 haplotypes among 16 individuals, and 15 of these shared the same haplotype (h = 0.125); however, here, mtDNA sequence data rendered 15 haplotypes among 26 individuals (h = 0.945) for that locality. Secondly, the parsimony network presented by Gricius [94] was not rooted against the outgroup (i.e., the Mummichog, F. heteroclitus), and consequently, the position of lineages relative to our data cannot be determined. Characterization of additional mtDNA sequences from Mobile Bay and Florida and also from northern Mexico, as suggested by Williams et al. [55], is needed to unravel the sequence of events that gave rise to the phylogenic history of F. grandis along the Gulf coast, as it diverged from the Mummichog.

By using mtDNA sequence data, the current study was able to further subdivide the northern F. grandis population to the west of Mobile Bay into three units, with phylogeographic breaks placed roughly at the mouth of the Mississippi River and the upper Laguna Madre region of the Texas coast, north of Corpus Christi. These findings are concordant with the Mississippi River acting as a barrier to gene flow that results in a phylogeographic break for many coastal marine species in the Gulf, including fish and elasmobranchs [98,99,100], and studies that demonstrate a phylogeographic break at the hypersaline Laguna Madre system of Texas for oysters [101] and fiddler crabs [102].

The high levels of haplotypic diversity of F. grandis samples are indicative of large population sizes [103], congruent with previous studies on this species [55], and estimates of contemporary N_ef in this study ranged from 10 s to 100 s of millions. To gain an understanding of how F. grandis reached such high numbers, the historical demographic data were analyzed in two different ways. The first approach was to obtain estimates for the groups that yielded the highest amount of among-group variance (F_CT) with SAMOVA (Table 5), which consisted of three populations: (1) S.Tex (Corpus Christi); (2) E.Tex (Matagorda, Reitan Marsh and Sportsman’s Road); and (3) CN.Gulf (Venice and Ocean Springs). Mismatch distributions [91,92,104,105] were used to provide insight into historical demography for each group, irrespective of the phylogeographic association of the two distinct mtDNA phylogroups (Figure 6). The corresponding mismatch distributions for S.Tex and E.Tex were multimodal [91,105], and the results of D and R2 tests (Table 6) were non-significant (p > 0.05), which is indicative of large populations at equilibrium [93,94]. Despite a large number of major tropical storms and oil spills in the Gulf in recent years, these populations appear not to have suffered dramatic bottlenecks recently, corroborating that the effects of catastrophic hurricanes, such as Katrina and Ike, may not be as influential on coastal fish assemblages as previously thought [106]. Such resilience may be partially associated with the F. grandis’ benthic feeding behavior, shared with other estuarine fish, which keeps them closely associated with their preferred habitat, rather than moving to avoid the presence of oil [107]. By contrast, the shape of the mismatch distribution and the associated significance of CN.Gulf to D and R2 (p = 0.003, 0.001, respectively) suggests a population bottleneck followed by sudden expansion [91,105].

The presence of divergent mtDNA lineages whose origin can be assigned phylogeographically elsewhere increases the risk of misinterpreting a local multimodal signature as a a large stable population [24]. It is thus recommended to analyze the signatures separately by phylogroups or by clade. Accordingly, pairwise distributions of F. grandis were estimated separately by phylogroup, and then regionally by phylogroup. For Phylogroup I, only E.Tex was reanalyzed in this way, since the CN.Gulf only contains members of a subgroup of Phylogroup I (Figure 3 and Figure 4), and those results have been presented (Table 6). For Phylogroup II, estimates were calculated for E.Tex and S.Tex separately. Regional estimates of N_ef by phylogroup were four times smaller than with all the data combined; still, estimates of 10 s to 100 s of millions of females for each population were obtained (Table 7). Mismatch distributions for Phylogroup I for E.Tex and CN.Gulf, or for the analysis of all the lineages from all localities (Figure 7) combined, are congruent with historical bottlenecks followed by population expansion (D and R2 tests, both p < 0.05) that occurred approximately 170–180 thousand years bp. By contrast, mismatch distributions for Phylogroup II, whether conducted for all lineages and localities together or separately for E.Tex and S.Tex, display multiple peaks, and non-significant D and R2 tests suggest long demographic histories at equilibrium for populations in this region. Williams et al. [55] also concluded that their microsatellite data suggest that F. grandis populations in the western Gulf may be at or near migration–drift equilibrium at a regional scale, but that dispersal barriers and potential historical signatures on population structure will need to be taken into consideration at larger spatial scales.

Several hypotheses on how Pleistocene and Holocene climatic events modeled the phylogeographic patterns of estuarine coastal species found in the northern Gulf have been porposed, many of which are discussed in great detail by Barnwell and Thurman [102]. Delcourt and Delcourt [108] hypothesized that at the peak of the Wisconsin glaciation (~18,000 years ago), a eustatic sea level drop of about 120 m below its current level occurred, causing most salt marsh habitats to disappear from the northern Gulf except along isolated patches of the Texas continental shelf and southwestern parts of Florida. Additionally, Florida’s land mass expanded and created a cooler, more arid environment around the southern tip of the peninsula, allopatrically separating Gulf and Atlantic populations [26,109]. The resulting isolated patches in Florida and Texas became refugia that eventually served as sources to colonize new marsh habitats formed when the continental glaciers retreated, causing the sea level to rise during the Holocene. The longer demographic history of Phylogroup II, as indicated by higher τ values and the stability of the S.Tex and E.Tex regions, could be explained by the resilience of isolated patches of salt marsh habitat throughout the last glacial maxima (LGM), but perhaps extending to the glacial maxima and minima over the past 1.6 million years. The sudden expansion of the CN.Gulf population after being subject to a bottleneck could be explained by a founder event, as individuals from refugia elsewhere arrived in this region.

These hypotheses only partially explain the relationships of haplotypes revealed by the MJN and the historical demography patterns of the two phylogroups of F. grandis. The lineages belonging to Phylogroup I appear basally relative to F. heteroclitus, but a phylogeographic tree (Figure S1) indicates this basal position in relation to the sister species used as the outgroup corresponds to the three lineages identified as intermediate in the MJN. This indicates that the oldest lineages found in the Gulf, from Mississippi to Texas, are found towards the west (S.Tex and E.Tex). This is concordant with some models of Pleistocene coastal features and currents, which indicate net movement starting along the Mexican coast, moving north towards Texas, and then east along the Gulf coast [102]. Williams et al. [55] found a significant negative relationship between genetic diversity and latitude, a pattern consistent with the presence of hypothesized refugia in the southern Gulf regions during the Pleistocene that later recolonized the northern Gulf. Mitochondrial sequences, however, do not support this pattern, as no correlation (positive or negative) was found between haplotypic diversity and latitude (not shown). An alternative explanation centers on the oceanographic properties of the upper Laguna Madre. Studies of other coastal species in the Gulf show strong phylogeographic breaks at the Laguna Madre, leading to the hypothesis that this hypersaline system is a barrier to gene flow [26,101,102,109]. Additional sampling to include locations in Florida and Mexico is necessary to resolve the patterns of historical gene flow for F. grandis.

4.2. Grass Shrimp P. pugio

According to F_ST values and SAMOVA, P. pugio populations are genetically structured in the north and west regions of the Gulf of Mexico (Table 9 and Table 10). However, at a regional level, the groupings that explain the largest proportion of variance make little sense geographically or biologically, as they do not correspond to a migration–drift equilibrium model. Flowers [59] characterized the distribution of genetic variation in P. pugio along the U.S. Atlantic Coast, from S. Carolina to northeast Florida, contained in 16sRNA using SSCP, and identified six haplotypes. The dominant haplotype A was present at frequencies between 66–100% across nine of the ten localities sampled. The exception was St. Mary’s River, GA, where it only accounted for 8%, with most P. pugio containing haplotype B (62%) or haplotype F (31%). Haplotype B was present only in the southern portion of the 350 km of coastline sampled, and its frequency increased from the Ogeechee River in the north towards St. Mary’s River in the south. Conversely, the frequency of haplotype A increased towards the north, and included two localities in which that haplotype was fixed. As a result of this cline, a significant Mantel test, concordant with IBD, was reported along the east coast by Flowers [59]. Such a pattern was not found along the more than 1200 km of Gulf coastline surveyed here. Further, twice as many haplotypes were found in the Gulf than on the east coast, and while this could be explained by the higher resolution of direct sequencing compared to SSCP, or by the much longer stretches of sequence or coastline surveyed here, but not by the geographic manner by which this variation is distributed. Specifically, no single locality along the east coast contained as much of the overall mtDNA variation as that recorded in the Gulf locality of Ocean Springs, MS (see below). In the Gulf, there is no single dominant haplotype, as exemplified by haplotype A along the east coast. However, the Gulf and the east coast are similar in having two localities each devoid of mtDNA variation (h = 0). However, the two depauperate east coast localities of Moon River and Lazarretto Creek in northern Georgia are in close proximity (<30 km) to each other, whereas Pt. Mansfield and Galveston Island (Sportsman’s Road) are >500 km apart, with Matagorda Bay (h = 0.284), the most variable location in Texas, between them.

While until now, the focus has been on explaining genetic variation as a function of the connectivity among locations, patterns of gene flow at the local level should also be explored. Using allozymes, Fuller [57] described population structuring at a small spatial scale when comparing P. pugio living in ponds on Galveston Island to semi-open and open systems connected to Galveston Bay. Specifically, P. pugio living in closed systems (ponds) were less diverse and displayed stronger signals of population differentiation compared to those in channels connected to the Bay, with the largest diversity in localities open to the Bay, where gene flow was expected to be more substantial. Here, none of the localities sampled were isolated ponds, and thus the reported patterns of variation would correspond to semi-open or open systems. In light of these findings, it is important to note that the overall low levels of haplotypic diversity reported here and by Flowers [59] may be largely due to the low levels of variability contained in 16sRNA, which may be similar to Penaeid shrimp, wherein low levels of genetic variation at both COI and 16sRNA loci have been reported [110]. Accordingly, the characterization of variation at the nuclear level needs to be investigated.

Despite reduced overall levels of genetic variation (τ = 1.435; π = 0.001) observed in P. pugio, the data suggest the presence of multiple demes along the northern and western portions of the Gulf of Mexico, each with a very distinct historical demographic signature. Ocean Springs, for instance, is the most variable locality surveyed (h = 0.811), with Venice, the second most variable (h = 0.667). When these two localities are pooled together, they virtually encompass all the variation documented by the MSN (Figure 8) along the entire region of the Gulf surveyed, from Mississippi to Texas. Such depauperate levels of variation in Texas described herein collectively represent only a snapshot of the total mtDNA variation observed for the northern Gulf, and could be explained by either founder events [111,112] or sweepstakes in reproductive success (i.e., SRS or the Hedgecock effect) [96]. P. pugio females are highly fecund, capable of producing 100 s of eggs per spawning event, which are repeated multiple times per season [44,45,46]. Additionally, while actual population densities of P. pugio have been reported to be as high as 1.2 million individuals per 0.01 km² in a single marsh [44], this study estimated N_ef to be 1.9–2.4 million for the entire range sampled, from Mississippi to Texas (Table 11). This dramatic reduction in N_ef in comparison to actual population numbers, combined with high fecundity and Type III life history of P. pugio, is concordant with expected outcomes based on the SRS hypothesis. Accordingly, the results of this study indicate that P. pugio in Texas have not reached the migration–drift equilibrium displayed by F. grandis in the same region.

4.3. The Planthopper P. marginata

We found that among 124 COI sequences, only four very closely related (τ = 0.258) haplotypes were discovered, which included a common haplotype found at frequencies of 85% or higher. Despite such low levels of variation, P. marginata displayed higher levels of genetic population structure along the north and west Gulf of Mexico than F. grandis, although not as pronounced as P. pugio. Within E.Tex, none of the localities differed from each other, but individually, they differed from the S.Tex localities of Corpus Christi and Matagorda, but not from South Padre, which lies farther to the south, thus explaining why SAMOVA groups South Padre with E.Tex, a grouping that makes little sense geographically or biologically. Using the same primers for COI, Denno [60] characterized 53 P. marginata individuals from 15 populations, 7 from the Atlantic US Coast and 4 from the Gulf, in addition to putatively introduced populations from California and Portugal, using P. dolus and Toya venilia as outgroups. Their phylogenetic analysis unraveled an extensive geographic structure among native North American populations of P. marginata, characterized by strong phylogeographic associations. P. marginata haplotypes clustered into two well-supported sister clades, one comprising the mid-Atlantic coast (Virginia to New York) lineages, and another of south Atlantic (South Carolina to northern Florida) and Gulf Coast lineages (Figure 11).

Within the South Atlantic–Gulf Coast clade, there is only one well-supported subgroup (bootstrap > 87) of haplotypes specific to the western Gulf (Louisiana and Texas). The rest of this clade consists of three closely related haplotypes, one found in South Carolina, Florida, and Mississippi, another one from Florida, and a third haplotype found in Mississippi, but also in Virginia Beach, VA, where the mid-Atlantic clade dominates. Accordingly, the variation characterized here matches the well-supported groups of western Gulf P. marginata haplotypes. This is relevant because rather than concluding that COI is not sufficiently variable to characterize P. marginata populations, it instead illustrates that the northern Gulf from Ocean Springs, MS to Pt. Mansfield, TX, contains only a small fraction of the total variation contained in that mtDNA gene across the range of this species (Figure 11), whose origins can be traced to the vicariant event that took place during the Pleistocene and which affected a host of coastal marine species along the US East Coast and Gulf of Mexico [26,102,108], thereby helping to model the patterns of genetic variation observed in P. marginata but also in F. grandis, as documented above. However, whereas allopatry during the Pleistocene resulted in anagenesis for Fundulus, resulting in F. grandis-F. heteroclitus sister species, in P. marginata, it gave rise to the mid-Atlantic and South Atlantic–Gulf Coast clades. Further, the reduced levels of genetic variation observed in the northern Gulf are likely the result of past losses associated with founder events, as the western portion of the Gulf of Mexico was colonized from the east. With time, these populations became isolated from the eastern end of this basin, and new variants evolved, giving rise to a strongly (>87%) supported western group, which in turn forms part of the strongly supported (94%) South Atlantic–Gulf Coast clade that separated from the mid-Atlantic clade during the Pleistocene.

While the above phylogeographic reconstruction of events may explain the establishment of P. marginata mtDNA lineages found in the western portion of the northern Gulf of Mexico, it does not account for the patchiness in the geographic distribution of such variation within this region. The population structure of planthoppers in general, and of P. marginata in particular, may be strongly influenced by the plasticity of their life histories, as P. marginata can switch between a fecund, wingless morph (brachypterous) and a less fecund winged morph (macropterous) with high dispersal potential [22,53]. While previous studies suggest that the brachypterous form is most common in salt marshes along the Gulf Coast [54,113], our collections consisted almost exclusively of the macropterous form. The higher dispersal potential of macropterous forms may increase gene flow between neighboring marshes, and weaken genetic signals of population structure [114,115]. While the presence of a dominant allele could result from high levels of gene flow, that same explanation does not necessarily account for the low levels of genetic diversity observed, which instead may be due to the slow rate of evolution of the COI gene in Hemiptera. A study of COI sequences in 344 species belonging to that order found low levels of intra-specific sequence divergence (<2%) in the majority of the species sampled [116]. Yet, while a low mutation rate may explain the reduced number of haplotypes found in P. marginata along the Texas coast, the local absence of genetic variation in two out of six localities surveyed suggests that like P. pugio, local extinctions followed by recolonization may be responsible for observed genetic patchiness in the overall distribution of variation.

5. Conclusions

Comparison of the patterns of genetic variation in three species inhabiting Spartina salt marshes in the northern Gulf of Mexico provides insight into the forces that modeled such patterns within species and differentially among species. While recognizing that direct comparison of the three species is biased using different mtDNA loci and sequence lengths employed to characterize each species, in all instances, the resolution was sufficient to examine the general patterns of population differentiation and to reconstruct phylogeographic associations and historical demography. F. grandis mtDNA sequence data generated for this study found evidence of IBD and of population structure in the Gulf west of Mobile Bay, Alabama, that had not been found in previous studies of this species. The E.Tex population, in particular, can be further subdivided using the phylogroup associations found in this study, resulting in two subgroups within E.Tex that have distinct patterns of historical demography. These different signatures may be attributed to a succession of historical breaks between the Gulf and Atlantic during periods of glaciation in the Pleistocene. An alternative hypothesis is a historic phylogeographic break between S.Tex and E.Tex, during which time the E.Tex population underwent an expansion and colonized the CN.Gulf. Following this colonization, contact between S.Tex and E.Tex was reestablished, generating the two subgroups now present in the E.Tex population. Different theories on how the Pleistocene climate may have affected phylogeography in the northern Gulf could support either hypothesis. Therefore, future sampling efforts to include populations from Mexico and Florida are necessary to resolve the patterns of historical demography described in this study.

The mtDNA sequencing data generated for both P. pugio and P. marginata revealed relatively low levels of haplotypic diversity within the region sampled. For both, analyses of the hierarchical distribution of variation indicated population structuring, although the suggested groupings that maximize such variance do not make sense geographically or biologically, with no evidence of IBD. For P. pugio, the levels of variation in the eastern end of the northern Gulf region sampled encompass almost all the variation characterized in this study for this species, including that found in the western regions (E.Tex and S.Tex), wherein individual localities each represent only a snapshot of all the potential genetic variation present throughout the entire range surveyed. This suggests that each locality along the Texas coast may be subject to extinction and colonization events, or severe fluctuations in population size, which when combined with variance in reproductive success result in individual pieces of a mosaic representative of the overall variation. Considering the life history details of P. pugio, the observed patterns are consistent with the SRS hypothesis. For P. marginata, haplotypic diversity is extremely low, with most individuals sharing a single haplotype. However, such depauperate levels of variation found in the western portion of the Gulf represent only a small fraction of the overall mtDNA variation contained in P. marginata, since towards the east Gulf and along the Atlantic coast, substantial levels of genetic variation exist. Accordingly, the genetic signature of P. marginata suggests colonization of the Texas coast from the east, but sufficient time has passed to establish a distinct set of lineages in the western portion of the northern Gulf. Further investigations targeting nuclear markers with a greater degree of polymorphism would be beneficial for determining a finer-scale population structure along this geographic range.