Estuaries form a series of unique wetland habitats isolated from each other, often facilitating genetic divergence among populations. The estuarine seablite Suaeda esteroa Ferren & S.A.Whitmore (Suaedoideae; Chenopodiaceae) is common in northwestern Mexican estuaries, where the spatial isolation from one another may promote diversification within this plant species. In this study, we created a novel nuclear ribosomal DNA genome skim dataset from 30 S. esteroa herbarium specimens collected from estuaries along the northwestern Mexican coast to assess genetic patterns within and among localities. We constructed maximum likelihood and Bayesian phylogenetic trees, and conducted individual- and estuary-level landscape genetics analyses. While our landscape genetics analyses provide evidence that more geographically distant individuals are more genetically distant, our phylogenetic tree and estuary-level analyses demonstrate surprising groupings of geographically distant estuaries. We hope our findings will encourage further investigation of gene flow among northwestern Mexican estuaries and promote conservation action within the region.
The spatially discrete nature of estuaries is known to promote population differentiation of the organisms that reside within them (Bilton et al. 2002, Palumbi 2003). Estuarine taxa often have low colonization ability due to the distance and inhospitable habitat between estuaries, reducing gene flow among populations and increasing genetic divergence due to genetic drift (Baggio et al. 2017). Phylogenetic analyses have been used to assess genetic relatedness among estuarine populations (Zhang et al. 2007, Sabry et al. 2013). However, depending on the dispersibility and colonizing ability of the organism, as well as how recently the estuaries formed, genetic differences among estuarine populations can be difficult to detect (Gold and Richardson 1998). Additional factors, including population size and mating systems, influence the rate of genetic differentiation generated within and among populations as a result of mutation and drift (Spieth 1974, Bawa 1992). Genetic differentiation among estuarine individuals and populations can be measured on a geographic scale to determine whether or not geographic proximity corresponds to genetic similarity (Chenoweth et al. 1998, Palumbi 2003). Isolation by distance analyses can be used to detect the accumulation of genetic divergences due to genetic drift, while calculating genetic differentiation metrics and conducting an analysis of molecular variance can indicate the degree of gene flow among populations and regions (Wright 1943, Meirmans 2006, Verity and Nichols 2014). Even when genetic differences among populations appear insignificant, combining phylogenetic analyses and tests for genetic differentiation can reveal patterns on a broad geographical scale (Bradburd et al. 2013).
Around 3.5 million years ago, the northwestern coast of Mexico saw the formation of over 100 estuaries. The formation of these estuaries has been linked to multiple occurrences of Pleistocene glaciation that caused significant changes in the sea level (Jacobs et al. 2004). Evidence of genetic differentiation among these estuaries has been observed in copepods and fishes, most likely because of spatial isolation and low migration rates among estuaries (Ganz and Burton 1995, Gold and Richardson 1998). While most estuarine studies have focused on aquatic taxa, the aridity of the terrestrial habitat between these estuaries may serve as a second barrier to migration and gene flow among populations of terrestrial plant taxa that are confined to estuary margins (Ferren and Schuyler 1980). The terrestrially and aquatically isolated estuaries along the coast of Baja California and Sonora, Mexico provide an opportunity to examine the distribution of genetic diversity among morphologically similar conspecific or con-generic plant taxa.
The genus Suaeda Forssk. ex J.F.Gmelin (Suaedoideae; Chenopodiaceae) has a worldwide distribution and is composed of approximately 100 species, with the highest species diversity occurring in temperate zones (Fisher et al. 1997, Schütze et al. 2003, Dehghani and Akhani 2009). Members of the Suaeda subgenus Brezia sect. Brezia occur in the New World, with several species inhabiting estuaries along the coasts of the Baja California peninsula and Sonora, Mexico (Brandt et al. 2015). Suaeda is one of the few halophytic genera present in these estuaries that is an obligate seed-producer (Hopkins and Blackwell Jr. 1977, Ferren and Roberts 2011). Most estuarine halophytic plant species reproduce vegetatively rather than sexually, due to the low rates of seed germination and seedling survival in hypersaline environments (Jefferies et al. 1979, Yuan et al. 2019). Sexual reproduction increases the chances of mutations and genetic drift, and therefore, the generation of genetic differentiation among populations (Eckert 2001). The sexual reproduction of Suaeda sect. Brezia species, combined with the terrestrial and aquatic isolation among estuaries, creates the potential for genetic differentiation among populations.
Over a 23-year (1978–2001) period, Wayne R. Ferren Jr. and colleagues collected more than 350 herbarium specimens of Suaeda from estuaries along the Pacific and Gulf coasts of Baja California and from the Gulf Coast of Sonora, Mexico. From this work, two new species, Suaeda esteroa Ferren & S.A.Whitmore and Suaeda puertopenascoa C.Watson & Ferren, were described (Watson and Ferren 1991). In his study of Suaeda sect. Brezia in northwestern Mexican estuaries throughout the 1990s, Ferren observed subtle differences in life history and morphology (i.e., branching patterns, leaf shape, leaf scar shape, seed type) among estuarine S. esteroa populations (Ferren and Roberts 2011). Ferren's natural history observations, coupled with the geographic isolation among the estuaries, provide strong justification for further study of the phylogeographic relationships of extant S. esteroa populations within northwestern Mexican estuaries.
In the present study, we assess intraspecific relationships among 30 S. esteroa specimens sampled in nine estuaries of Baja California and Sonora, Mexico. We hypothesized that S. esteroa specimens from the same or neighboring estuaries would have greater genetic similarity than those from estuaries further away. Additionally, we predicted that the Baja California peninsula acts as both a terrestrial and aquatic barrier to gene flow, resulting in a genetic break between populations on the Pacific Ocean side and within the Gulf of California. To explore these hypotheses, we use nuclear ribosomal DNA (nrDNA) sequence data to assess the phylogenetic structure of S. esteroa populations and to evaluate genetic differentiation among individuals and estuaries as a function of geographic location. We discuss the patterns of genetic relatedness among these estuarine populations and its implications for conservation actions.
Materials and Methods
Taxon Sampling
All samples (n = 30) were taken from specimens collected between 1985 and 2001 and curated by the University of California, Santa Barbara Herbarium (UCSB) at the Cheadle Center for Biodiversity and Ecological Restoration ( Appendix 1 (Motta_seablite_Appendix1.jbw.docx)). Each herbarium specimen contained one individual. We sampled one to four herbarium specimens apiece of S. esteroa collected from the nine estuaries, determined by Ferren to be morphologically distinct based on differences in branching patterns, leaf shape, leaf scar shape, and seed type. We consider that the nine estuaries contain nine distinct populations of S. esteroa. All 30 individuals are included in the phylogenetic tree and Mantel tests described below; while each specimen is associated with a specific estuary, geographic location was not considered in phylogenetic analysis of these individuals. By contrast, in the landscape genetics analyses conducted to calculate Nei's GST as a metric of genetic differentiation, individuals were grouped by estuary and by region to evaluate genetic differentiation among estuaries. While internal transcribed spacer (ITS) sequences of ten S. esteroa individuals from these estuaries were already available (Brandt et al. 2015), our use of genome skimming resulted in a much longer nrDNA sequence; therefore, previously sequenced samples were not included in our analysis. Between one and four herbarium specimens were selected from each of the following estuaries: Las Animas, Las Lisas, Los Angeles, San Carlos, San Felipe, San Gregorio, San Ignacio, Santa Cruz, and Santa Rosa (Fig. 1).
Fig. 1.
Map of localities where Suaeda esteroa was sampled in our study. Each estuary is associated with a unique shape while the colors denote the three main groups we observed: San Ignacio and San Gregorio (green); Santa Cruz, Santa Rosa, Las Lisas, San Carlos, and San Felipe (yellow); Los Angeles and Las Animas (purple). Shapes and colors are intended to facilitate identification of geographic placement of populations on the corresponding phylogenetic tree (Fig. 2) and heatmap of Nei's GST values (Fig. 3).

