Introduction

Sessile anthozoans form the basis of some of the most biodiverse and economically important ecosystems on Earth (Moberg and Folke 1999; Barbier et al. 2011). These communities, however, are increasingly threatened by anthropogenic disturbances (e.g., Hoegh-Guldberg et al. 2017). Of particular interest are Caribbean octocorals, commonly referred to as gorgonians, as their abundance has remained relatively stable or has increased in recent decades (Ruzicka et al. 2013; Lenz et al. 2015; Tsounis and Edmunds 2017; Lasker et al. 2020), while their scleractinian counterparts (i.e., stony corals) have suffered declines (e.g., Gardner et al. 2003; Côté et al. 2005; Roff and Mumby 2012). Our ability to conserve these ecosystems relies on assessing traits such as reproductive potential, dispersal, population connectivity, and abundance (Mumby and Steneck 2008), tasks that require an understanding of species boundaries and relationships among groups (Knowlton and Jackson 1994).

Despite their importance, there are clear gaps in our understanding of Caribbean octocorals (Kupfner Johnson and Hallock 2020; Lasker et al. 2020). The most comprehensive taxonomic treatment of the group was written six decades ago and relies on a mix of qualitative and quantitative distinctions in colony form and the morphology of sclerites (i.e., calcium carbonate spicules; Bayer 1961). While these characters have been useful in distinguishing some octocoral species (e.g., Aharonovich and Benayahu 2012), subtle differences in morphological traits often limits their utility for species delimitation. Phenotypic plasticity can also influence sclerite morphology across different environments, particularly size and shape (e.g., West et al. 1993; West 1997; Prada et al. 2008; Joseph et al. 2015). The interplay between genetic and environmental variation of these characters makes them difficult to use reliably in systematic studies.

Molecular approaches to classification have resulted in significant improvements in octocoral systematics (e.g., Sánchez et al. 2003), and much of this work has relied on mitochondrial genes. These markers have been useful in determining relationships between genera, but the slow rate of molecular evolution in Anthozoan mitochondrial DNA has limited their utility in distinguishing species (Shearer et al. 2002; McFadden et al. 2011). There are few single copy nuclear loci that reliably amplify across octocorals (McFadden et al. 2010) and those that are available may be low-resolution (e.g., Wirshing and Baker 2015), while other markers are multi-copy and show high intra-individual variation (e.g., Vollmer and Palumbi 2004). These obstacles have hindered researchers’ ability to draw boundaries between closely related taxa, but recent work indicates that high-throughput sequencing approaches can enhance our understanding of Anthozoan evolution and systematics, especially in recent radiations with low levels of genetic divergence among species (e.g., Pante et al. 2015; Combosch and Vollmer 2015; Quattrini et al. 2019; Erickson et al. 2021).

More recently, integrated approaches have proven fruitful in clarifying taxonomic and evolutionary relationships of octocorals using combinations of morphology, ecology, and molecular traits (e.g., McFadden et al. 2014; Wirshing and Baker 2015; Calixto-Botía and Sánchez 2017; Prada and Hellberg 2021). For example, McFadden et al. (2001) used molecular, morphologic, and reproductive traits to characterize reproductive trait evolution in Alcyonium and Prada et al. (2014) found host-symbiont specificity between symbiont phylotypes and distinct host lineages of Eunicea flexuosa. Aratake et al. (2012) similarly combined analyses of morphology, host genetics, symbiont genetics, and chemotypes to clarify relationships in the genus Sarcophyton, and Arrigoni et al. (2020) developed a robust taxonomic revision of the scleractinian genus Leptastrea using genomic (host and symbiont), morphometric, and distributional patterns. Studies such as these have greatly increased systematic resolution compared to analyses based on single lines of evidence, and they provide new perspectives in Anthozoan evolution.

An example of the challenges associated with octocoral species delimitation occurs with Plexaura homomalla (Esper 1792), one of the most common and readily recognized Caribbean octocorals. Plexaura homomalla forms colonies with dark brown branches and white polyps that make the species particularly distinctive (Bayer 1961, Fig. 1a). The species is widely distributed in forereef to backreef environments of the Caribbean and western Atlantic (Bayer 1961). The congener P. kükenthali Moser 1921 (the International Code of Zoological Nomenclature no longer recognizes diacritic marks in species names, therefore the accepted species name is P. kuekenthali, although we will use the original spelling in the species description by Moser) was originally described as a separate species on the basis of its bushier colonies, thinner branches (up to 2.5 cm vs. 4.5 cm diameter for P. homomalla), and light brown to tan coloration (Fig. 1b). Bayer designated the previously described P. kükenthali as a junior synonym of P. homomalla, recognizing the two groups as distinct morphotypes of P. homomalla due to the similarities in sclerite morphology (Bayer 1961, Fig. 1c, d). Plexaura kükenthali is distributed throughout the Caribbean and the Florida Keys (Bayer 1961), but has not been reported from Bermuda. The two species have been recorded in sympatry, for example, in Belize (Lasker and Coffroth 1983), Colombia (Velásquez and Sánchez 2015), Cuba (Espinosa et al. 2010), and Jamaica (Kinzie 1973). García-Parrado and Alcolado (1998) resurrected P. kükenthali, based primarily on the differences in branch thickness and color in conjunction with differences in the depth distribution of the two forms on Cuban reefs. Despite this revision, many researchers have continued to follow Bayer’s classification scheme, possibly because the publication associated with the taxonomic revision is not widely available (i.e., absent from international repositories). However, the reclassification is recognized in the World Registry of Marine Species (Cordeiro et al. 2020).

Fig. 1
figure 1

Colonies of Plexaura homomalla from the Florida Keys a and P. kükenthali from St. John b. Sclerites (spindles and clubs) isolated from P. homomalla c and P. kükenthali d

