Population genetics of Bursatella leachii (De Blainville, 1817) and implications for the origin of the Mediterranean population

The sea hare Bursatella leachii (de Blainville 1817) (Mollusca: Gastropoda: Heterobranchia: Aplysiidae) is a pantropical sea slug that has colonized the Mediterranean Sea in modern times. Because the initial records in the non-native range started in the Eastern Mediterranean, and its pattern of spread was relatively consistent with those observed in well-known Lessepsian invaders, B. leachii is commonly considered to be a migrant from the Red Sea. In this study, we investigate for the first time the origin of the Mediterranean populations of B. leachii inferring their population structure and assessing relatedness levels of different regional populations. Sequence data from the cytochrome oxidase I were used to conduct population genetic analyses on this species, particularly by investigating the genetic structure of Atlantic, Mediterranean, and Indo-Pacific populations. Our results confirm that B. leachii is a truly pantropical species which displays geographic structure among major ocean basins. More importantly, sequenced Mediterranean and Atlantic animals share similar or identical haplotypes, which are distinct in at least 5 substitutions from haplotypes recovered from Indo-Pacific specimens. The results suggest that the Mediterranean population of B. leachii here examined probably have an Atlantic origin, and for the first time casts doubts on the assumed primary pathway of migration into the Mediterranean Sea.

The Mediterranean Sea constitutes a particularly interesting system to study marine biological invasions. Because it is a relatively closed basin, connected to two markedly distinct biogeographic regions (the Atlantic Ocean and the Red Sea), the Mediterranean is susceptible to invasion by organisms adapted to different environmental conditions and having very different evolutionary histories. Marine alien species can potentially enter and settle in the Mediterranean through its three aquatic open access points (the Gibraltar Strait, the Suez Canal, and the Bosporus) by a wide range of different vectors [1]. Some species have been introduced by direct human activities such as aquarium trade, shipping, or aquaculture, which constitute the most common vectors of introduction of alien species in marine environments worldwide [2]. Others, the so called Lessepsian or Eritrean migrants [3], have arrived into the Mediterranean by “natural” dispersal during their larval stage from the Red Sea [1]. This is the most common pathway for Mediterranean invasions, and it is likely to become more prevalent in the future because of the recent expansion of the Suez Canal [3]. Moreover, the “tropicalization” of the Mediterranean [4] has resulted in species poleward distributional shifts, and has favored the arrival of tropical species to temperate regions [5]. Also, these shifts in environmental and ecological conditions have facilitated the successful establishment of species introduced by direct human activities in the Mediterranean [6]. To add to this complexity, the weakening of natural barriers for dispersal from the Atlantic Ocean into the Mediterranean, such as the Canary Current [7] and the Almería-Oran Front [8], caused by climate change, have facilitated the recent migration of tropical Atlantic species into the Mediterranean [9].

The ragged sea hare Bursatella leachii (de Blainville, 1817) (Gastropoda: Heterobranchia: Aplysiidae) is a benthic detritivore typically found grazing on the surface layers of muddy and sandy bottoms of protected bays and estuaries; it is distributed worldwide in tropical and warm temperate waters [10,11,12,13,14,15] and present in the Mediterranean as an established alien species [16, 17]. Because apparently B. leachii colonized the Mediterranean starting in the mid-20th century in an east-to-west direction (Fig. 1), almost all authors working on the Mediterranean biota have hypothesized a Red Sea origin of the local populations, categorizing B. leachii as a Lessepsian migrant [e.g. 17–21, among recent reviews]. However, molecular studies have never been used to more accurately determine the origin of these populations. In this study, we used a population genetics approach through mtDNA sequence data to infer the population structure of B. leachii for the first time across its entire range and to determine the genetic affinities of the Mediterranean population, therefore testing the origin of this invasion event.

Records of Bursatella leachii in the Mediterranean Sea, adapted from [17, 41]. Dots indicate the location of each sighting, reference list in Selfati et al. [17]. The decades in which the first sighting occurred for each location are coded by color (black: 1950’s and 60’s; dark blue: 1970’s; light blue: 1980’s; green: 1990’s; orange: 2000 to present)