Fig. 2.
Unrooted phylogenetic tree of Suaeda esteroa derived from the best tree of a maximum likelihood (ML) analysis of nrDNA. Bootstrap values are indicated before the slash (/); posterior probabilities from a consensus Bayesian analysis are shown following the slash (/) and any discrepancy between the two analyses is indicated by an asterisk (*). Each monophyletic group has been labeled as Clade 1–7 and individuals are associated with symbols, the shape of which denotes the estuary from which it was collected and the color denotes the three main groups we observed: San Ignacio and San Gregorio (green); Santa Cruz, Santa Rosa, Las Lisas, San Carlos, and San Felipe (yellow); Los Angeles and Las Animas (purple). Shapes and colors are intended to facilitate identification of geographic placement of populations and correspond to those used in Fig. 1 and Fig. 3.

Fig. 3.
Heatmap depicting pairwise values of Nei's GST of Suaeda esteroa among estuaries calculated using nrDNA sequences. The shape denotes the estuary from which it was collected and the color denotes the three main groups we observed: San Ignacio and San Gregorio (green); Santa Cruz, Santa Rosa, Las Lisas, San Carlos, and San Felipe (yellow); Los Angeles and Las Animas (purple). Shapes and colors are intended to facilitate identification of geographic placement of populations and correspond to Fig. 1 and Fig. 2.