In this study, we coupled morphological and reproductive data with molecular phylogenetics, genomic approaches, and symbiont diversity to revisit the status of P. homomalla and P. kükenthali. Using these species, we tested whether a multifaceted approach could clarify a taxonomic relationship where data from any single source are ambiguous. In particular we addressed: (1) differences in sclerite morphology between the two closely related groups, (2) reproductive status and timing of the two species using field observations and dissection of collected colony samples, (3) genetic differences across a wide array of molecular markers (i.e., mitochondrial, nuclear, and transcriptomes), and (4) differences in Symbiodiniaceae symbionts hosted by the two groups. This study highlights the importance of an integrative approach to octocoral biology, including the application of genomic resources, to understand the relationships between closely related species of ecologically relevant cnidarians.

Materials and methods

Sampling

Samples of P. homomalla and P. kükenthali were collected between 2013 and 2019 from three locations at different depths: St. John, US Virgin Islands, the Florida Keys, and Puerto Rico (Table 1, Appendix 1 in Supplementary Information). Samples were identified based on branch thickness and color following Bayer (1961) and Sánchez and Wirshing (2005). Although P. kükenthali is reported to occur in the Florida Keys (Bayer 1961), we did not find colonies after searching numerous patch reefs in the Middle Keys; therefore, only P. homomalla is represented in the Florida Keys collections. A total of 44 P. homomalla and 42 P. kükenthali colonies were sampled for genetic and morphological analysis. During 2018 and 2019, small branch samples were collected from 10 to 15 colonies of each species on seven occasions from St. John to assess reproductive status. We did not note pronounced differences among colonies collected at the depths in this study; the importance of depth varies between octocoral lineages (see Calixto-Botía and Sánchez 2017 and Prada and Hellberg 2021 for examples in Antillogorgia and Eunicea, respectively).

Table 1 Listing of sample collection sites, number of colonies sampled, depth, and sampling year(s) for colonies used in this study

Sclerite analysis

Sclerites from each colony were isolated by digesting 2–4 mm branch fragments taken several centimeters below the tip of branches with commercial laundry bleach, rinsing the sclerites with water followed by ethanol, and then mounting sclerites onto glass slides using Permount. Plexaura homomalla and P. kükenthali, like all Plexaura spp., contain a variety of sclerite types. In discussing the two species, Bayer noted differences in size of the clubs, and Lasker et al. (1996), in their description of the congener P. kuna, reported small differences in spindle lengths and widths of the two forms of P. homomalla. Thus, our analysis focused on clubs and spindles (Fig. 1c, d). Ten randomly selected spindles and clubs from each of 29 P. homomalla (ten from the Florida Keys, 16 from St. John, three from Puerto Rico) and 27 P. kükenthali (21 from St. John, six from Puerto Rico) colonies were measured for tip-to-tip length and maximum width (excluding tubercles) using a Wild compound microscope at 200 × magnification with an eyepiece micrometer or an Olympus BX51 compound microscope and Infinity camera at 200 × magnification. Comparisons of length, width, and length/width (L/W) ratios were made with generalized linear models (SPSS ver. 26) using a linear model and the natural log of lengths and widths and the square root of the L/W as the response variable. Colony averages of lengths and widths for both sclerite types (i.e., four predictor values per colony) were used in a stepwise discriminant analysis with validation (SPSS ver. 26) to determine if those traits alone could differentiate the species.

Length and width provide a simple characterization of the complex morphology of sclerites, whereas elliptical Fourier analyses provide a more detailed description of the shape, which also enables quantitative comparisons of form within and between colonies and species (Carlo et al. 2011; Joseph et al. 2015). Sclerites were chosen by scanning the slides and selecting sclerites that were clearly clubs or spindles and had minimal contact with other sclerites. Multi-focus composite images were created using photographs taken at ten different foci for each sclerite. A 50 μm scale was transposed onto each composite image, and the images were processed as in Carlo et al. (2011). Elliptical Fourier descriptors (EFDs) for clubs and spindles were independently determined using SHAPE ver. 1.3 (Iwata and Ukai 2002). EFDs were developed using 20 harmonics for 65 spindles from 18 P. homomalla colonies (eight from St. John, eight from Florida Keys and two from Puerto Rico) and for 46 spindles from 14 P. kükenthali colonies (all from St. John). EFDs were generated as above for 40 clubs from 17 P. homomalla colonies (eight from St. John, nine from Florida Keys) and 20 clubs from ten P. kükenthali colonies (all from St. John). We generated the EFDs both with and without normalizing by the first harmonic. The dimensions of EFDs were then summarized in a principal components analysis (PCA). The two species were compared with analysis of variance of the principal component (PC) scores followed by discriminant analyses (SPSS ver. 26) based on both the PCs and the EFDs to determine if elliptical Fourier analysis could differentiate the species. Discriminant analyses were repeated using EFDs that were averaged for each colony.

Reproductive analysis

Studies of the two species (Goldberg and Hamilton 1974; Behety-Gonzalez and Guardiola 1979; Fitzsimmons-Sosa et al. 2004) suggest peak reproduction from June through August. In the summers of 2018 and 2019, samples of P. homomalla and P. kükenthali branches were randomly collected from colonies that were at least 30 cm in height at 7–10 m depth in South Haulover Bay, St. John, US Virgin Islands. Collections were made at intervals of two to three weeks. A single 5 cm branch fragment was collected from 15 colonies of each species (with the exception of June 20, 2018, when only 11 P. kükenthali were sampled). Samples were stored in 70% ethanol.

