Invasion genetics of the mummichog (Fundulus heteroclitus): recent anthropogenic introduction in Iberia

Sampling

A total of 248 Fundulus heteroclitus individuals from 16 locations: 13 sites in the western Atlantic, one in the Mediterranean Sea and two locations in the eastern Atlantic (Table 1 and Figs. 1A, 1B) were obtained and stored in 96% ethanol and kept at −20 °C. The populations from the western Atlantic sampled above 40°N are hereafter referred to as northern locations (ID 1-8), while samples collected below that latitude are referred to as southern locations (ID 9-13). Samples collected from the Mediterranean and eastern Atlantic are referred to as Iberian (ID 14-16).

DNA extraction, PCR amplification and sequencing

Total genomic DNA was extracted from caudal fin tissue following a standard Chelex 100 protocol (Walsh, Metzger & Higuchi, 1991). Extraction results were checked by electrophoresis in 0.8% agarose gel stained with GelRed. Polymerase Chain Reactions (PCR) were conducted in a total volume of 25 µL, with 1X buffer, 10 mM dNTPs, 10mM of each primer, 1U Taq Advantage 2 Polymerase mix DNA polymerase (CLONTECH-TaKaRa), 2 µL of DNA and Milli-Q water to the final volume. A fragment of the cytochrome b (cyt b) gene (1,000 base pairs) was amplified with the forward primer GludG-L14724 (Palumbi et al., 1991) and the reverse primer cb6b.h (Martin & Bermingham, 1998). PCR amplification consists of an initial 4 min denaturation step at 95 °C, followed by 40 cycles of 1 min at 94 °C (denaturation), 1 min at 50  °C (annealing) and 1.5 min at 72 °C (extension), and a 5 min final extension step. When amplification was not successful the following profile was used: initial 3 min denaturation step at 95 °C, followed by 32 cycles of denaturation for 30 s at 95 °C, annealing for 30 s at 54 °C and extension for 1 min at 68  °C, and a final extension step for 4 min at 68  °C. PCR products were checked afterwards by electrophoresis in a 1% agarose gel stained with GelRed. Mitochondrial DNA was purified by ethanol/sodium acetate precipitation (Sambrook & Russel, 2001). Its purity and quantity were analysed using a NanoDrop1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Sequencing was performed on an ABI 3130xl capillary sequencer (Applied Biosystems –CCMAR, Portugal) using the forward primer from the PCR amplification (GludG-L14724).

Genetic analysis

Cyt b sequences were aligned and manually checked using the software Geneious v4.8.2 (Biomatters, Ltd., Auckland, New Zealand). The number of haplotypes (n), number of private haplotypes (n_p), and the haplotype (h) (Nei & Tajima, 1981) and nucleotide diversities (π) (Nei, 1987) were calculated for each location using the DnaSP v5.10.1 (Librado & Rozas, 2009). All sequences were uploaded to GenBank (MH809691–MH809938). To represent the phylogeographic relationships among haplotypes, a haplotype network was constructed using the Median Joining algorithm implemented in NETWORK v5.0 (Bandelt, Forster & Röhl, 1999, fluxus-engineering.com).

We used two approaches to infer the most probable source area of the introduced Iberian populations within the native range of F. heteroclitus: (1) the geographical distribution of haplotypes in native populations, and (2) the phylogeographical relationships among haplotypes.

Ecological niche modelling