DNA Extraction and Sequencing
Genomic DNA was isolated from herbarium tissue using a GeneJETTM Plant Genomic DNA Purification Mini Kit (Thermo Scientific). Library preparation was performed by Global Biologics LLC (Columbia, Missouri, USA). All samples were barcoded, pooled, and pair-end sequenced on an Illumina HiSeq 2500 for 100 rapid cycles. We used genome skimming as opposed to other high-throughput methods because DNA isolations from 30-year-old herbarium specimens yielded sufficient DNA quality and quantity for this approach, the bioinformatics were straight-forward, and subsequent analysis was possible at low cost (Ripma et al. 2014, Simpson et al. 2017, Nevill et al. 2020).
DNA Sequence Reads Quality Control
Quality control of raw reads was conducted using PRINSEQ (Schmieder et al. 2011) to filter out reads using the following parameters: exact sequence duplicates, reads with a mean quality Phred score below 30, and reads with more than one N were removed (Ripma et al. 2014). These post-quality control reads were imported into Geneious v.11.1.5 ( http://www.geneious.com/) in FASTQ format.
Nuclear Ribosomal Cistron Assembly
A reference sequence for the nuclear ribosomal DNA was constructed using a 629-bp sequence from S. esteroa (FJ449791) with complete internal transcribed spacer region 1 (ITS 1), 5.8S gene, and internal transcribed spacer region 2 (ITS 2) obtained from GenBank. A reference-guided assembly of the read pool from one specimen of S. esteroa (Specimen ID 994; Appendix 1 (Motta_seablite_Appendix1.jbw.docx)) was executed in Geneious v.11.1.5 ( http://www.geneious.com/) with medium-low sensitivity, default settings, and 100 iterations. The resulting consensus sequence (6,706 bp) was saved using the highest quality threshold where sequence quality was used to call the best base, and areas with less than 10× sequence coverage were masked with Ns. Because the consensus sequence of Specimen ID 994 was ten times the size of the original reference sequence, this sequence was used as a reference to assemble the remaining read pools into sequences using medium/low sensitivity, 25 iterations, and default settings ( Appendix 1 (Motta_seablite_Appendix1.jbw.docx)). All consensus contigs were saved using the highest quality threshold where sequence quality was used to call the best base, areas with less than 10× sequence coverage were masked with Ns, and IUPAC ambiguity codes were retained. Sequences were annotated using the “Annotate and Predict” feature on Geneious from a database created by searching GenBank for “Chenopodiaceae internal transcribed region”. The following annotations with 50% or greater similarity to relatives were copied onto sequences: Suaeda nigra Raf. (MF963955) for 5.8S and the two ITS regions, Dysphania ambrosioides (L.) Mosyakin & Clem-ants (KY968902) for the 26S gene, and Spinacia oleracea L. (SPIRG18S) for the 18S gene. Sequences were aligned using the Geneious MAFFT plugin (version 7.450; Katoh et al. 2002, Katoh and Standley 2013) with default settings. Polymorphism Information Content (PIC) of the final alignment was calculated using the Geneious GARLI plugin (version 2.0; Zwickl 2006), which displays variable characters and PICs in the “info tab”.
Phylogenetic Analyses
Unrooted phylogenetic analyses were performed using maximum likelihood (ML) and Bayesian inference. Maximum likelihood phylogenetic inference was implemented using the RAxML Geneious plugin (v8.2.11; Stamatakis 2014) with a GTR+GAMMA model for nucleotide evolution. Statistical support was assessed using 1,000 rapid bootstrap (BS) replicates. Bayesian inference analyses were performed using the MrBayes Geneious plugin (3.2.6) using a GTR substitution model with gamma rate variation (Huelsenbeck and Ronquist 2001). A Markov chain Monte Carlo was initiated with 1,100,000 generations, sampling trees every 200th generation. Each run consisted of four heated chains using a heating parameter of 0.2. A Bayesian consensus tree with posterior probability values of nodal support was constructed by Geneious. Resulting trees were viewed in Geneious in an unrooted tree layout and formatted in Adobe Illustrator CS (Adobe Systems, San Jose, California, USA). Posterior probabilities from a consensus Bayesian analysis were mapped onto the best ML tree. Bootstrap values (BS) and posterior probabilities (PP) are indicated next to the relevant branch and are separated by a slash (e.g., BS/PP) and any discrepancy in topology between the two analyses is indicated by an asterisk (*).
Landscape Genetics Analyses
Aligned sequences of nrDNA were imported into R (4.1.1) as a DNAbin object using the ‘adegenet' package (Jombart and Ahmed 2011, R Core Team 2021). All 30 individuals were included to conduct individual-level analyses, and while each was associated with an estuary, the individuals were not grouped together. Pairwise genetic distances (%) among individuals were calculated using the ‘raw’ model from the ‘ape’ package (Paradis and Schliep 2019). The ‘raw’ model calculates the proportion or the number of sites on the genome that differ between each pair of sequences while excluding columns with ambiguity codes; this model yielded results similar to those of other models included in the ‘ape’ package. Euclidean distances (km) among the localities from which specimens were collected were calculated in R and oceanic distances (km) as well as shoreline distances (km) were measured by hand in Google Maps (Google Maps 2021). We considered the shortest distance across water to be the “oceanic distance” while the “shoreline distance” was defined as the minimum shoreline distance required to connect two localities. Mantel tests using the Pearson method were conducted using the ‘vegan’ package for genetic distance as a function of each of the following: Euclidean distance, oceanic distance, and shoreline distance (Oksanen et al. 2016). Results were visualized as pairwise scatter plots using the ‘ggplot2’ package (Appendix 2; Wickham 2011).
To analyze genetic differentiation on an estuary level, ‘adegenet’ was used to transform our DNAbin object of nrDNA sequence data into a genind object. Individuals were grouped by estuary to evaluate genetic differentiation among different estuaries. Each estuary was assumed to represent a distinct population and assigned to one of two regions (Pacific Ocean or Gulf of California) to allow us to test for genetic differences among the estuaries as well as between the two regions. Pairwise values of genetic differentiation metric Nei’s GST of the estuaries were calculated using ‘mmod’ and visualized as a heatmap using ‘pheatmap’ (Nei and Chakravarti 1977, Winter 2012, Kolde 2018). An analysis of molecular variance (AMOVA) was conducted using the ‘poppr’ package to test for variation within and among estuaries, as well as between regions (Pacific Ocean or Gulf of California) (Kamvar et al. 2014).
Results
Assembly, Depth of Coverage, and Library Content
Total nuclear ribosomal DNA (nrDNA) sequencing depths were between 42.9× and 26,064.5×, with an average depth of 1,196.8× (± SE 893.69). Between 0.068% and 12.8% of the overall readpool of each sample consisted of nrDNA, with a mean of 1.28% (± SE 0.43%). The nrDNA dataset alignment comprised a partial external transcribed spacer (ETS), a complete 18S gene, ITS 1, 5.8S, ITS2, and 26S gene (size range 6,124–6,706 bp) containing 0.48% polymorphism information content (PIC).
Phylogenetic Analyses
ML and Bayesian phylogenetic analyses resulted in nearly identical topology (Fig. 2). The one conflict, present in Clade 7, is due to the best ML tree having one additional, weakly supported node that was collapsed in the Bayesian analysis. The Bayesian analysis included stronger support values overall, especially at deeper nodes. Strongly supported, biogeographically intuitive groupings containing individuals from a single or neighboring estuaries include Clades 4, 5, and 7. Clade 4 contains all S. esteroa specimens collected from Los Angeles and Las Animas, two neighboring estuaries within the Gulf of California, and is well-supported (BS 93, PP 1). Three of the four specimens collected from Las Lisas appear within Clade 5, for which there is mixed support (BS 67, PP 0.99). The fourth individual collected from this estuary (Specimen ID 992, Las Lisas), however, is genetically isolated from the other three. Clade 7 contains all four specimens sampled from San Felipe with fair support (BS 71, PP 0.99). Specimens from estuaries along the Pacific Ocean and within the Gulf of California group together. Individuals from San Carlos, located on the Pacific Ocean coast, consistently group with individuals from Gulf of California estuaries Santa Cruz and Santa Rosa, in well-supported Clades 2, 3, and 6 (BS 99, PP 1; BS 99, PP 1; BS 63, PP 0.99, respectively). Specimens from two of the three populations on the Pacific Ocean coast, San Ignacio and San Gregorio, are found together in well-supported Clade 1 (BS 99, PP 1). Clade 1 also includes an individual from a Gulf of California estuary, Santa Cruz (Specimen ID 1061).
Landscape Genetics Analyses
All Mantel tests resulted in significant (P < 0.05) correlations between genetic distance and geographic distance. Genetic distance is best explained by shoreline distance (r2 = 0.35, P = 0.0014), followed by ocean (r2 = 0.34, P = 0.002), and Euclidean distance (r2 = 0.18, P = 0.039) (Appendix 2). The dendrogram produced with the heatmap to visualize pairwise values of Nei's GST contains three major clades that support the trends observed in the phylogenetic tree (Fig. 3). Estuaries Las Animas and Los Angeles are weakly differentiated from each other and are the most genetically distinct from the other populations, including other populations within the Gulf of California. Two of the three populations along the Pacific Ocean coast, San Ignacio and San Gregorio, have little to no genetic differentiation between one another. As observed in the phylogenetic tree, the third Pacific Ocean coast population, San Carlos, is genetically more similar to Gulf of California populations and has a lower GST value compared to Santa Cruz and Santa Rosa populations than to San Ignacio and San Gregorio S. esteroa. The AMOVA test did not reveal a significant difference in pairwise comparisons of Nei's GST values when populations were grouped by region and defined as located on the Pacific Ocean coast or within the Gulf of California (P > 0.05).
Discussion
Analysis of nrDNA sequences in this study indicated significant phylogeographic structure among the 30 samples of S. esteroa collected by Ferren in northwestern Mexican estuaries. Our phylogenetic analyses revealed several strongly supported clades of specimens from the same or neighboring estuaries (i.e., Clade 7, Clade 4; Fig. 2). However, our results did not indicate a clear distinction between individuals from populations along the Pacific Coast and within the Gulf of California. Most surprisingly, we found that individuals from the Pacific Coast estuary San Carlos consistently grouped with individuals from the Santa Cruz and Santa Rosa estuaries located within the Gulf of California (Figs. 1–3). The well-supported S. esteroa clade of the remaining two Pacific Coast estuaries, San Ignacio and San Gregorio, also includes an individual from Gulf of California estuary, Santa Cruz. Overall, our results demonstrate that gene flow has occurred between estuaries on either side of the peninsula.
Our individual-level landscape genetics analyses using nrDNA revealed overarching genetic trends due to isolation by distance, while the estuary-level landscape genetic analyses provide further support for our phylogenetic tree findings. All three methods of measuring geographic distance (Euclidean, oceanic, and coastal) resulted in a positive, significant correlation between genetic difference and geographic distance (Appendix 2). However, aquatic geographic distances better explained genetic distance than terrestrial geographic distance. Oceanic and coastal distance resulted in a higher positive correlation (r2) with genetic distance than Euclidean geographic distance between sampled specimens. The results align with existing literature on Suaeda, which suggests that its dispersal occurs through water, rather than land or wind (Brandt et al. 2015). Consequently, genetic distance can be more accurately predicted by considering the distance between bodies of water. While our Mantel tests indicate that more geographically distant individuals are more genetically distant, the pairwise values of Nei's GST group San Carlos with Santa Cruz and Santa Rosa, rather than with the geographically closer, neighboring San Ignacio and San Gregorio (Figs. 1 and 3). The AMOVA test does not indicate significant genetic differentiation (P > 0.05) between populations on the Pacific Ocean coast and those within the Gulf of California. Although our individual-level landscape genetics data support the idea that the arid, terrestrial environment acts as a barrier to dispersal, our findings do not confirm that the peninsula of Baja California determines genetic distinctness among populations.
The patterns reported here suggest that gene flow has occurred among the populations sampled in this study. However, our results do not indicate whether the genetic similarities observed among individuals sampled from geographically distant populations are the result of historical vs. contemporary processes or events (Feng et al. 2014). The trends observed in our phylogenetic tree and our pairwise analysis of Nei's GST may be due either to the retention of ancestral alleles, or to incomplete lineage sorting (ILS; Zhou et al. 2017, Cai et al. 2021). While these Mexican estuaries formed 3.5 Mya., the separation of congeneric Suaeda linearis and S. esteroa occurred approximately 0.15 Mya. (Brandt et al. 2015). The recent age of S. esteroa as a species makes ILS a probable explanation for the groupings of geographically distant individuals and populations that we observe (Fig. 2). Alternatively, the observed genetic similarities may be due to recent long-distance gene flow, given that there are discernable patterns rather than the similarities being constant across all estuaries, as would be expected with ILS (Takayama et al. 2008, Tomizawa et al. 2017). It is also possible that the phylogenetic patterns we observe in S. esteroa are caused by both ILS and gene flow (Blanco-Pastor et al. 2012, Kleinkopf et al. 2019). Comparative analysis of other sources of genetic variation, such as chloroplast and mitochondrial DNA, could address whether the observed clades are the result of historical (ILS) or contemporary (gene flow) events, or both (McCauley 1995, Tomaru et al. 1998, Petit et al. 2005, Degnan and Rosenberg 2009). Genome-skimming data also contain potentially valuable information beyond nuclear ribosomal, chloroplast, and mitochondrial DNA in the form of low-copy genes (Wolf et al. 2015, Berger et al. 2017). The inclusion of additional gene regions combined with statistical analyses could reveal the origin of occasional, yet unexpected genetic similarities observed between samples from opposite sides of the Peninsula (Maddison and Knowles 2006).
If the observed trends are the result of contemporary, long-distance gene flow, then an investigation into the dispersal and colonization capability of S. esteroa seeds may inform the interpretation of genetic similarities between geographically distant samples. In particular, empirical observations of the dispersal and establishment capabilities of S. esteroa seeds are needed. Similar to some of its congeners, S. esteroa produces dimorphic seed types: a larger, coiled, dull-brown seed type and a smaller, bi-convex, black-shiny seed type (Watson and Ferren 1991). Studies of other Suaeda taxa with dimorphic seeds have found that the two morphs differ in dispersal characteristics as a result of differences in germination rate, dormancy, and saline tolerance (Wang et al. 2008). The larger, coiled, dull-brown seed type germinates rapidly, is not dormant, and has a lower saline tolerance, making it adept for local, terrestrial dispersal. Meanwhile, the smaller, bi-convex, black-shiny seed type germinates relatively slowly, is buoyant, dormant, and has a higher saline tolerance, suggesting that it is adapted for long distance, aquatic dispersal. The congener, S. maritima, produces black-shiny type seeds that remain buoyant and viable for up to 90 days; however, we do not know if S. esteroa seeds are similarly resilient (Chang et al. 2008). Whether just one or both seed types are produced may vary among individuals and populations; variation in seed types might therefore affect the dispersal potential of particular populations (Ferren and Roberts 2011). For example, S. esteroa individuals from the San Felipe region have been observed to produce only the larger, coiled seed type, which is also viviparous, potentially limiting their dispersal potential, and contributing to the genetic distinctness of this population (Clade 7, Figs. 2 and 3; Ferren and Roberts 2011). A systematic study of S. esteroa seed dimorphism, buoyancy, and germination frequency would inform whether contemporary, long distance dispersal events likely contribute to gene flow among populations.
If S. esteroa seeds are capable of long-distance seed dispersal followed by the successful colonization of sites inhabited by conspecific populations, then ocean currents or bird migration paths may explain the trends we observe (Marinone 2003, Weising and Freitag 2007). Depending on the time of year, ocean currents travel from north to south along the Pacific Coast of the Baja California Peninsula and round the tip of the peninsula to enter the Gulf of California (Collins et al. 1997, Valle-Rodríguez and Trasviña-Castro 2017). Ocean currents therefore offer a potential explanation for the consistent grouping of San Carlos with estuaries within the Gulf of California (Fig. 2, Clades 2, 3, 6; Fig. 3). Meanwhile, geographic features other than the Baja California Peninsula may be more relevant in preventing gene flow. Punta Hughes separates San Ignacio and San Gregorio from San Carlos and potentially obstructs ocean-mediated dispersal, offering an explanation as to why S. esteroa specimens from these Pacific Coast estuaries do not group together (Figs. 1 and 3). If material from San Carlos is arriving in the Gulf of California with enough frequency that it consistently groups with individuals from Santa Cruz and Santa Rosa populations, it is also possible that material from San Ignacio and San Gregorio occasionally arrives within the Gulf. This possibility offers a potential explanation as to why a sample from Santa Cruz groups with samples from San Ignacio and San Gregorio (Fig. 2; Clade 1). Islands may also restrict aquatic seed dispersal; isolation of the neighboring estuaries Las Animas and Los Angeles from other estuaries by the Isla Angel de la Guarda potentially causes these estuaries to be more genetically distinct (Figs. 1 and 3). Transport of seeds by migrating waterfowl has also been suggested as an explanation in other genetic studies of Suaeda in which geographically distant individuals have been found to group together (Weising and Freitag 2007, Brandt et al. 2015). The peninsula of Baja California is an important part of the Pacific shorebird flyway in the autumn and many bird species overwinter in the estuaries along both the Pacific and Gulf side of the peninsula (Carmona et al. 2004). Migratory waterbirds are known to be important dispersal vectors for plants and have the potential to transport seeds well over 1000 km (García-Álvarez et al. 2015). A closer study of ocean currents and waterfowl migration is also necessary to determine dispersal potential among geographically distant estuaries.
Understanding whether and how gene flow is occurring among estuarine populations in northwestern Mexico is crucial, considering the growing number of anthropogenic threats faced by this region (Camacho-Ibar and Rivera-Monroy 2014). These estuaries have already undergone drastic changes as the result of human activities since these specimens were collected nearly 30 years ago. While directly affected by the development of mariculture facilities, destination resorts, ports, salt production facilities, and airfields, these coastal ecosystems are also being indirectly affected by sea-level rise and reduced rainfall as the result of climate change (Ibarra-Obando and Escofet 1987, Ortiz-Lozano et al. 2005). The multitude of threats to the estuaries of Baja California and Sonora creates an urgency to investigate the cryptic sources of genetic variation we detected here. We encourage a more expansive study of additional gene regions, S. esteroa seed dispersal capabilities, and potential dispersal mechanisms among estuaries to help us understand how these estuarine populations will be affected by ever-growing anthropogenic threats.
Acknowledgments
We thank the University of California, Santa Barbara Herbarium (UCSB) and the UCSB Cheadle Center for Biodiversity and Ecological Restoration for providing the herbarium specimens used in this study. We would like to express our gratitude to Dr. Marina Corrêa Côrtes for her input and recommendations about the landscape genetic analyses. We are also grateful for the financial support provided by the Undergraduate Research and Creative Activities Office at UCSB. SJM is grateful to the Yale Institute of Biospheric Studies (at Yale University), which provided critical support during her 2019–2020 sabbatical, when this manuscript was completed.
Literature Cited
Appendices
Appendix 1
Appendix 1.
Thirty Herbarium Specimens Housed at University of California, Santa Barbara (UCSB), Collected From 1985 to 2001. In some cases, Wayne R. Ferren used the same collector ID for multiple individual specimens (noted by *). Each sequence used in this study was from a unique individual. Specimen ID is a unique identifier assigned during DNA extraction. **Specimen ID 994 was assembled and used as a reference sequence for remaining sequences. The Herbarium Catalog Number is equivalent to an herbarium accession number.

Appendix 2
Appendix2.
Mantel Tests of Genetic Distance as a Function of Geographic Distance or Presence of a Barrier. Each point represents a pairwise value of genetic and geographic distance between two individuals. Genetic distance (y-axis on each graph) is the percentage (%) of the number of sites on the nrDNA sequence that differ between each pair of individuals. Overlapping points are not distinguishable from one another. A. Isolation by Euclidean distance (km), the straight-line distance between localities (accounting for the curvature of the Earth). B. Isolation by oceanic distance (km), shortest oceanic distance between two localities. C. Isolation by shoreline distance (km), shortest shoreline distance between two localities.