Reproductive status was assessed from dissections of 15 non-adjacent polyps from each sample. Polyps were randomly selected for dissection starting 1 cm below the tip, to avoid rapidly growing areas that are not reproductive (Gutiérrez-Rodríguez and Lasker 2004). Each polyp was carefully dissected, the number of eggs or spermaries counted, and their diameters measured with an eyepiece micrometer using a Wild stereomicroscope at 25 × magnification. Gonad diameters were used to estimate volume, assuming gonads as perfect spheres. As it is difficult to differentiate small, developing spermaries and eggs, we only determined the sex of the colony if the gonads were > 200 µm in diameter. Reproductive status is reported as the average volume of the gonads per polyp in individual colonies, as well as the number of mature eggs per each sample. Histology of a subset of samples was conducted to assess reproductive status in greater detail and to verify the sex of samples that had been identified in dissections.

In an attempt to observe spawning, night dives (19:00–21:00) on a fringing reef in Lameshur Bay at 4.5–9 m depth were conducted three to five days after the full moon in July 2016 and July 2017. Colonies of both species were located and mapped prior to the night dives and the colonies were then repeatedly surveyed during each dive. Plexaura kükenthali released planulae in 2016 and planulae were collected from the water column immediately above the colony using 50 mL syringes. In July 2019, branches from 12 P. homomalla colonies that had visible gonads were collected at South Haulover Bay, transported to and maintained in flowing seawater at the Virgin Islands Environmental Research Station (VIERS). Colonies of P. kükenthali at the South Haulover site did not have visible gonads at that time and were not collected other than for gonad analysis. In June 2021, branches from 7 P. kükenthali colonies that had gonads were collected and maintained in flowing seawater at VIERS.

DNA extraction and PCR

For single locus analyses, DNA was extracted from 0.25 cm branch sections of P. homomorphically and P. kükenthali using the 2X CTAB method (Coffroth et al. 1992) or the Qiagen DNEasy Blood and Tissue Kit following the manufacturer’s protocol (Qiagen, Venlo, Netherlands). Extracted DNA was diluted to five to ten ng/μL and kept at − 20 °C until amplified via PCR.

The mitochondrial genes MutS and ND2 and the nuclear introns EF1a, CT2, and I3P were amplified using PCR following Sánchez et al. (2003) and Prada and Hellberg (2013), respectively. Amplicons were sequenced by TACGEN (Richmond, California, USA) via Sanger sequencing. Each locus was aligned using the MAFFT ver. 7.833 plugin (Katoh and Standley 2013) in Geneious Prime 2020.1.2 (https://www.geneious.com). We constructed phylogenetic trees for each gene using IQTREE2 ver. 2.1.2 with 1000 ultrafast bootstraps (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017; Hoang et al. 2018) and MrBayes ver. 3.2.6 (Ronquist et al. 2012). Phylogenies were built from a concatenation of all five genes for samples with at least four of the five regions sequenced.

Transcriptome isolation and sequencing

Samples used in the transcriptome analysis were preserved in at least a 1:4 volume solution of RNAlater (Thermo Fisher Scientific, Waltham, MA) at 6 °C for twelve hours, RNAlater replaced, and stored at − 80 °C. mRNA was extracted from ~ 25 mg of tissue from each colony using the NEB mRNA magnetic isolation kit (New England BioLabs, Ipswich, MA) following the manufacturer’s protocol. Library preparation was done with the NEBNext RNA Library Preparation Kit, following the manufacturer’s recommendations, performed at half volume. Libraries were examined with an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA) and quantified with a fluorescent plate reader. To avoid lane effects, all samples were barcoded with unique identifiers, pooled, and sequenced across six lanes of Illumina HiSeq 2000 (100 PE) at the Vincent J. Coates Genomics Sequencing Laboratory (GSL), University of California, Berkeley.

Transcriptome assembly and SNP analysis

We retained sequences that had a score of Q30 or higher, and the first six bases were trimmed to remove adaptors using Trimmomatic ver. 0.33 (Bolger et al. 2014). Sequences were examined for quality and successful adaptor removal using Fast-QC ver. 0.10.1 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc). Prior to assembly, bacterial and algal contamination were removed with DeconSeq-graph ver. 0.1 (Schmieder and Edwards 2011), using BWA ver. 0.5.9-r16 (Li and Durbin 2009). NCBI’s bacterial rRNA submissions, Breviolum minutum genome (Shoguchi et al. 2013), and Symbiodiniaceae transcriptomes (Bayer et al. 2012) were used as the contaminant database.

Clean reads from all collected individuals were used to assemble the transcriptome for each species using Trinity ver. 2.1.1 (Grabherr et al. 2011) with normalization and Trimmomatic flags. Both de novo transcriptomes were annotated using a three-step approach using blastx with e-value < 10–5 and > 100 bp sequence overlap. First, we annotated transcripts with two Cnidarians (Nematostella vectensis, Thelohanellus kitauei) available on Ensembl (release 47). Transcripts without hits were then annotated with the complete Ensembl Metazoa database and the Cnidarians in UniprotKB (SwissProt and Trembl; The UniProt Consortium 2019). We identified the corresponding gene ontology (GO) terms for each annotated transcript from blastx results. We assessed the completeness of the assembled transcriptomes using BUSCO ver. 3, with the eukaryote and metazoan odb 9 databases (Simão et al. 2015). From these results, we chose to use P. homomalla as our mapping reference, given that it had higher BUSCO completeness.