To explore the potential spread of F. heteroclitus within European waters beyond its currently known locations we used an ecological niche model. We selected the interactive AquaMaps system in FISHBASE (http://www.fishbase.org). Using observed locations of adults to construct a range of acceptable environmental conditions within which a species can exist, this model uses a c-squares distribution modelling approach (Rees, 2003) to predict the probability of occurrence of the adult fish for a particular location, possessing specific environmental parameters. While this is normally computed just within the native range of a species it can be extended within the interactive system to consider all possible locations, and the associated probability of occurrence linked to the location’s mean environmental state if the species was able to reach that region. The basic statistical approach is given in Kaschner et al. (2006) and the AquaMaps version is described in Kesner-Reyes et al. (2012). Another ecological niche model, based on maximum entropy principles (MAXENT, Phillips & Dudík, 2008), was also tested using environmental parameters of sea surface temperature (SST) and sea surface salinity (SSS), but the specialist ecological substrate niche of F. heteroclitus led to poor solutions for range prediction using ENM. This latter approach is therefore not considered further here.

A set of favourable and extreme environmental conditions compatible with the native occurrence of F. heteroclitus is automatically specified by AquaMaps, relying heavily on Page & Burr (2011). It includes parameter ranges for water depth, SST, SSS, primary production rates and sea-ice cover, leading to probabilities of occurrence exceeding 0.6 along the whole eastern seaboard of North America from South Carolina to the Canadian Maritime Provinces south of the Gulf of St. Lawrence. This matches well to the known range (Page & Burr, 2011). However, the long-term survival of F. heteroclitus in two Iberian estuaries has implications for extending the extreme salinity level that this species can tolerate beyond that automatically specified. The extreme salinity tolerance has therefore been raised to 39.6 psu, compared to the automatic 36.47 psu initially prescribed by AquaMaps, consistent with the values found for the Ria Formosa (Cristina et al., 2016). In addition, the sea-ice variable has been excluded from constraining the AquaMaps solution, as this is not relevant to the area being considered in this study, and also the primary production variable, as coastal environments satisfy the automatic requirement almost uniformly. The environmental variables used for the AquaMaps simulation shown here are given in Table 2.

AquaMaps was re-run with these environmental constraints and using the combined occurrence data from the native range, along the North American eastern seaboard, and the two Iberian sites of Ria Formosa and the Ebro River Delta. There is no change to the North American predictions, so here we concentrate only on those for European waters. The occurrence of F. heteroclitus was further constrained to coastal zones with the muddy benthic habitats in which F. heteroclitus is found. These are defined as the coastal fine muds to muddy sands benthic environments of the EMODnet Seabed Habitats project (http://www.emodnet-seabedhabitats.eu). These environments are formally defined by Long (2006), but basically include sediments with more than a 10% mud component.

Genetic diversity

The presence of a single haplotype common to all F. heteroclitus sampled in the Iberian Peninsula lends support to the hypothesis of an extremely recent introduction of the species which has not allowed the accumulation of mutations at the mtDNA level, and with a single introduction event composed by a very small number of individuals (Roman & Darling, 2007). Theoretically, invasive species are expected to suffer loss of genetic variation since the new established populations are often based on a few individuals, which by definition, have lower genetic diversity than the native source populations (Dlugosch & Parker, 2008). The single-haplotype characteristic can be found in other invasive species, such as Equulites elongatus, the slender pony fish (Sakinan, Karahan & Ok, 2017); Cercopagis pengoi, the fishhook waterflea, a planktonic cladoceran crustacean (Cristescu et al., 2001); Corbicula fluminea, the Asian clam (Gomes et al., 2016) and Didemnum perlucidum, a sea squirt (Dias et al., 2016). However, many successful invasive species do not display significant erosion of genetic diversity (Dlugosch & Parker, 2008). For example, in a recent review of the literature on European sea invasion genetics, in 54% of studies that compared the genetic diversity between introduced species and their native range, 74% reported comparable levels of diversity between them, while only 23% displayed a reduction in the genetic diversity of introduced species, and the remaining 2% showed an increase in diversity (Rius et al.,  2014).

Population sources

We identified a single Iberian haplotype present in all northern populations in high frequency (between 47 and 93%) and in two individuals in Suffolk, one of the southern locations. According to our results, the northern group is the most probable source of the founder individuals, which corresponds to the natural range of the subspecies F. h. macrolepidotus. However, we cannot definitely exclude Suffolk as a presumptive population source. The absence of genetic diversity in the Iberian Peninsula populations prevents the precise determination of the putative source population. The low spatial resolution of our data arises from the lack of genetic variability in the invasive range, rather than insufficient sampling of F. heteroclitus individuals in the Iberian Peninsula or in its native range (Muirhead et al., 2008).

Nevertheless, our findings are consistent with two previous studies on the origin of invasive individuals found in the Guadalquivir. First, based on mtDNA restriction fragment length polymorphisms (RFLP) (Fernández-Pedrosa, Latorre & González, 1996) reported the presence of two haplotypes: the most abundant corresponding to the northern haplotype 1, dominant between Maine and Nova Scotia in North America; the other haplotype did not match any of the sampled native haplotypes and we found no evidence of its presence in the present study, using a larger number of individuals. It was previously suggested to be either a native unsampled haplotype or an endemic haplotype from Iberia, which is rather unlikely due to its recent invasion (Fernández-Pedrosa, Latorre & González, 1996) and total absence of records in the area. Although there are no reports of hybridization between F. heteroclitus and any of the Spanish endemic species, the presence of a new haplotype could nonetheless be due to hybridization (Rius et al., 2014). Secondly, a study based on cyt b sequences (Bernardi et al., 1995) concluded the individuals from the Guadalquivir originated in the region between Maine and Nova Scotia. While we cannot discount the possibility of a few rogues from Suffolk being responsible for the invasion, the weight of probability falls heavily on the side of Maine to Nova Scotia origin.

Introduction vector

Since it was first recorded in Iberian saltmarshes, several studies linked F. heteroclitus introduction to different vectors. The aquarium trade has been suggested as the most important vector responsible for the introduction of this species (e.g., Gozlan, 2010 and references therein), followed by ballast water (Fernández-Delgado, 2010; García-Revillo & Fernández-Delgado, 2009), biological control (Gozlan, 2010) and unknown origins (Fernández-Delgado, 1989; Gutiérrez-Estrada et al., 1998). We posit that the introduction via the aquarium trade or animal acquisition for scientific purposes is the most probable scenario responsible for the establishment of the first individuals in Iberia, given the genetic results obtained.

It is well known that vessels can transport large numbers of organisms from several species at the same time in their ballast water (Gollasch, 2007). Since there are countless vessels active around the world (e.g., Kaluza et al., 2010), ballast water-mediated transport allows the possibility of multiple introduction events, each with large groups of individuals (Hulme, 2009). Considering such large groups usually comprise higher genetic diversity than fewer individuals alone, the assemblages transported are likely to display levels of genetic diversity similar to the levels found within their native range (Wilson et al., 2009). Our findings, however, suggest that this was not the introduction vector responsible for the spread of the species; in contrast to the expectation of similar levels of genetic diversity between the native and invasive range, the Iberian locations display a strong founder effect, with all samples sharing a single haplotype.

Similarly to ballast water, the aquarium and ornamentals trade transport many species at a global scale (Padilla & Williams, 2004). For instance, at least 19% of the invasive fishes found in the Iberian Peninsula were introduced via the aquarium trade (Maceda-Veiga et al., 2013). However, there are two main differences between introductions that follow ballast water or the aquarium trade. First, each introduction event after aquarium release is likely to comprise a small number of individuals (Duggan, Rixon & MacIsaac, 2006). Thus, even though this vector may be responsible for the establishment of several non-native species (Padilla & Williams, 2004), multiple introductions would be necessary for an invasive species to display high genetic diversity (Roman & Darling, 2007). Secondly, the individuals released by aquarists are usually adults of higher fitness, which makes them better adapted to survive in a natural environment (Padilla & Williams, 2004). Thus, not only are F. heteroclitus’s invasive genetic diversity and structure consistent with an introduction of a low number of individuals via the aquarium trade, it is also plausible that a few resistant individuals would manage to survive, reproduce and colonize the environment in which they were released.

Human mediated introduction in the Ebro Delta

Although the absence of genetic structure within Iberia limits possible insights into the invasion pathway, our data support the hypothesis of a human-mediated introduction episode responsible for the establishment of F. heteroclitus in the Ebro Delta, as previously suggested by Gisbert & López (2007) based on taxonomic identification. The hypothesis of a long-distance colonization via natural dispersal is quite unlikely, as is strongly shown by the ENM analysis (Fig. 4).

According to a review of the geographical distribution of Cyprinodontiformes along the northeastern coast of Spain by García-Berthou & Moreno-Amich (1991), no F. heteroclitus individuals were found at the Ebro Delta in 1989; the first record of the species was only registered 16 years later by Gisbert & López (2007). Thus, we estimate the date of establishment in the Ebro Delta ranges between 12 and 28 years ago. Assuming this estimate is correct, a natural colonization hypothesis implies individuals would have taken roughly two decades to travel more than 1,000 km from their southernmost limit located in the Guadalquivir saltmarshes (Gutiérrez-Estrada et al., 1998) to the Ebro Delta. However, when compared with the natural colonization of the Ria Formosa, this hypothesis seems quite improbable. While F. heteroclitus was never collected during sampling events that happened in the Ria Formosa between 1980 and 2006 (França, Costa & Cabral, 2009; Ribeiro et al., 2006; Ribeiro et al., 2008), analysis of prey remains left by Little Terns (Sterna albifrons) in the salt-pans and barrier islands revealed this prey species was present in the salt-pans and adjacent channels at least since 2002 (Catry et al., 2006). Although this may sound contradictory, F. heteroclitus could in fact have been present in the Ria Formosa in specific unsampled locations or at extremely low densities, avoiding capture. Nonetheless, assuming that colonization happened around 2002, it seems that F. heteroclitus took no more than ca. 20 years to travel around 50 km from the Guadiana Delta, where it was first detected in 1976 (Coelho, Gomes & Ré, 1976)}, despite both areas being on a stretch of coast with a favourable benthic habitat (Fig. 4). This estimate indicates that if the Ebro colonization happened via natural dispersal, it must have happened 20 times faster than the natural colonization of the Ria Formosa. Given that F. heteroclitus has very low dispersal abilities (e.g., Fritz, Meredith & Lotrich, 1975; Lotrich, 1975), the natural colonization scenario for the Ebro seems unlikely.

Furthermore, if establishment in the Ebro Delta followed a natural range expansion one would expect to find several established populations between the Barbate marshes in the Guadalquivir, its southernmost location (Gutiérrez-Estrada et al., 1998), and the Ebro Delta. However, no invasive individuals were found at Mar Menor, a coastal lagoon with suitable conditions for successful establishment (Pérez-Ruzafa et al., 2006), located midway between the Guadalquivir and the Ebro Delta.

Finally, the Strait of Gibraltar would represent a barrier to dispersal and gene flow, limiting F. heteroclitus natural range expansion towards the Mediterranean Sea (Doadrio, Carmona & Fernandez-Delgado, 2002) because of the strong currents prevailing in the area and the absence of suitable habitats. Although uncertain for F. heteroclitus, this movement has already been documented for two other Iberian toothcarps, where restricted gene flow in this region led to speciation of the Aphanius iberus in the Mediterranean Iberian coast and Aphanius baeticus in the southwestern Atlantic Spanish coast (Doadrio, Carmona & Fernandez-Delgado, 2002).

Ecological niche modelling

Conditions for the spread of F. heteroclitus are limited by the existence of benthic muddy saltmarsh environments compatible with the species’ ecological requirements. These habitats are only found near major estuaries or lagoons areas along the Atlantic and western Mediterranean coastlines. The exception to this is along the North African coast (Fig. 4), where such environments are more continuous, although other aspects may be less favourable there, as suggested by the rapid decline in occurrence probability eastwards. The consequence of this benthic habitat constraint is to make along-shore colonisation unlikely in most areas, suggesting that any such sudden expansion would need to be aided by human release.