Phylogeography and Diversity Among Populations of the Toxigenic Benthic Dinoflagellate Prorocentrum From Coastal Reef Systems in Mexico

Cembella, Allan D.; Durán-Riveroll, Lorena M.; Tarazona-Janampa, Ulrike I.; Okolodkov, Yuri B.; García-Sandoval, Ricardo; Krock, Bernd; Hörstmann, Cora; John, Uwe

doi:10.3389/fmars.2021.716669

ORIGINAL RESEARCH article

Front. Mar. Sci., 15 October 2021

Sec. Global Change and the Future Ocean

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.716669

This article is part of the Research Topic Harmful Algal Blooms in Marine Environments: From Biologically Active Compounds to Species Diversity View all 6 articles

Phylogeography and Diversity Among Populations of the Toxigenic Benthic Dinoflagellate Prorocentrum From Coastal Reef Systems in Mexico

$\r\nAllan D. Cembella*&#x;$ Allan D. Cembella^1*†

Lorena M. Durán-Riveroll^1,2*†

Ulrike I. Tarazona-Janampa^1,3,4

Yuri B. Okolodkov^5†

Ricardo García-Sandoval⁶

Bernd Krock¹

Cora Hörstmann^1†

Uwe John^1,7†

¹Alfred-Wegener-Institut, Helmholtz-Zentrum für Polar- und Meeresforschung, Bremerhaven, Germany
²CONACyT-Departamento de Biotecnología Marina, Centro de Investigación Científica y de Educación Superior de Ensenada, Ensenada, Mexico
³International Max Planck Research School, Max Planck Institute for Marine Microbiology, Bremen, Germany
⁴Laboratorio de Microbiología Molecular y Genómica Bacteriana, Universidad Científica del Sur, Lima, Peru
⁵Instituto de Ciencias Marinas y Pesquerías, Universidad Veracruzana, Boca del Río, Mexico
⁶Departamento de Biología Comparada, Facultad de Ciencias, Universidad Nacional Autónoma de México, Coyoacán, Mexico
⁷Helmholtz Institute for Functional Marine Biodiversity, Oldenburg, Germany

The marine dinoflagellate genus Prorocentrum Ehrenberg comprises many species occupying primarily benthic or epiphytic habitats, particularly in tropical and sub-tropical waters. Despite concerted efforts to establish phylogenetic associations, there remain unresolved issues in defining morphospecies and membership in species complexes. The study described herein addressed the inter- and infraspecific relationships of members of the Prorocentrum lima and Prorocentrum hoffmannianum species complexes (PLSC and PHSC, respectively) by applying multivariate approaches in morphotaxonomy, molecular phylogenetics and chemodiversity to establish affinities among multiple clonal isolates. Morphotaxonomic analysis showed consistency with classical morphospecies descriptors, and high variability in cell size and dimensions, but did not challenge current species complex concepts. Phylogenetic analysis of ITS/5.8S rDNA sequences from isolates from the Gulf of California, Caribbean Sea, and Gulf of Mexico coasts compared with archived global GenBank sequences served to define five consistent clades with separation of the PLSC and PHSC. Secondary structure modeling of ITS2 rRNA variation based on compensatory base changes (CBC) was effective in resolving details of the respective species complexes and even indicated putative incipient or cryptic speciation due to potential hybridization barriers. This study represents the largest (n = 67 isolates) chemodiversity analysis of polyketide-derived toxins associated with diarrheic shellfish poisoning (DSP) from a benthic dinoflagellate genus. Relative composition of some analogs (OA, OA-D8, DTX1, DTX1a, and DTX1a-D8), including two new undescribed isomers, distinguished P. lima from P. hoffmannianum sensu lato, but without clear associations with substrate type or geographical origin. Although all P. lima and most (one exception) P. hoffmannianum were toxigenic, the total cell toxin content could not be linked at the species level. This research demonstrates that clonal chemodiversity in toxin composition cannot yet be effectively applied to define ecological niches or species interactions within local assemblages. Phylogenetic analysis of the ITS/5.8 rDNA, particularly when combined with secondary structure modeling, rather than only a comparison of LSU rDNA sequences, is a more powerful approach to identify cryptic speciation and to resolve species complexes within benthic dinoflagellate groups.

Introduction

Benthic species of the dinoflagellate genus Prorocentrum Ehrenb. (Prorocentrales: Prorocentraceae) are globally distributed from polar latitudes to the tropics, but comprise a higher percentage of the epibenthic microflora in tropical and subtropical coastal waters (Durán-Riveroll et al., 2019). Although not all Prorocentrum species are toxigenic, the production of polyether toxins is a common trait among benthic species, such as Prorocentrum concavum Y. Fukuyo, Prorocentrum faustiae S. L. Morton, Prorocentrum hoffmannianum (syn. Prorocentrum belizeanum M. A. Faust), Prorocentrum lima (syn. Prorocentrum arenarium M. A. Faust), and Prorocentrum maculosum M. A. Faust, but these toxins are essentially absent or have not been confirmed from planktonic taxa. The polyether toxins found among Prorocentrum species include okadaic acid (OA) and at least two dozen dinophysistoxin (DTX) analogs, plus related polyketides of uncertain toxicity (Hu et al., 2010). The toxins from these benthic species can accumulate in suspension-feeding shellfish and thereby cause a complex syndrome known as diarrheic shellfish poisoning (DSP) in human consumers of contaminated shellfish.

Prorocentrum lima (Ehrenb.) F. Stein 1878 emend. Nagahama et al. (2011) is one of the most widely distributed species of the genus. In tropical and sub-tropical coastal ecosystems (Durán-Riveroll et al., 2019) the species is commonly found as an epiphyte upon macroalgae, seagrasses, and inanimate substrates (Tarazona-Janampa et al., 2020). There is limited and somewhat contradictory evidence of biotic substrate preferences for benthic Prorocentrum (Boisnoir et al., 2019), but little is known about the basis of the selective mechanisms, chemical ecological interactions or potential effect on toxin production or composition. Recent molecular genetic studies of P. lima have revealed a high genetic variability, apart from (and not always consistent with) great variability in some morphological features (cell shape, length and width, number and shape of marginal pores). Accordingly, these similar but unresolved taxa were assigned to the “P. lima complex” as proposed by Aligizaki (2009), and in subsequent literature [e.g., Moreira-González et al. (2019)]. Herein these taxa are referred to as members of the P. lima species complex (PLSC).

Prorocentrum species from coastal and reef systems of Mexico commonly include P. hoffmannianum M. A. Faust 1990, often coinciding with P. lima (Durán-Riveroll et al., 2019). As in the current study, specimens of P. lima and P. hoffmannianum sensu stricto are well distinguished morphologically within field populations and usually in clonal culture. Nevertheless, unresolved issues regarding the synonymy and affinities with other described morphospecies closely related to P. hoffmannianum remain under debate. Combined molecular genetic and morphological analysis for P. hoffmannianum and relatives such as P. belizeanum (Herrera-Sepúlveda et al., 2015) have demonstrated the existence of a P. hoffmannianum species complex (PHSC), but not resolved the association with the PLSC. The challenge addressed in the current study at the infraspecific level was to determine affinities within and between the PLSC and PHSC by multivariate application of morphological, molecular phylogenetic, and toxin chemodiversity techniques.

Prorocentrum rhathymum A. R. Loeblich et al. (1979) is morphologically distinct from both the PLSC and PHSC, but remains an enigma with uncertain species distinction, e.g., from Prorocentrum mexicanum B. F. Osorio (Steidinger and Tangen, 1996; Faust and Gulledge, 2002; Cortés-Altamirano and Sierra-Beltrán, 2003), based on both morphological and ecological (habitat) differences. Molecular phylogenetic analysis of P. rhathymum strains from Florida (An et al., 2010) vs. P. cf. rhathymum from Korea has served to highlight inconsistences in species level discrimination (Lim et al., 2013). The most recent molecular phylogenetic analysis of Prorocentrum from the eastern Caribbean Sea based on LSU rDNA sequences clusters P. mexicanum closely together P. rhathymum but basal to and distinct from the paraphyletic PLSC clade (Chomérat et al., 2019).