To call single-nucleotide polymorphisms (SNPs), reads from all samples were mapped to the P. homomalla transcriptome using Bowtie2 ver. 2.3.5.1 (Langmead and Salzberg 2012) with default parameters for single-end reads. Resulting files were sorted and indexed using SAMtools (Li et al. 2009), and SNPs were identified using the mpileup and call commands in BCFtools. SNPs were filtered using vcfutils.pl in SAMtools, to exclude variants that had < 30X coverage. VCFtools was used to remove sites with missing data and sites within 5,000 bp of each other to avoid biases due to linkage (Danecek et al. 2011). The final SNP dataset was analyzed using the R packages Adegenet ver. 2.1.2 (Jombart 2008) and Hierfstat ver. 0.4–22 (Goudet 2005) in R ver. 3.6.1 (R Core Team 2019). For phylogenetic reconstruction, all SNPs were filtered using default parameters in SNPhylo ver. 3.695 (Lee et al. 2014) and aligned with MUSCLE ver. 3.8.31 (Edger 2004). Phylogenies were built as above. We used VCFtools and the script vcf_tab_to_fasta_alignment.pl (Bergey 2012) to extract an alignment of SNPs with Weir–Cockerham FST values > 0.7 (henceforth “high FST SNPs”), which were then used to build a second phylogeny. We performed a GO term analysis to determine if high FST SNPs were enriched for specific functions using the R package topGO (Alexa and Rahnenfuhrer 2019). Reads from each colony were also mapped (as above) to the mitochondrial genome of Muricea crassa (NC029697); SNPs were called following the same pipeline for the transcriptomes with the exception that SNPs were not filtered based on location. IQTREE2 ver. 2.1.2 with 1000 ultrafast bootstraps (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017; Hoang et al. 2018) was used to construct a phylogeny using the mitochondrial SNPs.

Symbiont genotyping

To assess patterns of host-symbiont (Symbiodiniaceae) specificity, we used fragment length analysis of the chloroplast 23S-rDNA hypervariable region to identify symbionts to the genus level following Santos et al. (2003). All symbionts detected were in the genus Breviolum (former Clade B). We then genotyped colonies at five Breviolum microsatellite loci: B7Sym4, B7Sym15, B7Sym34, B7Sym36 (Pettay and LaJeunesse 2007) and GV2_100 (Andras et al. 2009); primers and amplification conditions followed published protocols with annealing temperatures 54–57 °C. Amplicons were tagged with fluorescently labeled primers and visualized following Santos et al. (2003) and then scored by eye against 50–350 bp size standards (LI-COR Biotechnology Division). In some colonies, two alleles were detected at more than one locus. Given that symbionts are haploid in the vegetative state (Santos and Coffroth 2003), these were assumed to represent different symbiont genotypes. As it was not possible to assign a single multilocus genotype to these samples, we developed a curtailed dataset, excluding colonies where two alleles were detected at more than one locus (curtailed data = 92.7% of total dataset) and coded alleles as absence/presence data. We assessed population structure of both the original and curtailed datasets using RST and ΦPT in Arlequin ver. 3.5.2.2 (Excoffier and Lischer 2010) with 10,100 permutations and GenAlEx ver. 6.51b2 (Peakall and Smouse 2012) with 9,999 permutations. We corrected for multiple comparisons using Bonferroni correction. We found no significant differences in these metrics between the total and curtailed datasets and performed the remaining analyses with the curtailed dataset. Bayesian clustering analysis was conducted in STRUCTURE ver. 2.3.4 (Pritchard et al. 2000) with 100,000 burnins and 100,000 replicates. We assessed the best K based with the Evanno method (Evanno et al. 2005) in Structure Harvester (Earl and vonHoldt 2012). STRUCTURE results were fed to CLUMPP ver. 1.1.2 (Jakobsson and Rosenberg 2007) and visualized in DISTRUCT ver. 1.1 (Rosenberg 2004).

To reconstruct phylogenetic relationships of symbionts within different populations, the flanking region of the B7Sym15 microsatellite was amplified as above and sequenced for a subset of the samples. Sequences for Breviolum dendrogyrum (MH727217), B. endromadracis (KT149354), B. faviinorum (MH727210), B. meandrinium (MH727212), B. minutum (JX263427), and B. psygmophilum (JN602461) were included to provide phylogenetic context. Sequences were aligned and a maximum likelihood tree was constructed following the methods outlined above.

Results

Sclerite morphology

Sclerite dimensions overlapped substantially between P. homomalla and P. kükenthali, but there were significant differences between species and collection sites and among the colonies within each species and site (Table 2). Plexaura homomalla sclerites differed in club length and width between the three sites: clubs from the Florida Keys were longer and wider than those from either St. John or Puerto Rico (GLM, comparison of marginal means), but had the same L/W ratio (Table 2). Plexaura homomalla spindles from the three locations were similar in size and L/W. For P. kükenthali, spindle length was greater among the St. John samples. As P. kükenthali was absent in Florida surveys, interspecific comparisons were restricted to the colonies from St. John and Puerto Rico, where we had sympatric collections. Interspecific differences were present in club length and width and spindle width and L/W ratio.

Table 2 Morphological traits of Plexaura homomalla and P. kükenthali sclerites. Ten of each type of sclerite were measured from each colony. Values are means (standard deviation). Significant statistical tests are noted in bold

Linear regression of length and width identified significant heterogeneity of slopes between the two species (ANOVA, slope x species interaction, F1,466 = 4.09, p = 0.044; Fig. 2a, b). Plexaura kükenthali spindles were narrower than P. homomalla spindles and L/W of the spindles differed between the two species. Discriminant analysis based on length and width of individual sclerites was able to correctly identify the source species of 64.5% of the clubs and 71.1% of the spindles (SPSS ver. 26; discriminant analysis). When colony means of sclerite lengths and widths were used, the differences between the species were clearer (Fig. 2c, d) and discriminant analysis based on average lengths and widths using both sclerite types in a single analysis correctly identified 83.9% of the colonies.

Fig. 2
figure 2

Linear regressions of length and width measures of individual sclerites of Plexaura homomalla (blue) and P. kükenthali (red) clubs a and spindles b. Length and width measures of Plexaura homomalla (blue) and P. kükenthali (red) clubs c and spindles d with each point representing the average value for a single colony