Full size image

Specimens

A total of 33 individuals from various localities across the global range of B. leachii were studied (Fig. 2). The samples were preserved in 70–90% ethanol; some consisted of full specimens, while others were small pieces of tissue cut from either the foot or the mantle of the animal. The specimens were collected in the field, mainly by scuba diving and snorkeling at low tide or obtained from the Cal Poly Pomona Invertebrate Collection (CPIC) and the Invertebrate Zoology collection of the California Academy of Sciences (CASIZ). A full list of specimens with locality and collection date is available in Table 1. The list also comprises the GenBank accession numbers for all of the sequences as well as those obtained from GenBank used in the study.

Map showing localities of specimens sequenced for this study. Multiple specimens from the same locality are merged into a single point

Full size image

DNA extraction, amplification, sequencing

A small piece cut from the tissue sample or from the foot or mantle of the full specimen was fragmented using sterile forceps and razor blades. DNA of more recently collected specimens was extracted using the Hot Chelex^® extraction method. For older specimens, the DNeasy Blood and Tissue Kit (Qiagen) was used. The Hot Chelex^® extraction was conducted as follows. After maceration, the tissue was transferred in 1 mL of 1X TE Buffer, contained in a 1.7 mL microcentrifuge tube. To allow the tissue to rehydrate and the cells to lyse, the tube was then placed on a rotating block for at least 20 min. After this time, the tube was vortexed for 5 s and then centrifuged for 3 min at 21,000g. Careful not to disturb the tissue pellet, 975 µL of 1X TE Buffer were removed from the tube, and 175 µL of 10% Chelex^® solution were then added to the tube. After being vortexed, the tube was put in a 56 °C hot water bath for at least 20 min. After 5 s of vortexing, the tube was then put in a 100 °C heat block for exactly 8 min. Finally, the tube was vortexed and centrifuged for 3 min at 21,000×g. The supernatant resulting from the centrifugation, was then used as a template for amplification through Polymerase Chain Reaction (PCR). DNeasy extraction was conducted following the manufacturers’ protocol for tissue samples.

A portion of the cytochrome c oxidase 1 (COI) mitochondrial gene was amplified by PCR. For the amplification of COI, both the universal primers (LCOI490 5′-GGTCAACAAATCATAAAGATATTGG-3′, HCO2198 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) identified by Folmer et al. [22], and specific primers designed for B. leachii (Forward 5′-GGGATGTGATGTGGTCTAGTAGGAAC-3′, Reverse 5′-ATAGGATCTCCTCCTCCCGCAGGATC-3′), were used in this study. Each reaction tube contained a total of 50 µL of reaction mixture, which included: 0.25 µL of 5 u/µL taq polymerase, 5.00 µL of 10× PCR buffer, 5.00 µL of 25 mM MgCl₂, 1.00 µL of 40 mM dNTPs, 1.00 µL of each 10 mM primer, 34.75 µL of H₂O and 2.00 µL of template DNA. The thermocycler was set as follows: 3 min at 95 °C for an initial denaturation of the DNA; 35 cycles of denaturation at 94 °C for 45 s, annealing at 45 °C for 45 s and elongation at 72 °C for 90 s; a final elongation step of 10 min at 72 °C.

A 10% agarose gel was prepared and stained with ethidium bromide, and a gel electrophoresis of the PCR products was run, in order to verify the presence of DNA fragments. PCR products that showed appropriate size fragments were then purified using the GeneJET PCR Purification Kit by Thermo Scientific, following the manufacturers’ protocol. Concentrations of the purified DNA fragments were then checked using a NanoDrop 1000 spectrophotometer (Thermo Scientific) and sent to the Source Biosciences (San Diego, CA) company, which performed Sanger sequencing of the fragments. Assemblage and editing of the sequences was done with the software GENEIOUS PRO 4.7.4 [Biomatters Ltd., 23]. Consensus sequences were aligned and concatenated for the following analyses using the same software.