Benthic dinoflagellate species are quite well described morphologically from Latin America (reviewed in Durán-Riveroll et al. (2019) but there are many gaps in known phylogeographical distribution, few confirmed records of cell toxicity or toxin composition and a paucity of molecular confirmation of species identities and phylogenetic relationships. This poor state of knowledge also applies to benthic Prorocentrum species although members of the genus have been reported and described from both the Pacific [for references, see Okolodkov and Gárate-Lizárraga (2006); Gárate-Lizárraga et al. (2007)], and Gulf of Mexico and Caribbean coasts of Mexico (Okolodkov et al., 2007, 2014; Barón-Campis et al., 2014; Almazán-Becerril et al., 2015; Aguilar-Trujillo et al., 2017).

The work presented herein comprises a comparative analysis of cell morphotaxonomy, toxigenicity and molecular diversity of clonal isolates from populations of benthic Prorocentrum from three coastal marine ecosystems, representing the Gulf of California, Gulf of Mexico, and Caribbean Sea. The aim was to identify diversity patterns within and among geographical populations and substrate affinities and to define coherences and discrepancies (e.g., putative cryptic species) based on these different characteristics.

Materials and Methods

Isolation and Culture of Prorocentrum

Samples for benthic dinoflagellate isolation were collected from targeted locations along the Gulf of California, Gulf of Mexico, and Caribbean coasts of Mexico (Figure 1). Primary collection sites were within Bahía de La Paz (Baja California Sur (24°09′30.01″N 110°19′12.10″W), the Veracruz Reef System (VRS), Veracruz (19°11′54.10″N, 96° 4′0.70″W), and Puerto Morelos, Quintana Roo (20°50′48.55″N, 86°52′30.53″W) (Table 1). Epibenthic substrates were collected from seagrass beds (Thalassia testudinum Banks ex König) and from polypropylene ropes indicating swimming areas by snorkeling in shallow water along sandy shores. Specimens were also collected along rocky shores, from attached and floating macroalgae, Ulva (Chlorophyta: Ulvaceae), Laurencia (Rhodophyta: Rhodomelaceae), Sargassum and Padina spp. (Phaeophyceae: Sargassaceae and Dictyotaceae, respectively), and detritus on sediments and buoys.

FIGURE 1

Figure 1. Map of the sampling localities for Prorocentrum isolates from the Gulf of California, Gulf of Mexico, and Caribbean coasts of Mexico. (1) Bahía de La Paz, Baja California Sur. (2) Veracruz Reef System, Veracruz. (3) Puerto Morelos, Quintana Roo.

TABLE 1

Table 1. Cultured clonal isolates (n = 69) of Prorocentrum species from locations on the Gulf of California, Gulf of Mexico, and Caribbean coasts, indicating locality of origin and substrate type, and provisional species assignment according to morphological characteristics.

Live samples were transported with site water in 50 mL conical plastic centrifuge tubes in an ice chest with ice packs to maintain ambient temperature around 24°C during the transport to the laboratory at Instituto de Ciencias del Mar y Limnología, Universidad Nacional Autónoma de México (ICMyL-UNAM), Mexico City. Substrate specimens and surrounding seawater were examined for epibenthic dinoflagellates in Petri plates under a stereo-dissecting microscope (Discovery.V8, Zeiss, Göttingen, Germany). Substrates were gently brushed and single-cells of epibenthic dinoflagellates were isolated by micropipette into sterile 96-well microplates containing 300 μL 50%-strength GSe growth medium (Blackburn et al., 2001) modified without soil extract and prepared from autoclaved seawater filtered through sand, activated carbon and 1 μm-cartridge-filters. The growth medium, supplemented with GeO₂ (final concentration: 2.5 mg L^–1) to inhibit diatom growth, was prepared from heat-sterilized seawater stock at salinity 36. Clonal isolates were cultured by incubation at 25 ± 1°C on a 12:12 h light:dark cycle and illumination of 50 μmol photons m^–2 s^–1.

Once dinoflagellate cell division was observed (around 6–8 cells per well) after 10–12 days, preliminary cultures were transferred into 24-well microplates, each well containing 2 mL growth medium, and after subsequent growth, into 60 × 15 mm sterile plastic Petri plates with 15 mL 50%-strength GSe seawater medium. Clonal isolates were maintained at ICMyL-UNAM, Mexico under the same conditions for morphological species identification by microscopy.

Subcultures of the collection of Prorocentrum isolates (n = 69) were transferred to the culture facilities at the Alfred Wegener Institute (AWI), Bremerhaven, Germany, for production of cells for rDNA and toxin analysis. These stock cultures were initially grown in plastic Petri plates on 50%-strength GSe, and then acclimated to K medium (Keller et al., 1987) added as enrichment to 0.2 μm-filtered and autoclaved North Sea seawater at salinity 33. All reference and experimental cultures were maintained at 24 ± 1°C on a light:dark cycle of 14:10 h and illumination of 86 μmol photons m^–2 s^–1.

Exponentially growing Prorocentrum cells were harvested from Petri plates for rDNA analysis, and after cell enumeration by microscopy (target ca. 10⁶ cells), were collected by centrifugation at 3,220 × g for 15 min at 4°C (Eppendorf 5810R, Hamburg, Germany). After removal of the aqueous supernatant, the pellet was transferred to a 2 mL microcentrifuge tube for total DNA extraction from fresh pellets.

Harvest details of actively growing Prorocentrum cultures and cell enumeration for DSP toxin analysis followed previously described methods (Tarazona-Janampa et al., 2020). In brief, cultured cells were harvested from plastic Petri plates by multiple centrifugation steps. The final cell pellets in 2 mL FastPrep tubes were placed in a thermomixer (Eppendorf, Hamburg Germany) at 100°C for 5 min to inactivate esterase enzymes that could modify the toxin profile. Cell pellets were stored frozen (−20°C) for later extraction.

Morphological Description of Prorocentrum Isolates

Specimens of Prorocentrum isolates from early phase stable cultures were examined for morphological characteristics by both light and scanning electron microscopy (SEM). Teratological forms rather common in cultures were avoided. Provisional assignment to morphospecies was conducted by light microscopy (Motic BA310E, Hong Kong, SAR China) at 1000× magnification and from micrographs taken with a digital camera (6 MP, Moticam S6, Motic, Hong Kong, SAR China). Further details of the thecal plates were studied after staining cells with 0.2% Calcofluor White M2R in aqueous solution (Fritz and Triemer, 1985). Cells were observed under epifluorescence microscopy (Axio Scope.A1, Carl Zeiss, Oberkochen, Germany) with filter set 18 shift free EX BP 390–420 (excitation), BS FT 425 (optical divider), and EM LP 450 (emission), using the Plan-Neofluar 40×/0.75 and 63×/0.95 Korr objectives (total magnification 400× and 630×, respectively). Photomicrographs were taken with a Carl Zeiss Axiocam 506 color camera (6 MP) using the ZEN 2012 SP2 program (Carl Zeiss Microscopy GmbH, Göttingen, Germany).

Cell dimensions (n = 30) for morphometric analysis were determined from cells either live or freshly fixed in Lugol’s iodine or neutralized 4% formaldehyde and examined with a light microscope (Motic BA310E, Hong Kong, SAR China) at 400× magnification. Cell size was estimated from digital camera (6 MP, Moticam S6, Motic, Hong Kong, SAR China) images captured with Motic Images Advanced 3.0 software, previously calibrated with a micrometer.

Cultured dinoflagellate cells for SEM were fixed with 2% glutaraldehyde for 90 min. Specimens were washed in 1.5 mL distilled water (5°C) and centrifuged at 1,200 × g at 5°C for 6 min. The wash procedure was repeated four times. Samples (1.5 mL each) underwent a graded ethanol dehydration series (10, 30, 50, 70, 90, and 99%), with centrifugation as for the sample wash at each dehydration step. Finally, 250 μL of hexamethyldisilazane air-drying agent was placed upon the sample on aluminum SEM stubs, and specimens were gold sputter-coated for 5 min. Dinoflagellate cells were observed with a JEOL JSM6360LV SEM equipped with a backscattered electron detector under 8 kV voltage acceleration and at 15 mm working distance.

DNA Sequencing and Molecular Phylogenetic Analysis

DNA Extraction and Sequencing

Total DNA was extracted and processed with the Genomic DNA PowerSoil kit (Macherey-Nagel, Düren, Germany), following the manufacturer’s protocol. The cell pellet for each isolate was extracted in a 2 mL microcentrifuge tube with 400 μL SL1 lysis buffer (Macherey-Nagel, Düren, Germany) by homogenization in a FastPrep FP120 (Bio 101, Thermo Savant, Illkirch, France) at 6.5 m s^–1 for 45 s. The DNA purity and quantity were assessed by UV-spectroscopy with a NanoDrop ND-1000 system (Peqlab, Erlangen, Germany) and DNA integrity was confirmed by 1% agarose gel electrophoresis.

The D1/D2 region of the 28S large subunit (LSU) rDNA and the internal transcribed spacer (ITS) region (including the ITS1, 5.8S subunit, and ITS2 sequences; referred to hereafter as ITS/5.8S) were amplified from each total DNA extract by polymerase chain reaction (PCR). The forward and reverse primers for LSU amplification were: D1R-F (5′-ACC CGC TGA ATT TAA GCA TA-3′) and D2C-R (5′-CCT TGG TCC GTG TTT CAA GA-3′), respectively. The forward and reverse primers for ITS/5.8 amplification were: ITS a (5′-CCA AGC TTC TAG ATC GTA ACA AGG (ACT)TC CGT AGG T-3′) and ITS b (5′-CCT GCA GTC GAC A(GT)A TGC TTA A(AG)T TCA GC(AG) GG-3′), respectively. Reaction conditions were as follows: HotMaster Taq^® (5 Prime, Hamburg, Germany) buffer 1×, 0.1 mM of dNTPs, 0.1 mM of each forward and reverse primer and 1.25 units of Taq polymerase were added to 10–30 ng of the extracted genomic DNA in total reaction volumes of 50 μL. For the 28S rDNA amplifications, the reactions were subjected to the following thermocycling conditions: one cycle at 95°C for 7 min; 35 cycles at 94°C for 45 s, hold at 54°C for 2 min, 70°C for 1.5 min; and a final extension at 70°C for 5 min. The thermocycling conditions for the ITS/5.8 amplifications were: one cycle at 94°C for 4 min; 9 cycles at 94°C for 50 s, hold at 60°C for 40 s and then at 70°C for 1 min; 29 cycles 94°C for 45 s, then at 50°C for 45 s, 70°C for 1 min; and a final 5 min extension step at 70°C.

Samples were kept at 4°C until analysis on 1% agarose gel to ensure the expected amplification products were present. After PCR amplification LSU amplicons were sequenced and ITS/5.8S amplicons were subsequently cloned into the vector provided with the TOPO TA Cloning^® kit (Invitrogen, Carlsbad, CA, United States); three to five clones per amplicon were sequenced using the M13 vector primers supplied with the kit. Sequencing was conducted with standard cycle sequencing chemistry ABI 3.1 (Applied Biosystems, Darmstadt, Germany). Sequencing products were analyzed on an ABI 3130 XL capillary sequencer (Applied Biosystems, Darmstadt, Germany) and the generated sequences were analyzed and assembled with CLC main workbench version 12.0.¹ The resulting sequences (Supplementary Table 1) were submitted to GenBank (Accession numbers: MZ308606-MZ308617 for the ITS, and MZ310155-MZ310171 for the LSU rDNA sequences).

Molecular Taxonomic and Phylogenetic Analysis

Taxonomic sampling for phylogenetic analysis included representatives of P. lima sensu lato from most of its known global distribution. Sequences from strains labeled as P. lima, P. hoffmannianum, and P. belizeanum from the Americas, Asia, and Europe were retrieved from GenBank². Sampling was restricted to LSU and ITS rDNA and most of the sequences were gathered from GenBank to be able to assemble a data set that included each locus from the same species. The respective rDNA loci of cultured isolates collected in Mexico, provisionally identified by light microscopy as referable to P. lima sensu lato and P. hoffmannianum (Table 1), were amplified, sequenced, and included in the alignment, yelding a full data set for phylogenetic analysis. Two data sets with 13 species were assembled, one for the LSU locus with 68 sequences, and another for the ITS locus with 72 sequences (Supplementary Table 1). Taxonomic sampling mainly followed Chomérat et al. (2019), with the addition of recently described species in the genus (Luo et al., 2017; Nascimento et al., 2017; Lim et al., 2019). Each data set was analyzed with maximum likelihood and Bayesian inference approach under the same conditions as described below.

For each locus, multiple sequence alignment was performed independently in MAFFT 7.471 (Katoh and Standley, 2016) using the option auto, and visually inspected and corrected in Mesquite 3.61 (Maddison and Maddison, 2019), following suggestions by Morrison (2015). A Bayesian inference analysis was conducted in MrBayes 3.2.7 (Ronquist et al., 2012) through CIPRES web portal (Miller et al., 2010). Each locus was analyzed independently under the following conditions: nucleotide substitutions were modeled with a reversible jump approach (Huelsenbeck et al., 2004), allowing a mixture of models weighted by its posterior probability (PP) (Xie et al., 2011), rates heterogeneity was modeled with a gamma distribution (Yang, 1994; Huelsenbeck and Rannala, 2004), with four categories. The Monte Carlo Markov Chain was set to 10 million generations, with burning of 25%. For the specific case of ITS region, temperature was set to allow appropriate mixing. For both analysis, mixing and convergence were monitored in the output from MrBayes and in Tracer 1.7.1. (Rambaut et al., 2018).

Character support for clades was estimated with a non-parametric bootstrap (BS) (Felsenstein, 1985) in IQ-TREE 1.6.12 (Nguyen et al., 2014) using a ultrafast BS approach (Minh et al., 2013) with 1000 replicates. A Maximum Likelihood (ML) analysis was performed but only for support values. The ML analysis not only recovered the same branches, but the tree is the same for the relevant nodes. The PP and BS values are quite similar; nodes with high PP have high BS support, yielding virtually the same results.

Secondary Structure Modeling of Internal Transcribed Spacer Region 2 of Ribosomal RNA

The compensatory base change (CBC) approach was applied to investigate putative speciation among dinoflagellates based on potential hybridization barriers among isolates from the different clade deduced from the ITS/5.8S phylogeny (Gottschling and Plötner, 2004; John et al., 2014). The pattern of CBCs was determined from sequences of the representative isolates from Mexico from this study for each phylogenetic ITS/5.8S clade. These sequences were aligned with those of other taxa belonging to the PLSC with strains identified via BLAST searches in GenBank. The beginning and end of the ITS2 region in these strains were identified via homology with sequences of P. lima (GenBank accession FJ823582) retrieved from the ITS2 database (Koetschan et al., 2010), and were confirmed by comparison with a broader alignment of dinoflagellates. A multiple sequences alignment was performed from isolates collected from three regions in Mexico, as well as for representatives of clades III and IV. Sequences were aligned manually in Mesquite 3.61 (Maddison and Maddison, 2019). The ITS2 region was identified by molecular homology, and excised in a different file; the ITS2 segment lengths range from 156 to 186 bp.

The ITS2 secondary structure was predicted from the existing structural “a” model for the Prorocentrum ITS2 region (Koetschan et al., 2010), in agreement with the general dinoflagellate structures predicted by Gottschling and Plötner (2004). An identity matrix with a threshold of 20% was selected, with a sequence of P. lima retrieved from the same database as reference. The predicted structures were exported in an extended FASTA format; the CBCs in the predicted secondary structure were inspected in 4SALE (Seibel et al., 2006), and visualized using the Bruccoleri algorithm (Bruccoleri and Heinrich, 1988).

Diarrheic Shellfish Toxin Extraction and Analysis

Toxin Extraction From Prorocentrum Cells

Cell pellets were extracted in 500 μL 50% methanol in FastPrep tubes, and after adding 0.9 g FastPrep lysing matrix D, were homogenized by reciprocal shaking in a FastPrep FP120 (Bio 101, Thermo Savant, Illkirch, France) at 6.5 m s^–1 for 45 s. Homogenized samples were centrifuged at 16,000 × g for 15 min at 11°C. The supernatant from each sample was transferred into a 0.45 μm pore-size spin-filter (Millipore Ultrafree, Eschborn, Germany) and centrifuged at 11,000 × g for 1.5 min at 11°C. Finally, filtrates were transferred into 2 mL LC-autosampler vials (Agilent, Waldbronn, Germany) for LC-MS/MS analysis.

Toxin Analysis by Liquid Chromatography Coupled With Tandem Mass Spectrometry

The analytical system consisted of an SCIEX-4000 Q-Trap (Sciex, Darmstadt, Germany), triple quadrupole mass spectrometer equipped with a TurboSpray^® interface coupled to an 1100 liquid chromatograph (LC) (Agilent, Waldbronn, Germany). The LC equipment included a solvent reservoir, in-line degasser (G1379A), binary pump (G1311A), refrigerated autosampler (G1329A/G1330B), and a temperature-controlled column oven (G1316A).

Separation of toxins was achieved following previous protocols (Krock et al., 2008; Nielsen et al., 2013) after injection of 5 μL extract onto an analytical column packed with 3 mm Hyperclone^TM 3 μm BDS C8 130 Å 50 × 2 mm (Phenomenex, Aschaffenburg, Germany), maintained at 20°C. The flow rate was 0.2 mL min^–1 and gradient elution was performed by eluents A (water, formic acid, and ammonium formate) and B (acetonitrile, formic acid, and ammonium formate). Initial conditions were 12 min column equilibration with 5% B, followed by a linear gradient to 100% B in 10 min and isocratic elution until 16 min with 100% B. The program was then returned to initial conditions until 19 min (total run time: 31 min).

Detection of DSP toxin analogs was performed by LC-MS/MS by multiple reaction monitoring (MRM) experiments carried out in positive-ion mode with selected mass transitions (Supplementary Table 2). The following parameters were applied: curtain gas: 10 psi, CAD gas: medium, ion spray voltage: 5500 V, temperature: no heating, interface heater: on, declustering potential: 50 V, entrance potential: 10 V, exit potential: 15 V and dwell times of 100–200 ms per transition. Due to detection of a putative novel DTX1 isomer and associated diol-ester, provisionally dubbed DTX1a and DTX1a-D8, respectively, a collision-induced dissociation (CID) spectrum of DTX1 and the novel DTX1a compound was recorded in an enhanced product ion (EPI) mode of m/z 836 (mass range m/z 150–800) in the positive mode. Mass spectrometric parameters were the same as in the respective SRM experiments (Nielsen et al., 2013).

Analytical standards of okadaic acid (OA) (1 ng μL^–1), dinophysistoxin 1 (DTX1) (500 pg μL^–1), and dinophysistoxin 2 (DTX2) (500 pg μL^–1) from the Institute for Marine Biosciences, National Research Council, Halifax, NS, Canada were used to identify and quantify DSP toxins in extracts of Prorocentrum cells. Due to lack of standards for derivatives such as OA diol-ester (OA-D8) and dinophysistoxin 1-diol ester (DTX1-D8), cell quotas were expressed as OA and DTX1 equivalents, respectively, considering the following detection limits for OA (47 pg μL^–1), DTX1 (35 pg μL^–1), and DTX2 (25 pg μL^–1). Data acquisition and processing was performed with the Analyst Software Version 1.5, Sciex.

Toxin Data Analysis

Statistical analysis of the toxin data was conducted with the software Infostat (V2019). Normality of toxin composition (OA, OA-D8, DTX1, DTX1-D8, DTX1a, and DTX1a-D8), cell toxin content, geographical localities (BLP, AAR, IV, and PLV) and different substrate types (Tracheophyta, Ochrophyta, Chlorophyta, Rhodophyta, and detritus) were examined by the Shapiro–Wilk test. Non-parametric analysis was performed on all the data sets because they did not show a normal distribution. The Wilcoxon-Mann-Whitney test was performed to determine differences in toxin composition among P. lima and P. hoffmannianum defined groups. Significant differences in cell toxin content among groups (e.g., geographical localities and substrate types) were assessed with the Kruskal–Wallis test.

Ordination by principal components statistical analysis was performed using R version 4.1.0 (R Core Team, 2020; RStudio-Team, 2020). Prior to ordination analyses, relative toxin compositions were Hellinger-transformed using the deconstand() function of the package vegan (V2.5.7). Principal components of the dataset were calculated using the prcomp() function and plotted using the ggbiplot package (V0.55). Betadispersion (homogeneity of variance) of individual groups (i.e., based on locality of origin, clade and substrate type) was tested with the betadisper() function followed by permutational ANOVA (PERMANOVA) (Anderson, 2001) using the adonis2() function in vegan. Boxplots of differences in toxin compositions were visualized with the ggplot2 package (V3.3.5).

Results

Several global biogeographical groups were defined based upon the location of origin of the isolates used for the morphological comparison and phylogenetic analysis of LSU and ITS rDNA sequence data for the PLSC and PHSC. These provisional groupings were established with reference to distributional regions for toxigenic benthic dinoflagellates from Latin America (Durán-Riveroll et al., 2019): Eastern Asian Pacific (China), Northeast Atlantic (Iberia, France), Mediterranean (Greece), Neotropical Southwest Atlantic (Brazil), Neotropical West Atlantic (Gulf of Mexico, Caribbean, Florida), Neotropical Pacific (Gulf of California, Mexico), and Tropical Australia. Not all biogeographical groups formed distinct clusters in the current analysis herein.

Morphospecies Description and Identification

The illustrated Prorocentrum isolates from the coasts of Mexico (Figures 2, 3) can be categorized into one of three morphotypes: members of the PLSC, PHSC, or P. rhathymum. Cells of all PLSC isolates (Figures 2A–E, 3A–E) have a smooth thecal surface and the central area free of trichocyst pores, with slight differences in cell outline in lateral view and size as well as in shape and size of pores (round, oval, or bean-shaped). Slight differences in cell outline and size among PLSC isolates show infraspecific variation, rather than distinct characteristics of an individual strain. The PLSC isolates vary mainly in cell outline in thecal view. For example, PA78 cells are more rounded compared to other isolates, with a lower L/W ratio (1.21 ± 0.04) (Table 1). Cells of PA78 cells are the longest (46.36 ± 1.33 μm), but also show maximum width at or slightly behind the median line (38.36 ± 1.36 μm), whereas other illustrated PLSC isolates exhibit maximal width at 1/3 to 2/5 of the cell length toward the posterior end. Isolates PA48 (Figures 2D, 3D) and PA72 (Figures 2E, 3E) show a cell outline similar to PA8.

FIGURE 2

Figure 2. Light micrographs of selected isolates of Prorocentrum species as representative for major clades: (A–E) P. lima morphotypes (PLSC) (PA1, PA8, PA9, PA48, and PA72, respectively) in right thecal view; (F) P. rhathymum (PA20) in left thecal view; (G,H) P. hoffmannianum (PA89 and PA70) in right thecal view. In all isolates a central ring-shaped pyrenoid can be distinguished. The posteriorly located nucleus with visible chromosomes can be seen. In most illustrated cells from the PLSC (A,B,E,F) and P. hoffmannianum (G,H) the maximum cell width is situated below the median line, closer to the posterior end, whereas for PLSC isolate PA9 (C) and P. rhathymum (F) it is at the median line.

FIGURE 3

Figure 3. Scanning electron micrographs of selected isolates of Prorocentrum species as representative for major clades: (A–E) P. lima morphotypes (PLSC) (PA1, PA8, PA9, PA48, and PA72, respectively) in right view; (F) P. rhathymum (PA20) in left thecal view; (G–I) P. hoffmannianum (PA89 and PA70); panels (G,H) in right thecal view; panel (I) in right lateral-apical view. Panel (I) shows a typical platelet arrangement and the two pores in the periflagellar area; numbering/lettering system for the platelets and pores in the periflagellar area: (e.g., ap–accessory pore, fp–flagellar pore) is according to Hoppenrath et al. (2013). Scale bars: 10 μm in Panels (A–H) and 5 μm in panel (I).

Cells of P. hoffmannianum also exhibit considerable morphological variation within the same clonal isolate. Variation in pore pattern similar to that among PLSC isolates can be also observed for P. hoffmannianum cells (Figures 2G,H), with more densely scattered pores across the thecal surface and the central area free of pores; cells may be without (Figure 2G) or with (Figure 2H) a marginal row of pores. Numerous subcircular depressions are visible around the theca, creating a rough thecal surface. Some PLSC isolates and P. hoffmannianum cells reveal tiny depressions along the valve margins. Isolates PA70 and PA89 that belong to P. hoffmannianum show the cell outline (Figures 2G,H) and L/W ratio (1.21–1.33) (Table 1), characteristic of this species (Faust, 1990b). In Figures 3G–I, P. hoffmannianum isolates are depicted in right thecal and in right lateral-apical view, with Figure 3I showing a typical platelet arrangement and the two pores in the periflagellar area.

Prorocentrum rhathymum cells (PA20) have sparsely scattered double-rimmed trichocyst pores on a smooth thecal surface, tending to form radial rows at both sides of the posterior end and lacking in the subcentral area of the major thecal plates (Figures 2F, 3F). PA20 cells coincide with the original description of P. rhathymum (Loeblich et al., 1979), although in the illustrated cells in lateral view, the two major plates narrow slightly more to the posterior end than anteriorly. In many cells of the isolate this feature is absent, reflecting infraspecific variation.

Molecular Genetics and Phylogeny

Separate multiple sequence alignments were produced for the ITS and LSU loci because there was little overlap between the sequences produced for both loci. We were only successful in amplifying and sequencing both loci for two strains from among the isolates from Mexico. We decided therefore not to perform a concatenated analysis with the two loci, but rather used the ITS alignment for further analysis. The ITS region is considered more informative in Prorocentrum at the infraspecific level and because it provides the opportunity of the opportunity for secondary structure analysis. In the present study, the higher intraspecific phylogenetic resolution and greater availability of reliable sequences focuses further discussion primarily on the ITS/5.8S region (Figure 4). Analysis for the LSU rDNA locus is also presented, but as Supplementary Material (shown in Supplementary Table 1 and Supplementary Figure 1), and only for broad context and comparison with published literature.

FIGURE 4

Figure 4. Majority rule consensus of the trees sampled in the Bayesian inference phylogenetic analysis of the ITS/5.8S rDNA sequences of Prorocentrum species, represented by isolates from the Gulf of California, Gulf of Mexico, and Caribbean coasts of Mexico incorporated into a global clade framework from geographical regions. The ITS2 variants of secondary structure (labeled as 1, 2, and 3) correspond with secondary structure of the specimens sampled in this study on each clade. Variant 1 corresponds with clade I, variant 2 with clade II, and variant 3 with P. hoffmannianum. Arrows and deep red dots indicate the distinctive compensatory base changes (CBCs) for each secondary structure for each lineage in P. lima. Values above branches correspond to Bayesian posterior probabilities/bootstrap frequencies, for main branches.

A multiple sequences alignment, 670 bp in length, was produced for the ITS locus, including 72 sequences for the ITS/5.8 rDNA alignment, and 68 sequences for the LSU rDNA alignment (Supplementary Table 1). After 10 million generations of the Bayesian inference analysis, the average standard deviation of split frequencies was 0.003; chain swap values were satisfactory, with effective sample size 231 or above for all parameters. Convergence for both chains as well as mixing was considered satisfactory after inspection. Posterior probabilities (PP) and bootstrap (BT) values were drawn on the majority rule consensus produced with MrBayes.

Based upon ITS/5.8S rDNA sequences, representatives of P. lima sensu lato were recovered in five well supported clades (Figure 4), and the overall relationships between the five clades were well supported by PP and BT. Most isolates from Mexico were recovered in Clades I and II. Six isolates from this study are in Clade I, while three from PLV, Veracruz fall into Clade II with isolates from Brazil. In any case, defined on the basis of ITS/5.8 sequences, Clade I does not exhibit global biogeographical consistency, as it includes strain representatives from China, Florida, and three distinct localities from Mexico.

On the ITS/5.8S rDNA tree, strains from Mexico in Clade II are all from PLV, Veracruz, but this clade also includes those from geographically distant Brazil. Clades III and IV, both well supported, do not include any Prorocentrum isolates collected in Mexico, but exhibit higher congruence with geographical origin. Clade III comprises strains from China, whereas Clade IV includes those from the Northeast Atlantic (Iberia) and Mediterranean Sea. The remaining three isolates of the P. lima sensu lato complex from PLV, Veracruz (PA86, PA89, and PA91) occur together with P. hoffmannianum/P. belizeanum morphotypes from Florida and Ibiza (western Mediterranean Sea), respectively, within Clade V, with two P. lima strains from China in its basal branches.

A multiple sequences alignment of 707 bp for 45 strains was produced from the LSU rDNA data set. After 10 million generations of the Bayesian inference analysis, the average standard deviation of the split frequencies was 0.004; effective sample size for all parameters was 234 or above. Convergence for both chains as well as mixing was considered satisfactory after visual inspection in Tracer 1.7.1.

A comparison of the respective phylogenetic trees for ITS/5.8S rDNA (Figure 4) with LSU rDNA (Supplementary Figure 1) shows minor conflicts, but both trees are basically congruent. The topology of the phylogenetic tree based only on LSU rDNA sequences (Supplementary Figure 1) also recovers the same five clades as from ITS analysis, but relationships between clades are slightly different and poorly supported. The clade structure of the phylogenetic tree based only on LSU rDNA sequences also groups strains from Brazil together with those from the Gulf of Mexico (primarily PLV, Veracruz), but those from China are widely distributed among clades. The association of P. lima PL6 from Isla El Pardito, Baja California Sur, Mexico forms a special case. On the ITS/5.8 tree (Figure 4) it belongs in Clade 1 with the other PLSC members from Mexico but is situated with strains from China and a few others from Mexico in another clade on the LSU tree (Supplementary Figure 1).

The ITS sequences of strains from Mexico and the respective phylogenetic ITS2 clades exhibit a distinctive pattern of CBCs for their given ITS2 secondary structure (Figure 4 and Supplementary Figure 2). The ITS2 region consists of 156–186 bp in the Prorocentrum isolates from the coastal locations in Mexico (Supplementary Table 1). Within Clade I and Clade II, the strains from Mexico also fail to exhibit a clear geographical or morphological segregation, but compared to the basal Clade V, they do share a distinctive pattern of CBCs in their ITS secondary structure. Sequences from Clades I and II have distinctive structural patterns: Clade I members have a distinctive CBC at helix 3 (alignment positions C142 and G107), and sequences from Clade II have two CBCs at helix 2 (alignment positions A58 and T64, and C54 and G67) and one CBC at helix 3 (alignment position T95 and G155) (Figure 4 and Supplementary Figure 2). This distinctive pattern may be indicative of early stages of speciation.

Toxin Composition and Cell Toxin Content

High variation in DSP toxin composition was exhibited among toxigenic clonal isolates (n = 67 considering that 2 isolates were non toxigenic) from three disjunct geographical locations on the Gulf of California, Gulf of Mexico, and Caribbean coasts of Mexico. The only isolate of P. rhathymum (PA20) contained no detectable DSP toxins, but seven of the eight P. hoffmannianum isolates (all except PA71) were also toxigenic.

All isolates provisionally assigned to the PLSC (Table 1) produced substantial quantities of DSP toxins, readily detectible in methanolic extracts (ca. 6.7 × 10⁵ cells mL^–1) by LC-MS/MS (Figure 5). Confirmed toxin spectra by LC-MS/MS revealed the presence of six major analogs, including four well-known derivatives (OA, DTX) and their respective diol-esters (OA-D8, DTX1-D8), plus two novel DTX1 analogs (DTX1a, DTX1a-D8) of unconfirmed structure. Among all Prorocentrum isolates, the following analogs were detectible in decreasing order of mean relative abundance: OA > OA-D8 > DTX1 > DTX1a > DTX1-D8 > DTX1a-D8. Toxin composition showed analogs OA, OA-D8, and DTX1 as the dominant toxins among PLSC and PHSC isolates. The relative toxin composition of OA, OA-D8, DTX1, DTX1a, and DTX1a-D8 were significantly higher for P. lima than for P. hoffmannianum isolates, but relative content of DTX1-D8 could not be distinguished between these species (total n = 67, p < 0.05, Wilcoxon-Mann-Whitney test).

FIGURE 5

Figure 5. Relative composition (% total) of DSP toxins (OA and DTX analogs) in isolates (n = 67) of the P. lima complex (PLSC) and P. hoffmannianum (indicated by ^∗asterisk) from the Gulf of California, Gulf of Mexico, and Caribbean coasts of Mexico. BLP, BCS = Bahía de La Paz, Baja California Sur; AAR, VRS = Anegada de Adentro Reef, Veracruz Reef System, Veracruz; IV, VRS = Isla Verde, Veracruz Reef System, Veracruz; PLV, VRS = Piedra de La Virgen, Veracruz Reef System, Veracruz; PM = Puerto Morelos, Quintana Roo. Toxin acronym codes: OA = okadaic acid; OA-D8 = okadaic acid diol-ester; DTX1: dinophysistoxin 1; DTX1-D8 = dinophysistoxin 1 diol-ester; DTX1a = undescribed isomer of DTX1; DTX1a-D8: undescribed isomer of DTX1a-D8.

The two novel DTX analogs (DTX1a, DTX1a-D8) were stably produced among toxigenic isolates and readily detectable by LC-MS/MS analysis of toxin composition. Even without confirmed structures, they were therefore considered valid chemotaxonomic markers in the principal components analysis (PCA) (Figure 6) and PERMANOVA box-and-whisker plots of statistical significance (Figure 7). According to the PCA (Figure 6), the relative cell abundance (% molar) of DTX1 and minor novel variants DTX1a and DTX1a-D8 are strongly associated with isolates from Piedra de la Virgen, Veracruz Reef System (PLV, VRS), whereas high OA is characteristic of Isla Verde (IV, VRS) and Bahía de La Paz, Baja California Sur (BLP, BCS). The diol-ester analogs OA-D8 and DTX1-D8 are more prominent in isolates from Anegada de Adentro Reef (AAR, VRS) and BLP, BCS. Permutational multivariate analysis of variance (PERMANOVA in R, 999 permutations) comparison of relative toxin composition also shows differences between groups defined by locality of origin, ITS-defined clade, and substrate type, based on a pseudo F-statistic (^∗∗∗ ≤ 0.0001; ^∗∗ = 0.001; ^∗ = 0.01) (Figure 7).

FIGURE 6

Figure 6. Principal component analysis and ordination of Prorocentrum isolates (n = 67) along the first and second PC axes representing 70.1% of the observed variance with relative toxin composition as explanatory variables. Samples are color-coded to reflect grouping by locality of origin.

FIGURE 7

Figure 7. Boxplot of Prorocentrum isolates according to relative toxin composition (% total) (n = 67) in different groups. (A) Locality of origin; (B) Clade, (C) Substrate type. Group differences in toxin composition were tested with PERMANOVA (999 Permutations); asterisk indicates significant difference between groups based on pseudo-F-statistic (** < 0.01; * < 0.05) between individual groups.

Only isolates for which verified cell counts were available were selected for calculation of cell toxin content (fmol cell^–1) (Table 2 and Supplementary Figures 3A,B). For statistical analysis, isolates collected from detritus as dead Ochrophyta (Phaeophyceae) (PA46 and PA51) were considered to belong to both detritus and Ochrophyta groups, but the PM location was excluded from the geographical comparison due to lack of replicate isolates (n = 1). Among isolates morphologically belonging to the PLSC, the mean total cell toxin content among geographical localities (BLP n = 4, AAR n = 2, IV n = 10 and PLV n = 2) showed no significant differences (n = 18, p < 0.05, Kruskal–Wallis test). Among epiphytic isolates collected from different substrate types, Tracheophyta (n = 7), Ochrophyta (n = 6), Chlorophyta (n = 3), Rhodophyta (n = 3), and detritus (n = 2), there were no significant differences in toxin content (n = 21, p < 0.05, Kruskal–Wallis test).

TABLE 2

Table 2. Cell content of DSP toxins (OA and DTX plus analogs) from cultured benthic Prorocentrum isolates assigned to the P. lima complex (PLSC).

Discussion

Membership in the PLSC has been defined originally on the basis of consistent morphological features, such as the smooth thecal surface, the structure of the periflagellar area and the presence of a central pyrenoid (Aligizaki, 2009). Within the genus Prorocentrum there remain several similar but unresolved taxa provisionally assigned to the PLSC as proposed by Aligizaki (2009) and subsequently revised with more molecular genetic support based on rDNA sequencing (Hoppenrath et al., 2014; Chomérat et al., 2019). Classification at the species level, however, is not always consistent between morphological vs. molecular genetic criteria. Early molecular studies of 18S rDNA gene sequences (Grzebyk et al., 1998) indicated the strong similarity between P. lima and P. arenarium, with only 11 nucleotide substitutions between them and several shared synapomorphies. These congruencies led to classification of P. arenarium as a round P. lima morphotype (also later supported by Zhang et al. (2015). P. lima, including P. arenarium, was proposed as monophyletic (Nagahama et al., 2011), but divided into several subclades correlated broadly to sample collection locations. No specific morphological characters were found to be related to any one subclade. Based on LSU and ITS/5.8S rDNA sequence phylogeny (Zhang et al., 2015) found very large infraspecific genetic differences in P. lima among five morphotypes from Hainan Island, South China Sea. All morphotypes fell within the morphological description of P. lima and P. arenarium, thereby leading to the conclusion that the morphological and molecular genetic data are contradictory.

Among and within defined benthic Prorocentrum species the high degree of phenoplasticity has certainly led to misidentifications. The presence of aberrant cells and multiple morphotypes varying in cell size, linear dimensions, number, and distributional patterns of thecal pores are frequently encountered in benthic Prorocentrum cultures, particularly in long-term ones (Y. Okolodkov, L. Durán-Riveroll, and A. Cembella, pers. observ.). The current analysis of cell size and shape variation among clonal cultures indicates that these are not generally reliable characteristics to distinguish among closely related species (e.g., P. lima vs. P. hoffmannianum) nor among infraspecific clones isolated from different substrate types or geographical locations within Mexican coastal waters. Faust (1990b), who originally described P. hoffmannianum from naturally collected specimens notes that cells of this species are typically larger and broader than those of P. lima. The cell size range (45–55 μm long and 40–45 μm wide) of P. hoffmannianum from Twin Cays, Belize is consistent with the dimensions of this species in cultures from Mexico, but the range is too broad to be diagnostic.

That such morphological variation can occur even in clonal cultures under standardized growth conditions in an apparently homogeneous environment indicates low constraints on expression of such key morphological features for taxonomic description based on autapomorphic criteria. If extrapolated to the natural benthic environment this may have led over-classification of morphotypes in some cases. On the other hand, the overall shape and structure characteristic of all benthic Prorocentrum, flattened heavily walled cells often adhering by mucus together as colonies upon substrates, are particularly adaptive for maximizing adhesion as epiphytes even under turbulent conditions and may provide barriers to grazing and space-competitors.

Such apparent contradictions between morphology and molecular genetic characteristics highlight the critical role of assessing cryptic speciation in defining diversity. There is indeed a high likelihood of cryptic speciation within “P. lima complex,” as originally suggested by Aligizaki (2009). In the current study of isolates from the three sub-regions of Mexico, all morphotypes consistent with P. lima are provisionally referred to the PLSC. In general, this follows the proposal by Aligizaki (2009) and recommendation by Hoppenrath et al. (2014), but with a minor amendment. The PLSC should be recognized as a complex of individual species and not of infraspecific taxa within P. lima sensu lato. Hoppenrath et al. (2014) consider Prorocentrum consutum Chomérat et Nézan among species similar to P. lima. The lack of the P. arenarium morphotype from the sampling in Mexico does not allow us to resolve its distinction from P. lima. We therefore consider them synonymous, following Nagahama et al. (2011), Hoppenrath et al. (2014), and Lassus et al. (2016). A recently described species from Brazil, Prorocentrum caipirignum S. Fraga, M. Menezes et S. Nascimento, is similar to P. lima in its morphology, and appears close to the clades of P. lima and P. hoffmannianum in previous phylogenetic trees based on ITS and LSU rDNA sequences (Nascimento et al., 2017). No morphotypes consistent with P. caipirignum were found among populations studied from Mexico, but this species was included in the ITS and LSU-based phylogenetic alignments herein, and the close relationship with the PLSC was confirmed.

Noting the uncertainty in the species boundaries between P. hoffmannianum and P. belizeanum, Hoppenrath et al. (2013) would also include Prorocentrum sabulosum M. A. Faust and Prorocentrum tropicale M. A. Faust in the PHSC. According to Herrera-Sepúlveda et al. (2015), P. lima sensu lato is a sister to the P. hoffmannianum/P. belizeanum/P. maculosum group. This phylogenetic analysis concludes that P. belizeanum M. A. Faust is conspecific with P. hoffmannianum, forming a highly supported monophyletic clade within the PHSC, but with three subclades broadly corresponding to the sample collection regions. Our current phylogenetic analysis support the paraphyly of the PLSC, wherein the PHSC with its autapomorphic feature forms a monophyletic subclade of Clade 5 (Figure 4 and Supplementary Figure 1; and see Chomérat et al. (2019).

Based on SSU rDNA analysis, P. hoffmannianum, P. belizeanum, and P. maculosum comprise a highly supported monophyletic group (Herrera-Sepúlveda et al., 2015). P. maculosum M. A. Faust is also morphologically very similar to members of the PLSC and P. hoffmannianum (Hoppenrath et al., 2013), and hence was concurrently misidentified as P. concavum (Glibert et al., 2012). This species is considered to be a sister group to the P. hoffmannianum/P. belizeanum clade.

Prorocentrum rhathymum is morphologically distinct from both the PLSC and PHSC, but still remains an enigma. Often considered a synonym of P. mexicanum B. F. Osorio (Steidinger and Tangen, 1996; Faust and Gulledge, 2002), P. rhathymum and P. mexicanum have also been recognized as distinct species (Cortés-Altamirano and Sierra-Beltrán, 2003), based on both morphological and ecological (habitat) differences. Phylogenetic analysis of strains from Florida (An et al., 2010) supports this separation. On the contrary, molecular characterization of Korean strains and their comparison with other P. rhathymum and P. mexicanum strains from different geographical regions did not allow Lim et al. (2013) to classify them as separate species. Molecular phylogenetic analysis of P. rhathymum strains from Florida (An et al., 2010) vs. P. cf. rhathymum from Korea has served to highlight inconsistences in species level discrimination (Lim et al., 2013) and has led to embedding these groups in the paraphyletic PLSC clade (Chomérat et al., 2019).

A comparison of the respective phylogenetic trees for ITS/5.8S rDNA (Figure 4) with LSU rDNA (Supplementary Figure 1) shows minor conflicts but both are basically in congruence with the literature on species associations within Prorocentrum [e.g., Chomérat et al. (2019)]. Nevertheless, comparative analysis of the LSU sequences (Supplementary Figure 1) among new strains from geographically distant sub-tropical locations (Gulf of Mexico, Caribbean and Gulf of California) integrated with global strains of benthic Prorocentrum species does not clearly resolve the PLSC at the species level. This remains paraphyletic and does not show clear geographical separation among the clades. The topology of the ITS/5.8S rDNA phylogenetic tree (Figure 4) recovered is generally congruent with previous molecular phylogenetic analyses of the PLSC and association with P. hoffmannianum and allies, but with minor differences. A Maximum Likelihood (ML) analysis based on an LSU rDNA matrix of Prorocentrum species (Chomérat et al., 2019) shows the PHSC in a tight cluster within the PLSC node/clade. Again, in all phylogenetic analyses the PLSC appears to be paraphyletic and the P. hoffmannianum group (PHSC) branches after some basal P. lima isolates (Figure 4 and Supplementary Figure 1) and is consistent with the scheme published in Chomérat et al. (2019).

Although we are aware that there is still limited and perhaps not representative sampling coverage of genotypes from the different target regions, reference to the ITS/5.8 sequences rather than LSU successfully resolves some of the issues regarding the scope of the PLSC and cryptic speciation. Nevertheless, there are also intriguing anomalies to consider from the perspective of biogeographical origins. Benthic dinoflagellates, including epiphytic Prorocentrum are typically rather sessile and tend to remain associated with the local substrate (Durán-Riveroll et al., 2019). In principle, this suggests that endemic populations have reduced potential for genetic exchange with geographically disjunct populations of emerging species. Knowledge of the sexual life history stages of benthic Prorocentrum is extremely limited, and recombination events are rarely observed (unlike for the phylogenetically distinct planktonic species). Putative sexual (resting) cysts have been reported for P. lima (Faust, 1993; Faust et al., 1999) and another benthic species Prorocentrum marinum (Cienkowski) Steidinger et Williams ex Head (Faust, 1990a, 1993), but these reports are typically inferences from morphological analysis of field specimens rather than confirmed from mating and germination studies. In fact, true fossilizable sexual resting cysts have just recently been confirmed for the first time and in only one species, Prorocentrum leve M.A. Faust, Kibler, Tester et Litaker (Mertens et al., 2017). Circumstantial evidence therefore suggests that sexual recombination events may be rather rare within endemic benthic populations.

Potential hybridization barriers appear within the PLSC (Clade I–V) including the PHSC (Clade V). Indeed, the segregation of P. hoffmannianum in Clade V indicates clear speciation and is consistent with the morphological differentiation and divergence from P. lima. Patterns of micro-ecological niche differentiation cannot, however, be inferred from substrate type or local biogeographical barriers. Arguments for cryptic speciation could be interpreted from the existence of hybridization barriers by the presence of CBCs by pairwise comparisons in the ITS2 region between the clades (Figure 4), as supported by several independent studies (Fabry et al., 1999; Coleman, 2003, 2009; John et al., 2014). Though not mechanistically linked to the development of either pre- or post-zygotic incompatibility, when present, CBCs or hemi-CBCs consistently correlate with species-level divergences even among closely related lineages (Gottschling and Plötner, 2004; Gottschling et al., 2005; John et al., 2014). Not all speciation events, however, lead to the development of CBCs or hemi-CBCs. Consequently, the lack of CBCs between clades should not be interpreted as evidence against their being separate species. Direct evidence of speciation events based on cyst unviability are missing, but the accumulation of CBCs within Clades I–III indicates at least cryptic speciation is in progress.

An additional caveat in interpretation of “speciation” in Prorocentrum is related to morphotaxonomic and nomenclatural uncertainties, particular with the older literature on species descriptions. For example, in the current analysis, P. belizeanum is included in Clade V with P. hoffmannianum and some basal P. lima. As the name implies, P. belizeanum was described originally from Belize; its current distribution appears to be restricted to the Caribbean region (GBIF.Org, 2020). The inclusion of P. belizeanum Dn139EHU from Spain (as reported in GenBank) may therefore simply represent an original misidentification or more likely the fact that they are synonyms (Herrera-Sepúlveda et al., 2015), as also supported by the tight association in Clade V.

Biogeographical distributional patterns and substrate associations are key to understanding population dynamics and speciation among benthic Prorocentrum. No definitive association have been identified in temperate areas (Foden et al., 2005), whereas in tropical environments this genus has been preferentially associated with both filamentous turf (Parsons and Preskitt, 2007) and Phaeophyceae (Delgado et al., 2006) as substrates. The current analysis of morphotypes, rDNA-defined genotypes and toxin chemotypes failed to yield clear patterns dependent upon substrate types or biogeographical locations in Mexican coastal waters. Nevertheless, the clustering of Prorocentrum genotypes belonging to the PLSC from the eastern North Atlantic coast and the western Mediterranean Sea (Ibiza) within Clade II on the ITS/5.8S rDNA tree, and similar clustering on the LSU rDNA tree, do indicate a biographical barrier to genetic exchanges with Caribbean and Gulf of Mexico populations. On the other hand, benthic Prorocentrum species are notorious for and commonly found rafting upon detritus and detached macroalgae (Bomber et al., 1988; Faust et al., 1999) –a viable means of long distance dispersal by ocean currents. Based upon ITS/5.8S rDNA analysis, strains in Clade II are from one locality in Mexico (Veracruz, PLV), but this clade also includes a strain from geographically distant Brazil. Admittedly speculative, but mixing of genotypes may have occurred via rafting from the Brazilian coast by entrainment in the northwestward flow of the Atlantic South Equatorial Current linked to the Guiana Current and eventually transfer to the Gulf of Mexico in the Caribbean Current. Similarly, the close ITS genotypic association of P. hoffmannianum from the Gulf of Mexico and the Florida coast in Clade V may be plausibly explained by prevailing current transport of detritus and macroalgae bearing Prorocentrum epiphytes.

High chemodiversity in toxin composition was evident among Prorocentrum species and isolates differing in geographical origin, substrate type and clade assignment. In fact, clonal variation in toxin composition was substantially greater than as defined by morphological or molecular genetic criteria such as ITS/5.8S sequences. Since these toxin profiles are stable yet distinct under homogeneous culture conditions and throughout the culture cycle, this provides the opportunity to define genetically fixed chemotaxonomic markers. To our knowledge the current study represents the largest comparison (n = 67 isolates) of toxin compositional profiles among populations for any genus of marine benthic dinoflagellate. Despite this large data set, we were unable to confirm clear distinctions among substrates or geographical origins to define populations within a species. Between species P. lima and P. hoffmannianum there were apparent tendencies for toxin compositional differences for some analogs (OA, OA-D8, DTX1, DTX1a, and DTX1a-D8) but not for others (DTX1-D8). There were no statistical differences between these species in total toxin content. The toxin compositional data support the contention that production of DSTs is likely an obligate trait within the PLSC sensu stricto (Chomérat et al., 2019). Close relatives of P. lima, such as P. caipirignum (previously as P. cf. maculosum) (Yasumoto et al., 1984; Hu et al., 2017; Luo et al., 2017) from diverse populations from Malaysia (Lim et al., 2019), the South China Sea (Luo et al., 2017), and Brazil (Nascimento et al., 2017), have also been found as consistently toxigenic. Nevertheless, at this stage it is premature to conclude that this high chemodiversity can be effectively applied to distinguish genotypes or “species,” in the absence of knowledge on the biosynthetic gene clusters for these polyketides.

This multivariate study of diversity of benthic Prorocentrum yields valuable insights into the processes of incipient speciation and definition of appropriate taxonomic markers and criteria for distinguishing species. Sexual life histories are lacking for most benthic dinoflagellates, and particularly for Prorocentrum species. Molecular genetic markers are usually capable of discriminating “benthic,” “tychoplanktonic,” and “planktonic” taxa, but provide only a static view of behavior and species interactions. Only by combining comparative approaches to morphotaxonomy, breeding affinities, inheritance of toxigenicity and molecular diversity within and among global population will it be feasible to define the ecological niches and adaptive strategies to account for the survival and success of benthic dinoflagellate species. Such multidisciplinary approaches are necessary to distinguish the mechanisms among those of phylogenetically close relatives in planktonic assemblages.

Dedication

This publication is dedicated to the memory of our dear friend and colleague Maria A. Faust (21 April 1930–24 April 2021), recently deceased after a 40-year career at the National Museum of Natural History within the Smithsonian Institution (Washington, DC, United States). She was a true pioneer in biodiversity research, particularly on tropical marine dinoflagellates, discovering and describing dozens of new benthic species including many Prorocentrum discussed herein.

Data Availability Statement

The datasets presented in this study can be found in online repositories. Alignments availabe upon request. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author Contributions

LD-R and AC conceived the project, designed the field and laboratory studies in consultation with YO, extracted and prepared the rDNA for LSU and ITS sequencing, contributed equally, and led and coordinated the drafting of the manuscript. LD-R, AC, and YO conducted the field sampling for epibenthic dinoflagellates. LD-R prepared the preliminary isolates for culture and maintained the clonal strains at ICMyL-UNAM, Mexico, and LD-R, UT-J, and AC at AWI, Germany. UJ supervised the rDNA gene sequencing and together with RG-S conducted the phylogenetic and molecular taxonomic analysis. LD-R and YO conducted the morphotaxonomic analysis of Prorocentrum isolates. UT-J performed the toxin extractions for LC-MS/MS of DSP toxins and data analysis under supervision from BK. CH with AC and UJ performed the statistical analyses of environmental, toxin, and phylogenetic data. All authors contributed actively to the preparation of data content and writing of this manuscript.

Funding

This research was funded primarily by the Consejo Nacional de Ciencia y Tecnología (CONACyT, Mexico) Basic Science Project, number A1-S-8616 awarded to LD-R. Financial contribution for the preparation and publication of this manuscript and also associated research activities was provided to AC, BK, and UJ within the Helmholtz research program “Changing Earth, Sustaining our Future” (Sub-topic 6.2 Adaptation of marine life) of the Alfred-Wegener-Institut, Helmholtz-Zentrum für Polar- und Meeresforschung, Germany. The Deutscher Akademischer Austauschdienst (DAAD) made possible the research stay of LD-R at AWI, with focus on the molecular biological aspects. UT-J acknowledges the funding support from the International Max Planck Research School of Marine Microbiology for her master studies and research opportunities.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Manuel Victoria-Muguira of Dorado Buceo for boat access and logistical support with field sampling from Veracruz, and Manuel Rodríguez-Gómez from the Acuario de Veracruz for donating filtered seawater. Mario Fernando Sánchez Bernal was responsible for imaging and curation. Laura Elena Gómez-Lizárraga from the Scanning Electron Microscopy Laboratory at Instituto de Ciencias del Mar y Limnología, Universidad Nacional Autónoma de México (ICMyL-UNAM) provided micrographs; Manuel E. Martínez-Cruz, Hugo Pérez-López, and Apatzingán Palomino-Hermosillo assisted with SEM sample preparation and cell measurements; Andrea García-Casillas (UNAM) performed statistical data analysis. Technical assistance of Annegret Müller and Thomas Max for LC-MS/MS analysis of DSP toxins and Lena Nieke for Prorocentrum culture maintenance at AWI is much appreciated.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.716669/full#supplementary-material

Footnotes

References

Aguilar-Trujillo, A. C., Okolodkov, Y. B., Herrera-Silveira, J. A., Merino-Virgilio, F. D. C., and Galicia-García, C. (2017). Taxocoenosis of epibenthic dinoflagellates in the coastal waters of the northern Yucatan Peninsula before and after the harmful algal bloom event in 2011–2012. Mar. Pollut. Bull. 119, 396–406. doi: 10.1016/j.marpolbul.2017.02.074

PubMed Abstract | CrossRef Full Text | Google Scholar

Aligizaki, K. (2009). “Spread of potentially toxic benthic dinoflagellates in the Mediterranean Sea: a response to climate change?,” in Phytoplankton Response to Mediterranean Environmental Change, ed. C. Workshop (Tunisia: CIESM), 57–61.