Inclusion of detailed shape information did not provide substantially more information than sclerite lengths and widths alone (Fig. 3). When EFDs of spindles were normalized to the first harmonic, the discriminant analysis found that PC3 differed between the two species (SPSS ver. 26; discriminant analysis, Wilks-Lambda, p < 0.001) and identified the species of the spindles in 78% of the cases in the cross-validated groups procedure implemented in SPSS ver. 26. Limiting the analysis to only the St. John samples reduced identifications to 65% of the cases. Reintroducing sclerite size into the analysis by including length and width did not change the proportion of sclerites correctly identified (78%). Discriminant functions for individual sclerites included PC 1 and 3 and spindle width. Among cases in which there was more than a single sclerite from a colony, the species identification based on the call for the majority of the sclerites was correct for 87% of the colonies, and repeating the discriminant analyses using colony means instead of individual sclerites led to the correct identification of 87% of the colonies. Analysis of clubs yielded results that were similar to those for spindles. EFDs that were normalized on the first harmonic identified significant differences (Wilks-Lambda, p < 0.001) between the species in PC1 scores. The discriminant analysis identified the source of 65% of the clubs, 70% for P. kükenthali and 62.5% for P. homomalla. Inclusion of club length and width data in the discriminant analysis led to the identification of 75% of the clubs, 85% for P. kükenthali, and 70% for P. homomalla. Exclusion of the Florida Keys P. homomalla colonies reduced the proportion of correct calls to 50% when only the EFD data were included in the discriminant analysis and reduced the correct calls to 63% when length and width data were included. In cases where more than a single club was analyzed per colony, the majority of the clubs were correctly identified for 58% of the colonies. Discriminant analyses based on colony means, which could combine both clubs and spindle data, increased the proportion of correct identifications to 93%.

Fig. 3
figure 3

Principal component (PC) analyses for a clubs and b spindles of elliptical Fourier descriptors for Plexaura homomalla (blue) and P. kükenthali (red). PC1 and PC2 describe 31.0% and 15.8% of the total variance in club shape, respectively. PC1 and PC3 describe 23.8% and 11.1% of the total variance of spindle shape, respectively

Reproduction

The most striking difference between the two species was the observations of broadcast spawning of eggs and sperm by P. homomalla and the release of planulae by P. kükenthali. Release of planulae by four P. kükenthali colonies was observed in Lameshur Bay after sunset approximately 20:00 h on July 22–23, 2016 (Lunar Days 18–19, 3–4 days after the full moon). Planulae that were collected from the P. kükenthali colonies settled and metamorphosed on ceramic tiles several days after being collected. Plexaura homomalla was not observed spawning at that time. Neither spawning nor planulation was observed during the 2017 dives. A single colony of P. kükenthali maintained in the seawater system at VIERS released planulae on the evening of June 29, 2021, 5 d after the full moon.

Broadcast spawning of eggs and sperm was observed among P. homomalla colonies that were maintained in sea tables at VIERS on July 19–22, 2019 (lunar days 17–20, 2–5 days after the full moon). In sea tables, P. homomalla started releasing eggs between 21:00 and 21:30, with different colonies (males and females) starting at different times over a 60 min period. Release continued for approximately 2 h. Eggs collected during the spawning events developed into planulae which successfully settled (see Tonra et al. 2021 for more details).

Samples of colonies from South Haulover Bay did not exhibit reproductive synchrony; many of the colonies that were collected did not have visible gonad, and in many cases visible gonads were not mature (Fig. 4). Among P. homomalla colonies that developed identifiable gonads, we observed an increase in the number of large eggs (> 450 µm) and the presence of identifiable male colonies in July of both years. This was followed by a decrease in the number of large eggs in most female colonies and the disappearance of identifiable male colonies in August. Among P. kükenthali colonies, a small number of samples contained large eggs in June and early July 2018. Males were not observed in dissections of the P. kükenthali samples, but we did not conduct histology on all of the samples and that result can only be interpreted as an absence of mature males.

Fig. 4
figure 4

Average gonad volume per polyp (with 95% confidence interval) of 15 polyps from each colony collected on St. John, Virgin Islands. Not shown are zero values for 40 samples in which no gonad tissue was visible at 25X magnification. Solid vertical line separates species, dark vertical line separates years

Only 17 of the 60 female P. homomalla and 11 of the 36 female P. kükenthali colonies sampled contained large eggs. Despite the small number of reproductive colonies in our samples, there were distinct differences between the two species. Numbers of eggs found in each polyp differed between dates and between species (Log-linear analysis: eggs per polyp by collection date by species, p < 0.001; eggs per polyp by collection date, p < 0.001; eggs per polyp by species, p < 0.001). If only large/mature eggs are considered, there was a significant difference between the two species and between collection times (χ2 analysis, p < 0.001).

Molecular approaches

Sanger sequencing of individual loci resulted in a matrix of 1,337 bp from the mitochondrial genes MutS and ND2, and 1,177 bp from the nuclear introns EF1a, CT2, and I3P. For the samples that had at least four regions sequenced, the matrix was composed of 23 samples (nine P. kükenthali, 14 P. homomalla). For this subset of samples, nearly all sites (2,479 bp, 98.6%) were invariant and only 22 sites (0.88%) were parsimony informative, with 12.12% missing data. The maximum likelihood tree from the concatenated gene set (log-likelihood -3774.36, Fig. 5a), had moderate support for P. homomalla and P. kükenthali clades (80 ML bootstrap support [BS], 0.99 posterior probability [PP]), while individual gene trees were largely unresolved or contained branches with low support (Fig. S1). The only exception was EF1a (Fig S1e), which provided moderate support for the differentiation of the two species.

Fig. 5
figure 5

Phylogenies constructed from a concatenated mitochondrial genes (MutS, ND2) and nuclear introns (EF1a, CT2, I3P) and b mitochondrial SNPs. Colony identification numbers provided in parentheses. Support values given as maximum likelihood bootstraps/posterior probabilities. Support values < 75 BS not shown