Population genetics

Population genetics analyses were conducted using partial sequences of the cytochrome c oxidase subunit I (COI) mitochondrial gene. In order to visualize and interpret genetic data at a population-scale level, a haplotype network was constructed using the TCS method [24] implemented in the software PopART (http://popart.otago.ac.nz). In order to assess the relatedness level of the different haplotypes, a matrix of pairwise genetic distances was constructed using the software MEGA 7.0.26 [25], using the Kimura 2-parameter model [26]. Distance values were plotted in a histogram.

The haplotype network analysis reveals there is geographic structure in B. leachii populations (Fig. 3), with two distinct haplogroups separated by four mutational changes. One includes haplotypes from the West Pacific, India, Japan, and South Africa. The other haplogroup includes samples from the Caribbean and the Mediterranean. In the Atlanto-Mediterranean haplogroup, the most common haplotype is shared among Atlantic and Mediterranean specimens. Most of other Mediterranean and Caribbean haplotypes only differ from that common haplotype by one nucleotide substitution, and the most distant ones differ by three mutational changes.

COI haplotype network produced with the TCS method in PopART. Geographic region in which each haplotype was found is coded by color. The relative size of the circles is proportional to the number of sequences of that same haplotype. For haplotypes found in only one specimen, the isolate and location is indicated. The groups of specimens sharing the same haplotype are as follows: H1 (isolates: GB03, GB04, EN016, EN203, EN205, KP25, KP27, KP59, BURIC1, HG88), H2 (isolates: EN201, EN202, EN213, HG124, HG125, HG127), H3 (isolates: HG121, HG123, HG126) and H4 (isolates: EN105, EN113, EN115, HG112, HG113)

Full size image

A histogram of the pairwise genetic distances is shown in Fig. 4. Most of the haplotypes are either identical to one another or very similar, with low genetic distance among them. Higher genetic distances among haplotypes are less common. There appears to be a linear inverse relationship between distance values and their frequency, without gaps in intermediate distances.

Histogram of pairwise genetic distances among COI haplotypes calculated with the Kimura-2 model in MEGA

Full size image

The main results from the present paper revealed that although Pacific and Atlantic populations of B. leachii are genetically distinct, the divergence between populations from these two geographic regions is relatively small (5 substitutions), suggesting B. leachii is likely to be a true pantropical species, experiencing relatively recent gene flow across ocean basins. This is supported by the results of the genetic distance calculations (Fig. 4), which show that no barcoding gap is present between the frequency of higher distances (intraspecific) and lower distances (interspecific). The finding of a single species contrasts with the majority of similar examples studied to date. In fact, most pantropical species of sea slug have not survived the advent of molecular systematics [see 27–32], with the sole exception of the pelagic aeolid nudibranch Glaucus atlanticus (Forster, 1777). A possible difference between B. leachii and G. atlanticus with the other species, which could explain their unique range, is that both species also occur in cold and temperate waters in South Africa, and this region could serve as a gateway that allows the maintenance of gene flow between the South Atlantic and the Indian Ocean [33]. Studies on the reproduction and population dynamics of B. leachii conducted by Paige [34] and Clarke [35] may provide additional support for such a broad range and dispersal ability. The planktonic larvae of B. leachii can survive for more than 3 months after hatching, and adults have relatively short life cycles, reaching maturity in only 2–3 months [34]. At the same time, when compared to other sea hares, B. leachii is able to reach relatively large population sizes even when lacking access to optimal food sources [35], a life history trait advantageous for a species colonizing a new territory [36].

Another result obtained from this study, is that no shared haplotypes between the Indo-Pacific region and the Mediterranean were found, despite of a relatively large Mediterranean sample size and sampling effort covering both the eastern and western Mediterranean. In addition to the membership to the same haplogroup, the existence of identical haplotypes strongly suggest a shared demographic history between populations of B. leachii of the Atlantic and Mediterranean regions. On the contrary, all Indo-Pacific haplotypes sequenced for this study, including specimens from the Indian Ocean and South Africa, are genetically distinct by at least 5 substitutions from all of the Atlantic and Mediterranean haplotypes. Therefore, although most authors have considered the Mediterranean populations of B. leachii to have an Indo-Pacific origin (see above), the data presented here suggest otherwise, that at least a part of the B. leachii populations in the Mediterranean Sea have migrated from the Atlantic Ocean. One important caveat is that for this study we had no access to specimens from the Red Sea, which prevent us from rejecting with certainty a Lessepsian origin. However, a Red Sea origin hypothesis would require the existence of Red Sea populations genetically similar to those present in the Caribbean and Mediterranean and different from those in other Indo-Pacific regions studied herein, including other Indian Ocean populations. If such population were found, it would suggest an anti-Lessepsian migration into the Red Sea, rather than a Red Sea origin of B. leachii.

The data here presented opens the possibility of an alternative pattern of colonization of the Mediterranean, but the mechanisms responsible for the invasion remain unclear. One hypothesis is that the migration into the Mediterranean could be human-mediated. Transport by ship ballast water has been shown to be one of the most common vectors by which alien marine species can disperse into new habitats [2]. However, B. leachii has planktotrophic larvae, which typically encounter sub-optimal growth conditions while traveling in ballast water tanks [37], reducing their ability to survive long journeys. Moreover, the spatial patterns of dispersal shown during the colonization of the Mediterranean are inconsistent with a ballast water mediated introduction event. For example, there is no clear link between the first sightings and proximity to the main Mediterranean harbors, where ballast water should had been discharged. An alternative hypothesis is that the colonization of Mediterranean could be the result of non-human mediated dispersal of B. leachii through the Gibraltar Strait. Although the Mediterranean Sea has been connected with the Atlantic Ocean through the Gibraltar Strait for more than 5 million years [8], recent oceanographic changes caused by global warming have altered the permeability of barriers to dispersal. For example, the weakening of the cold Canary Current and Saharan upwelling [7], which historically prevented the dispersal of tropical West African species into the Mediterranean, might have resulted in milder winters throughout the last century [38]. This should have facilitated the spread of species, which maintain populations in West Africa, like B. leachii. Nevertheless, the east to west direction of colonization observed in B. leachii appears to contradict this hypothesis. At least two explanations can be advanced to justify this counter-intuitive pattern. A similar pattern has been locally observed in another tropical sea hare, Aplysia dactylomela (Rang 1828). Valdés et al. [39] suggested that larvae of A. dactylomela, upon entering in the Mediterranean, were transported by the Algerian current directly into the Central and Eastern Mediterranean, where the first records were registered, therefore bypassing the Western Mediterranean. The observed pattern of Mediterranean colonization of B. leachii is not identical to the that observed for A. dactylomela, but it is comparable, and a similar case could be made to explain the first reports of B. leachii in the Eastern Mediterranean. Another possible explanation is that B. leachii, after entering the Western Mediterranean from the Atlantic, could have migrated towards the Eastern Mediterranean by stepping-stone dispersal. Permanent populations might not have established on the way to the Eastern Mediterranean most probably because of sub-optimal temperature conditions for survival, hence no sightings in that area have been registered. More recently, the warming of the Mediterranean has resulted in more favorable conditions for settlement of warm-water species [4, 6]. Differences in winter sea surface temperatures have been hypothesized as creating a division between distributional ranges of cold and warm water biota in the Mediterranean [35]. The February sea surface temperature isotherm of 15 °C, in particular, represents the possible northern limit to the distributional range of different warm water species [40]. With temperatures increasing during the most recent decades [38], these “barriers” have also shifted [40], and B. leachii might have found more suitable conditions for establishment and colonization of the Western Mediterranean (see Fig. 1).

The results of this study suggest that B. leachii is likely to be a true pantropical species, and provide for the first time evidence that at least part of the Mediterranean population of B. leachii has originated in the Atlantic Ocean, contrary to what was previously inferred purely from the sequence of sightings. Different hypotheses regarding the mechanism behind the migration into the Mediterranean were formulated, considering both man-made and natural dispersal possibilities. The uncertainty of B. leachii native range and relatively small sample size of the study, limit our ability to test these additional hypotheses.

EB extracted and amplified DNA of part of the specimens, conducted sequence alignment and analysis (haplotype network), prepared Figs. 1, 3 and 4, and was a major contributor in writing the manuscript. FC co-supervised the project, helped prepare Fig. 2, and was a major contributor in writing the manuscript. KEP and HG extracted and amplified DNA of part of the specimens. HB, SG, AJ, SC and JLCC contributed to specimen acquisition, and corrected and improved the final manuscript. DP contributed to specimen acquisition, prepared Fig. 2, and corrected and improved the final manuscript. AV supervised the project, contributed to specimen acquisition, conducted sequence alignment and analysis (Φ_ST pairwise analysis), and was a major contributor in writing the manuscript. All authors read and approved the final manuscript.

Acknowledgements

Hessam Ghanimi, Jennifer McCarthy, and Sabrina Medrano (California State Polytechnic University, USA) greatly helped with the laboratory work. Elena Tricarico (University of Florence) helped with the writing of part of the manuscript and supervised EB’s work during his Master thesis. Biju Kumar (University of Kerala, India) assisted with the specimens from India. Liz Kools and Terry Gosliner facilitated the loan from the CASIZ. Dr R. Ravinesh greatly assisted during field collection in India. An anonymous reviewer provided precious insight and suggestions that greatly helped to improve the manuscript.

Competing interests

The authors declare that they have no competing interests.

Availability of data and materials

All sequences analyzed in this study are available in GenBank (https://www.ncbi.nlm.nih.gov/genbank/). Accession number for each sequence can be found in Table 1.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Funding

The TROPIMUNDO Steering Board provided funding for traveling expenses through the TROPIMUNDO Thesis Grant. The Department of Science and Technology (Government of India) provided funding through the INSPIRE fellowship.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Affiliations

1. Department of Biological Sciences, California State Polytechnic University, 3801 West Temple Avenue, Pomona, CA, 91768, USA

   • Enrico Bazzicalupo
   • , Ka‘ala Estores-Pacheco
   • , Haleh Golestani
   •  & Ángel Valdés

2. Department of Integrative Marine Ecology, Stazione Zoologica Anton Dohrn, Villa Comunale, 80121, Naples, Italy

   • Fabio Crocetta

3. BioBio Research Center, Faculty of Sciences, Mohammed V University in Rabat, 4 Avenue Ibn Battouta, B.P. 1014 RP, Rabat, Morocco

   • Hocein Bazairi

4. ChiBioFarAm Department, Messina University, Viale Ferdinando Stagno d’Alcontres 31, 98166, S. Agata-Messina, Italy

   • Salvatore Giacobbe

5. Center for Marine Research, Ruđer Bošković Institute, G. Paliaga 5, 52210, Rovinj, Croatia

   • Andrej Jaklin

6. Foundation for Research and Technology - Hellas (FORTH), Institute of Applied and Computational Mathematics, N. Plastira 100, Vassilika Vouton, 70013, Heraklion, Greece

   • Dimitris Poursanidis

7. Department of Aquatic Biology and Fisheries, University of Kerala, Thiruvananthapuram, Kerala, 695581, India

   • B. K. Sneha Chandran

8. Departamento de Biología, Facultad de Ciencias del Mar y Ambientales, Campus de Excelencia Internacional del Mar (CEIMAR), Universidad de Cádiz, Avenida República Saharaui, s/n, Ap.40, 11510, Puerto Real, Cádiz, Spain

   • Juan Lucas Cervera

9. Instituto Universitario de Investigación Marina (INMAR), Campus de Excelencia Internacional del Mar (CEIMAR), Universidad de Cádiz, Av. República Saharaui, s/n, Ap.40, 11510, Puerto Real, Cádiz, Spain

   • Juan Lucas Cervera

Corresponding authors

Correspondence to Enrico Bazzicalupo or Ángel Valdés.

Reprints and Permissions