A total of 260 mitochondrial SNPs were retained after filtering for depth; 26 of these sites were parsimony informative with a total of 3.32% missing data. In the phylogeny (log-likelihood − 521.17, Fig. 5b), there was strong support for the differentiation of P. homomalla and P. kükenthali (100 BS, 1 PP). Interestingly, one P. homomalla colony (HF3) nested within the P. kükenthali clade (Fig. 5b). By all phenotypic aspects, HF3 is unequivocally P. homomalla (Fig. S8); it was therefore surprising to see this individual nested within P. kükenthali.

Transcriptome sequencing resulted in 73,197,074 reads after filtering of contaminants (average 9,149,634 reads per individual, SE ± 615,899 reads). The de novo transcriptome assembly of P. homomalla resulted in 66,781 contigs, 23,535 (35.24%) of which were annotated: 17,911 from Ensembl Cnidarians, 4,137 from Ensembl Metazoa, and 1,487 from UniprotKB Cnidarians. The P. kükenthali transcriptome resulted in 130,138 contigs, 41,228 (31.68%) of which were annotated: 25,645 from Ensembl Cnidarians, 7,433 from Ensembl Metazoa, and 8,150 from Uniprot KB Cnidarians (Table 3, Fig. S2-4). We identified 26,779 SNPs, which had an average FST of 0.2196 between the two species. Pairwise comparisons of individual genes allowed us to identify 319 SNPs with Weir–Cockerham FST > 0.7 (ca. 1% of identified SNPs). In total, there were 15,956 genes with successful GO annotation in the P. homomalla transcriptome (4,865 biological process [BP], 5,162 cellular component [CC], 5,929 molecular function [MF]). Within the high FST genes, there were 15 GO terms that were significantly enriched (nine BP, five CC, one MF) using Fisher’s exact test (p < 0.01) (Table S2; Fig. S5–7). No terms were significantly enriched after correcting for multiple comparisons using the Benjamini and Hochberg method. Enriched GO terms are related to metabolic processes (e.g., GO:1,901,575 organic substance catabolic process), with genes such as Derlin and proteasome subunits (e.g., 26S proteasome non-ATPase regulatory subunit 13) associated with enriched terms. The protein products of these genes are involved in the endoplasmic reticulum-associated protein degradation (ERAD) pathway and other catabolic processes. While further work is needed to validate these differences, our transcriptome-wide analysis here suggests that genes involved in protein catabolism may be some of the first to differentiate in this system.

Table 3 De novo transcriptome assembly statistics for Plexaura homomalla and P. kükenthali

A total of 16,859 SNPs were retained after filtering by SNPhylo (ca. 63% of all 26,779 SNPs). Nearly all sites (16,410, 97.34%) were parsimony informative with only 449 (1.67%) invariant sites. The maximum likelihood (log-likelihood -277,075.64) and Bayesian trees from the full SNP dataset resulted in the same topology as the mitochondrial genome tree (Figs. 5b, 6a). In the phylogeny constructed from the full SNP dataset, one colony identified as P. homomalla (HF3) nested within P. kükenthali (Fig. 6a, Fig. S8), whereas the same individual nested in the P. homomalla clade in the high FST SNP phylogeny in both maximum likelihood (log-likelihood − 2729.44) and Bayesian analyses with strong support (100 BS, 1 PP; Fig. 6b).

Fig. 6
figure 6

Phylogenies constructed from the entire SNP dataset a and high FST SNPs b. Support values given as maximum likelihood bootstraps/posterior probability. Support values < 75 BS not shown. Admixture analysis with Ohana with all SNPs c and high FST SNPs d for optimal K = 2

In order to assess the possibility of hybridization between the two species, as suggested by the mitochondrial genome and SNP phylogenies for the putative hybrid HF3, we performed an admixture analysis with Ohana (Cheng et al. 2017) using default parameters. The analysis was run using both the entire SNP dataset as well as the high FST SNPs (Fig. 6c, d). Both analyses resulted in an optimal K = 2 (all SNPs log-likelihood − 129,558.17; high FST SNPs log-likelihood -807.28) but the analysis with the high FST SNPs revealed that HF3 could be a potential hybrid individual (Figs. 5b, 6a, S8). While all SNPs would likely be affected by admixture, the statistical power to detect admixture is much stronger for highly differentiated SNPs than when using the complete dataset (e.g., Ma et al. 2019). In other words, recent hybrids are easier to spot in datasets where SNPs are highly differentiated. Given the placement of the mitochondrial genome of HF3 (Fig. 5b), if this individual is a hybrid, as suggested, the maternal parent is P. kükenthali and the paternal parent is P. homomalla. The possibility of cryptic hybrids could complicate taxonomic classification, and further highlights the importance of using genome-wide analyses for species identification.

Symbiont specificity

Of the 77 colonies analyzed for symbiont genotype, six harbored multiple alleles at more than one locus (three P. kükenthali from St. John; one P. kükenthali and two P. homomalla from Puerto Rico) and were excluded to create a curtailed dataset. Population pairwise RST and ΦPT were significantly different for all populations for both datasets (full and curtailed, Tables S4 & S5), but were not significantly different between the full and curtailed datasets. The removal of those six samples (to generate the curtailed dataset) only reduced the total number of alleles by three (full = 36, curtailed = 33). Thus, the curtailed dataset was used to compare symbionts within P. homomalla and P. kükenthali.

Among the five analyzed symbiont microsatellite loci, between three and ten alleles were detected at each locus; a total of 33 alleles were recovered (Table S6). Of these, eight alleles were unique to symbionts from P. homomalla colonies and ten were unique to symbionts from P. kükenthali colonies (Table S6). A total of 55 unique symbiont genotypes were identified among the 71 colonies of P. homomalla and P. kükenthali with the number of genotypes in a given region ranging from five to twenty-one (Table S7). Genotypes were not shared between regions and were rarely shared between host species. Other than four cases, genotypes were not shared between reefs or species.

Bayesian clustering analysis identified two symbiont populations (K = 2, Fig. 7a). Plexaura kükenthali colonies harbored a single symbiont population across all regions (P. kükenthali phylotype) (Fig. 7a). A second symbiont population, which differed from the P. kükenthali phylotype, was found in P. homomalla colonies from Puerto Rico and Florida (P. homomalla phylotype). However, in St. John while most P. homomalla colonies harbored the P. homomalla phylotype, the P. kükenthali symbiont phylotype were found in eight P. homomalla colonies (Fig. 7a).

Fig. 7
figure 7

a Results of STRUCTURE analysis for symbionts genotyped at five microsatellite loci for K = 2. b Maximum likelihood phylogeny generated from the B7Sym15 microsatellite flanking region for several P. homomalla and P. kükenthali symbionts and additional Breviolum species. Symbionts from P. homomalla colonies that fall in the P. kükenthali phylotype are bolded and denoted with asterisks. Branch supports were generated from 1000 ultrafast bootstrap replicates; only BS > 75 are shown

Phylogenetic relationships among symbiont genotypes provide robust support for two clades, representing the two host species (Fig. 7b). One group included the majority of the P. homomalla symbionts (P. homomalla phylotype), including those in the host HF3 that had grouped with P. kükenthali colonies based on host genetic data (Figs. 5 and 6a). Symbionts in the P. kükenthali phylotype formed a second clade that included symbionts from two P. homomalla colonies. It is noteworthy that symbionts from these two P. homomalla grouped with St. John P. kükenthali symbionts in the STRUCTURE analysis (Fig. 7a) and one of them harbored the same symbiont genotype as a P. kükenthali colony.

Discussion

We used a multifaceted approach to clarify the species boundaries of two widely recognizable Caribbean octocorals, whose taxonomic history has been uncertain for decades. Although García-Parrado and Alcolado (1998) officially resurrected P. kükenthali as a distinct species over two decades ago, they based their conclusion on depth distribution patterns of the two groups, without analyses of morphological, reproductive, or molecular divergence. Here, multiple lines of evidence reinforce the designation of these two groups as separate species, highlighting the utility of an integrated approach of multiple data sources in octocoral systematics. Furthermore, this study presents the assembly and annotation of two new reference transcriptomes of Caribbean octocorals.

Sclerite morphology

Traditional octocoral taxonomy relies heavily on sclerite morphology. However, these traits exhibit high levels of variability in form (Kim et al. 2004; Joseph et al. 2015). Our morphometric analyses identified differences between the two species in both the size (Table 2) and to a lesser extent shape of the clubs and spindles. The differences were evident when average lengths and widths of the two species were compared (Fig. 2c, d). Clubs of P. homomalla were thicker and more robust in appearance than those of P. kükenthali (Fig. 1c, d). However, neither EFDs nor length/width measurements unambiguously identified the source of any single sclerite. We found high levels of morphological variability of these characters, with individual colonies within a species significantly differing in length and width measures of both spindles and clubs (Table 2). The failure of the discriminant functions of both length/width measures and EFDs to fully distinguish between the two species highlights that morphological plasticity in these characters limits their utility in taxonomic studies. EFDs have been useful in analyses of other taxa (e.g., Carlo et al. 2011) but in this case quantifying sclerite form was not very different than the analyses based on length and width alone. There will be many cases in which size alone cannot differentiate species, but this study suggests the far more readily obtained size metrics should be analyzed first. It is also likely that EFDs are less effective in characterizing the shape of complex sclerites such as clubs, which exhibit little symmetry, beyond having a dominant axis. In these cases, positioning of the sclerite on the slide will dramatically affect the profile and mask similarities and differences between sclerites.

Reproductive cycles

Plexaura homomalla had been observed spawning in the field (Bastidas et al. 2005) and prior to that broadcast spawning was inferred by the absence of planulae within polyps during times when mature eggs and sperm were present (Goldberg and Hamilton 1974; Martin 1982; Fitzsimmons-Sosa et al. 2004). Behety-Gonzalez and Guardiola (1979) observed a similar cycle of gonad development in P. kükenthali. Thus, we were surprised to observe planulae being released by four P. kükenthali colonies in 2016. While differing reproductive systems are known within a single octocoral genus (McFadden et al. 2001), we are not aware of any reports of both broadcast spawning and planulation within a single octocoral species. We were unable to replicate this observation during 2019, but planulation of P. kükenthali was observed in 2021.

The clear pattern of peak reproduction during the summer that has been observed elsewhere for both species (Goldberg and Hamilton 1974; Behety-Gonzalez and Guardiola 1979; Martin 1982; Fitzsimmons-Sosa et al. 2004) was not evident in the samples from South Haulover Bay. The simultaneous sampling of both species revealed few mature colonies in both 2017 and 2018. Still, the greatest gonad content in P. homomalla was observed in July, whereas P. kükenthali gonad content was greatest in June (Fig. 4). It is possible that in both years we sampled either after a substantial number of colonies had already spawned, or months before spawning. Complete sampling across the entire summer and possibly year is needed to fully characterize the species’ reproductive behavior, but this study suggests there are differences in their reproductive cycles.

The overlap in the timing of reproduction may be sufficient to allow hybridization between the two species. Sperm from P. kükenthali and eggs and sperm from P. homomalla could be present in the water column at the same time, and this could allow for cross-fertilization, even if one of the species broods its eggs. Differences in gamete-reception proteins have been suggested to prevent hybridization during mass spawning (Babcock 1995). While other marine invertebrates such as sea urchins and abalone have become model systems for studying the role and evolution of bindin and lysin gamete recognition proteins (Swanson and Vacquier 2002), such receptor proteins remain poorly understood in Anthozoans. Future work may capitalize on the resources and data presented here to study the role of gamete recognition in octocorals with different (albeit overlapping) reproductive cycles and modes of reproduction.

Molecular insights

With the exception of Ef1a, gene trees constructed from single mitochondrial and nuclear loci did not provide enough resolution to confidently separate our two taxa (Fig. S1). However, the concatenation of two mitochondrial genes and three nuclear introns resolved two clades, corresponding to the two species (Fig. 5a). SNPs across the transcriptome showed differentiation between P. homomalla and P. kükenthali (FST = 0.2196). The most differentiated SNPs had FST values of ≥ 0.85 (210 out of 318 high FST SNPs), and 30 SNPs were within loci which had annotations at the gene level. Some of these genes encode proteasome subunits (26S proteasome regulatory subunit RPN9, proteasome subunit beta and alpha type), ubiquitin-associated proteins (e.g., E3 ubiquitin-protein transferase), transport and signaling proteins (e.g., COP9 signalosome complex subunit 6), and nucleotide replication and transcription proteins (e.g., DNA-directed RNA polymerases I and III subunit RPAC1). Interestingly, inositol-3-phosphate synthase (I3P) was one of these highly differentiated genes, although the gene tree generated from the I3P intron did not resolve the two species (Fig. S1d). Further investigation of the role these genes will be required to elucidate if they play a role in octocoral speciation.

We found a putative hybrid in our samples, indicating the potential for hybridization between the two species (Fig. 6). Hybridization is not uncommon in octocorals (McFadden and Hutchinson 2004; Quattrini et al. 2019; Erickson et al. 2021) and scleractinians (Willis et al. 1997, 2006). Although hybridization may blur and complicate species boundaries, it may also facilitate diversification on reef systems by leading to novel morphologies and ecological niches (e.g., Vollmer and Palumbi 2002). The possible overlap in the timing of reproduction of P. homomalla and P. kükenthali (Fig. 4) suggests that pre-zygotic isolation might be permeable in these groups. In mass-spawning corals, where many species spawn synchronously on the same reef, other isolating mechanisms such as differing gamete-reception proteins have been suggested to prevent hybridization (Willis et al. 1997). Although our admixture analysis suggests that P. homomalla and P. kükenthali may form hybrids (Fig. 6), our results show that hybridization might be sporadic as it has not resulted in the complete swamping of the morphological and genetic differences of the two groups. Therefore, other than the putative hybrid HF3, we suggest that gross morphology based on branch thickness and color remain reliable for field identification.

Symbiont specificity

Studies suggest that specific host-symbiont pairings are more common than previously thought, and species within the symbiont genus Breviolum have specific host groups (e.g., Prada et al. 2014; Lewis et al. 2019). Plexaura homomalla and P. kükenthali tended to harbor distinct Breviolum lineages, with symbiont populations in P. homomalla and P. kükenthali generally clustering according to host species rather than geography (Fig. 7). For example, P. homomalla symbionts from Puerto Rico clustered with individuals from Florida rather than clustering with P. kükenthali symbionts from Puerto Rico (Fig. 7). Furthermore, the presence of unique symbiont alleles for each host species suggests a lack of gene flow between the symbionts of P. homomalla and P. kükenthali (Table S6). Finally, phylogenetic analysis based on the flanking region of the B7Sym15 microsatellite locus also supports two host-specific symbiont phylotypes (Fig. 7b).

Despite these clear clustering patterns, further study is needed into the presence of shared symbiont phylotypes between the host species in St. John. It is unlikely that changes in symbiont types due to switching or shuffling (sensu Baker 2003) would be responsible for these shared phylotypes as bleaching is very rare in Plexaura (Lasker 2003; Prada et al. 2010). It is possible that symbiont loss does occur and is not observed due to the host coloration; however, in most cases where symbiont genotypes were followed across bleaching events or stress in octocorals, the symbionts did not change (Goulet and Coffroth 2003, Kirk et al. 2005, Coffroth unpubl. data; but see Lewis and Coffroth 2004 for an exception). Expanding the number of loci examined might differentiate these, however, sequence analysis of the flanking region of the B7Sym15 microsatellite suggests a close relationship between some St. John P. homomalla symbionts and the P. kükenthali symbiont phylotype (Fig. 7b). Alternatively, local microenvironments might have selected for this pairing as all P. homomalla colonies that harbored the P. kükenthali symbiont phylotype were collected in St. John, USVI; however, all P. homomalla with the P. kükenthali symbiont phylotype were collected from the same reef and at the same depth where P. homomalla with the P. homomalla symbiont phylotype were also found.

Host specialization is recognized as an important factor in the speciation of these symbionts (e.g., LaJeunesse 2005; Finney et al. 2010; Thornhill et al. 2014) and is likely a result of selection for a given host habitat or resource, where symbionts distinguish between the host species, based on differences in the intracellular environments (Thornhill et al. 2014; Lewis et al. 2019). The factors that influenced the observed symbiont genotypes shared between the host species are unclear and an area for investigation. Regardless, the separation of P. homomalla and P. kükenthali symbionts by host species over a broad geographic range further delimits these two octocoral species.

Conclusions

Using an integrative approach that included morphological analyses, reproductive observations, molecular evaluations, and symbiont diversity, we were able to reinforce that P. homomalla and P. kükenthali are, in fact, distinct species. While there are key external morphological characters that differentiate the two species, differences in reproductive timing and microscopic morphological traits can also be used to distinguish the species. Molecular techniques used to identify and delimit taxa have become commonplace in octocoral systematics and were an invaluable tool in this study. This work also presents two new transcriptomes of Caribbean octocorals that can be used to further the understanding of topics ranging from octocoral phylogenomics to fine-scale plexaurid population structure. Finally, it is important to emphasize the need to take a multifaceted approach to investigate octocoral systematics as individual data sources can fail to identify sufficient evidence to inform taxonomic classifications.