Evidence of genetic isolation between two Mediterranean morphotypes of Parazoanthus axinellae

Coralligenous assemblages are among the most species-rich and vulnerable habitats of the Mediterranean Sea. Nevertheless, data on connectivity patterns on species inhabiting these habitats, crucial to define management and protection priorities, are largely lacking. Moreover, unreliable species-level taxonomy can confound ecological studies and mislead management strategies. In the northwestern Mediterranean two Parazoanthus axinellae morphotypes differing in size, color and preferred substrate are found in sympatry. In this study, we used COI and ITS sequence polymorphism to assess (1) the genetic divergence between the two morphotypes, (2) their connectivity patterns and (3) their phylogenetic position within the Parazoanthidae. Specimens of P. axinellae were sampled in 11 locations along the northwestern Mediterranean; in 6 locations, samples of the two morphotypes were collected in sympatry. Small genetic diversity and structure were found within morphotypes, while marked and consistent differentiation was detected between them. Moreover, the less widespread morphotype appeared to be closer to Pacific species as P. juanfernandezii and P. elongatus. Our findings confirmed the limited knowledge on Parazoanthus species complex, and how this gap can have important implication for the conservation strategies of this widespread and valuable genus in the Mediterranean Sea.

Scientific RepoRtS | (2020) 10:13938 | https://doi.org/10.1038/s41598-020-70770-z www.nature.com/scientificreports/ The presence of species complexes in the sea can confound ecological studies and mislead management strategies. For example, connectivity patterns 12 may be erroneously inferred when cryptic species are ignored 13,14 , which may provide biased conclusions about the overall capacity of the studied populations to resist stress 15 . To disclose cryptic species before estimating gene flow and indeed evaluate connectivity among populations, we should (1) use more than one independent molecular marker, (2) sample individuals all along the distributional range of the investigated species, and (3) sample closely related taxonomic species to understand the spectrum of population to species divergence 13 .
To our knowledge, three species of Parazoanthus exist along the Mediterranean and North-East Atlantic coasts: Parazoanthus axinellae (Schmidt, 1862), P. anguicomus (Norman, 1869) and another Parazoanthidae yet to be described in the Macaronesian waters 16 . In the Mediterranean Sea, Parazoanthus axinellae (Cnidaria, Hexacorallia, Zoantharia) species complex 17 is one of the most common cnidarians inhabiting coralligenous assemblages. The simplicity of the zoanthid body plan makes morphology-based species identification quite challenging. Numerous morphological identification criteria have been used to identify phylogenetic signal ranging from colour, sphincter muscle anatomy 18 , tentacles number 19 , type and distribution of nematocysts 20 . In addition, differences in substrate preference, overall ecology 17 , biochemical profiles [21][22][23] and sequence divergence using various gene markers were also proposed as zoanthid identification criteria [24][25][26] . Original species descriptions were thereby based on different traits depending on the authors 27 . Due to the aforementioned lack of morphological characters to identify zoanthids, a comprehensive taxonomic identification of the P. axinellae species complex is still missing (but see 17 for a detailed description of P. axinellae species complex). Historically, there were four described subspecies of P. axinellae: P. axinellae mülleri Pax, 1957, common across the whole Mediterranean coastline; P. axinellae liguricus Pax, 1937, in the Ligurian Sea; P. axinellae brevitentacularis in the Gulf of Lion and P. axinellae adriaticus Pax, 1937 in the Adriatic Sea [28][29][30] . However, the morphological descriptions of these four subspecies are quite imprecise and overlapping, and clear diagnostic characters to discriminate them are still lacking. Recently, Ocana et al. 17 stated that these four subspecies can be grouped in two morphotypes based on their chemical profiles 31 and ecological requirement 17 . These two morphotypes live in sympatry in the northwestern Mediterranean 29 . One of them displays an elongated trunk, longer and thinner tentacles, and a light yellow color. This morphotype ("Morphotype 1" in 17 , "Slender" in 31 ; Fig. 1a) is found across the whole Mediterranean, including the Adriatic Sea, and could correspond to the descriptions of the two subspecies P. axinellae mülleri and P. axinellae adriaticus (Fig. 1a). A more pronounced orange color, and shorter and thicker trunk and tentacles characterize the second morphotype (Morphotype 2 in 17 , "Stocky" in 31 ; Fig. 1b). This form is restricted to the north-western Mediterranean and could correspond to P. axinellae brevitentacularis and P. axinellae liguricus (Fig. 1b). On top of their morphological differences, the "Slender" form is mainly characterized   32 and has never been observed on Axinella spp 17 .
In the present study, we aim to describe the patterns of diversity and connectivity for populations of the two Parazoanthus axinellae morphotypes along the Northwestern Mediterranean and the Adriatic Sea. Specifically, we use mitochondrial (COI) and nuclear (ITS) sequences polymorphism to asses: (1) the genetic diversity and structure among populations within and between morphotypes and (2) phylogeny and differentiation in the Parazoanthidae Family.

Results
Genetic diversity and structure among P. axinellae populations. A total of 322 sequences of 402 bp of COI were obtained (Table 1), and 5 polymorphic sites and 7 haplotypes were detected. Haplotype and nucleotide diversity were low within all the localities. The "Slender" and the "Stocky" morphotypes showed 3 (I, IV, VI) and 4 (II, III, V, VII) different haplotypes, respectively (Fig. 2, Supplementary Figure S1, Supplementary Table S1). The unrooted haplotype network exhibited two star-like patterns (haplogroups: "Slender" and "Stocky") connected by 1 mutational step. The two haplogroups were surrounded by private and low frequency haplotypes differing in 1 bp (Fig. 2, Supplementary Figure S1, Supplementary Table S1). Table 1. Locality, geographical coordinates, morphotype, depth, and code of Parazoanthus axinellae samples; number of individuals sequenced (n), haplotypes (H), haplotype (Hd) and nucleotide (π) diversity values per location (with their standard deviation-SD) for both molecular markers.

Locality
Coordinates Morphotype Depth Code n COI H COI Hd COI SD π COI SD n ITS H ITS Hd ITS SD π ITS SD www.nature.com/scientificreports/ As for the ITS fragment, a total of 20 haplotypes were obtained from 232 sequences (698 bp) and "Stocky" and "Slender" morphotypes were differentiated by a minimum number of 68 mutation steps (between haplotypes II and VI) (  Table S2). The two morphotypes showed different ITS haplotypes, with the exception of 9 individuals sampled in Olbia, West Tyrrhenian Sea (codes SRY and SRO) that shared the same haplotype with the morphologically distinct phenotype. In particular, 5 "Stocky" individuals at Olbia (code SRO) had the identical ITS sequence with the haplotype IV (mostly composed out of "Slender" individuals), and 4 specimens with the "Slender" phenotype (code SRY) shared the same sequence with haplotype I, composed mostly of "Stocky" Parazoanthus specimens (see Fig. 2, Supplementary Figure S1).
Genetic differentiation (pairwise Φ ST values) between populations ranged between 0 and 1 for both markers (See Supplementary Table S3, Supplementary Table S4 Figure S3). The samples from Olbia (SRY and SRO codes) showed a specific pattern of genetic differentiation, with ITS haplotypes being intermixed within morphotypes, as previously elaborated.

Figure 2.
Haplotype network analysis for COI (above) and ITS (below) gene markers in Parazoanthus axinellae species complex. Haplotype network was built in Rstudio version 1.1.453 (https ://rstud io.com/) with an R package pegas (https ://cran.r-proje ct.org/web/packa ges/pegas /index .html), using an infinite site model based on Hamming distances of DNA sequences. Each haplotype is presented as a separate pie chart, and size of the pie charts directly correlates to the number of individuals within the haplotype. Pie chart colors for geographic location were chosen to intuitively associate the reader to Parazoanthus axinellae morphotypes, with orange tones being attributed to the "Stocky" morphotype, and yellow and grey tones to the "Slender" morphotype. Distances between haplotypes correspond to the genetic differentiation observed for each of the markers, and each mutation between the haplotypes is shown as a hyphen. Alluvial diagrams were made using an online platform RAWGraphs (https ://rawgr aphs.io/) to further clarify haplotype-morphotype assignments, and final graph modifications were done in Inkscape version 0.92.4 (https ://inksc ape.org/it/). Scientific RepoRtS | (2020) 10:13938 | https://doi.org/10.1038/s41598-020-70770-z www.nature.com/scientificreports/ The Analysis of Molecular Variance (AMOVA) between groups of samples according to morphotype and biogeographic areas showed that for the "Slender" morphotype most variability is explained by differences between all populations (COI: 81.79%, P < 0.01; ITS: 68.87%, P < 0.01; Table 2). For the "Stocky" morphotype, the same is true according to COI, but ITS shows some variability related to the biogeographic pattern (Table 2).
AMOVA between the two groups defined by the morphotype attributed most of the variation to the differences between morphotypes for both markers (COI: 96.62%, P < 0.01; ITS: 89.52%, P < 0.01; Table 2).
Neutrality tests showed no departure from neutrality according to COI in both morphotypes, but a positive and significant departure according to Fu and Li test 33 for the ITS marker in both morphotypes (1.77, P < 0.02 for the "Slender"; and 2.28 P < 0.02 for the "Stocky"), which might indicate a decrease in population size and/or balancing selection acting in this marker. Mismatch distributions in both markers showed a bimodal distribution 34 Table S5, Supplementary  Table S6). These haplotypes were reconstructed from 344 (COI) and 275 (ITS) same-length sequences (474 and 734 nucleotides for COI and ITS genes) and were used for Maximum Likelihood (ML) treebuilding.
Likelihood ratio test (performed with modelTest function in phangorn R package), based on minimum AICc values, identified HKY + I and HKY + G + I models as best nucleotide evolution models for COI and ITS markers, respectively. The Maximum Likelihood trees were reconstructed for each marker using F81 nucleotide substitution model as this was the most similar model to HKY + I and HKY + G + I available in phangorn R package.
Although indels in ITS gene marker created large gaps, the overall network topology and the number of mutation steps between species did not vary greatly between the different MAFFT alignments performed (only one result is shown). Moreover, the phylogenetic trees with MAFFT alignment (without gap removal) and with MAFFT alignment + GBlock gap removal also give similar topology. We only show results for P. axinellae + Genbank ITS sequences aligned in MAFFT using L-INS-i algorithm with 1,000 iterations, and with Gblock gap removal parameters set to: b1 = 0.8, b2 = 0.9, b4 = 2, b5 = "h".
Bayesian Inference and maximum likelihood trees reflected the same topology for both COI and ITS genes but with different statistical support. It is interesting to note that various haplotypes in both markers shared the same sequence across more than one species. Particularly for COI gene, this was the case for haplotypes I (ANG, PAR and "Slender" P. axinellae), and V (ALI and "Stocky" P. axinellae), whereas for ITS marker this occurred only in haplotype VI (ELO and JUA). Overall "Stocky" and Slender" morphotypes cluster in two well-separated clades (although with low bootstrap support). In the "Stocky" clade all the "Stocky" morphotype haplotypes showed high similarity with P. elongatus (COI and ITS), and P. juanfernandezii (ITS), both species from the southeastern Pacific. The "Slender" morphotypes grouped with P. anguicomus (COI and ITS) from the eastern Atlantic, P. parasiticus (COI), P. swiftii (COI) both Caribbean species and P. aliceae (COI) from the Azores Archipelago (Fig. 3). Finally, it is interesting to note that despite clear genetic differences between the studies morphotypes, www.nature.com/scientificreports/ certain individuals sampled in Olbia, West Tyrrhenian Sea (codes SRY and SRO) shared the same sequence for ITS marker with the morphologically distinct phenotype (Fig. 3).

Discussion
In the present study we evidenced for the first time that (1) the two Mediterranean morphotypes of Parazoanthus axinellae differ more between them that in comparison to other allopatric zoanthid species and (2) that the "Slender" morphotype is characterized by a lower genetic differentiation among populations compared to the "Stocky" morphotype. COI and ITS evidenced clear genetic isolation between the two morphotypes. Differences in chemical profiles 31 and in morphological features 17 between the two morphotypes are consistent with our results, suggesting that the Parazoanthus axinellae species complex includes two different taxa. Both COI and ITS haplotype networks showed the presence of two highly abundant haplotypes surrounded by low frequency haplotypes that clearly discriminate the two morphotypes. This low variability can be explained by the geographical distribution of the morphotypes. In fact, the "Slender" morphotype is widespread all along the Mediterranean Sea while the "Stocky" morphotype 17 has a more restricted distribution. Although we cannot discard the existence of more populations of the "Stocky" morphotype in other unsampled areas of the Atlantic Ocean, information on the geographic distribution of this species in the Eastern Atlantic Ocean is scarce and refers only to the "Slender" morphotype (Boavida J., personal communication). The widespread distribution of the "Slender" morphotype can also explain the low genetic structuring observed among populations.
Despite the limited geographical distribution compared to the "Slender" morphotype, the "Stocky" morphotype shows genetic structure between northeastern Tyrrhenian populations (Portovenere, Portofino, and Giannutri) and the northwestern ones (Banyuls -sur -mer and Alassio) with a break around Portofino. Portofino area represents a barrier to gene flow for several species of sessile invertebrates (Corallium rubrum 36 ; Paramuricea clavata 37 ; Patella caerulea 38 ) related to the presence of a large-scale separation of currents that occurs in the region 39 . Biological and reproductive characteristics of the species can also contribute to explain the observed pattern. Pax and Muller 40 found hermaphroditic colonies of P. axinellae in the Adriatic Sea, while observed that colonies of P. axinellae from the British Isles (North Atlantic) are gonochoric 41 . Moreover, Previati et al. 42 observed that the "Stocky" morphotype collected in the Ligurian Sea, where the two morphotypes live in sympatry, exhibited both asexual and sexual reproduction (Cerrano C. personal communication). Although not much information is available on the ecological and biological traits of this morphotype, it seems to be defined by a lower dispersal capacity compared to the "Slender" morphotype in this area.
Sardinia samples (codes SRY and SRO) present an intermediate pattern of differentiation, with some samples from both "Slender" and "Stocky" morphs sharing the same sequence with the morphologically distinct phenotype. These identities on ITS-rDNA sequences were already observed within another genus within Parazoanthidae (Antipathozoanthus 43 ). These Authors found that, despite being highly variable compared to mitochondrial markers (due to presence of indels), ITS marker can leave some uncertainties regarding the specific status of different morphotypes. ITS markers are therefore unable to fully resolve the phylogenetic relationships of this group 43 , which coincides with what was already described for other groups of cnidarians (for instance Eunicella 9,44 ). Hypervariable ITS sequences can lead to alignment challenges (difference in both length and nucleotide identity), and intragenomic variation, that are widely acknowledged as being the primary obstacles to successfully using these sequences in phylogenetic inference in Zoanthidea 45 . These difficulties stress the need to increase the number of analyzed individuals to detect even small differences that can unlock the relationships among species, particularly in Zoanthidae phylogenetic studies 46 .
Despite these molecular issues, from the phylogenetic point of view, our results are in accordance with the latest accepted phylogeny of Parazoanthidae 24,25 suggesting paraphyly of the genus Parazoanthus. As observed previously by 46 , the "Slender" P. axinellae is genetically similar to P. anguicomus and, to a lesser extent, to the Atlantic-Caribbean counterparts such as Parazoanthus capensis 47 . These are Atlantic species which are able to colonize sponges just like P. axinellae, but do not depend on such associations to survive 46 . The "Stocky" morphotype showed a genetic similarity with the shallow Pacific water P. elongatus, P. juanfernandezii 24,[46][47][48] and with the deep Atlantic P. aliceae sp. n. 24 . All these species are mainly found on rocks rather than on sponges, similarly to the "Stocky" morphotype, and conversely to the "Slender" morphotype. However, despite consistency between our molecular results with those by other Authors 24,43 , genetic data should be taken with caution since discrepancies between genetic and morphological data were already documented. From a morphological point of view, P. elongatus and P. aliceae sp. n. are quite different from the P. axinellae species complex, and do not seem to share similar morphological features, which contradicts the low COI and ITS gene divergence we observed. Specifically, Ocaña and Brito 16 recently evidenced that Parazoanthus elongatus from Chile, should be placed into a different genus based on morphological data (e.g. the absence of special spirulae, the scarce presence of mineral particles in the ectoderm, and the large spirulae found in tentacles). This confirms once more that uncovering the origin of the "Stocky" morph in the western Mediterranean is a challenge. Future research using an integrative taxonomy approach 13 could complement the results of the present study regarding the status of the Mediterranean Parazoanthus genus. In fact taxonomy using only molecular tools should not be a replacement of classical taxonomy, but findings derived from one approach should be used as a driver for further investigation for the other. A multidisciplinary approach should include (1) reproductive biology of both morphotypes in sympatry, (2) field experiments to evaluate substrate specificity 26 , (3) studies of chemodiversity 49 , (4) naturalistic studies 16 and (5) more variable molecular markers (e.g. RAD sequencing 50 ; Terzin et al. unpublished) and/or sequencing of the complete mitochondrial genome 51 .
The results of this work stress the importance of studies on species delimitation and connectivity for the implementation of effective plans for the conservation of coralligenous species. In fact, other studies have observed Scientific RepoRtS | (2020) 10:13938 | https://doi.org/10.1038/s41598-020-70770-z www.nature.com/scientificreports/ that in closely related species, with similar biological features, differences in the phylogeographic patterns can occur 10,52 . Finally, our results highlight occurence of hidden diversity within the flag-genus Parazoanthus, even in well-studied geographical areas, calling for a carefull taxonomic reevaluation of other key-species of the fragile endemic Mediterranean coralligenous ecosystem 53 .

Methods
The original observation of the occurrence of two clearly distinguishable morphotypes of Parazoanthus was done by SCUBA diving in Portofino where, during the same dive, the "Slender" morphotype was observed in the deeper and in the shallower layers of a cliff, while the "Stocky" morphotype was found at an intermediate depth.
Based on these preliminary results we developed a sampling design at a Mediterranean scale identifying sites where the occurrence on Parazoanthus was recorded (Fig. 1c). We identified 11 locations, and in six of them both morphotypes occurred in sympatry (Table 1). In each location up to 30 polyps per morphotype of Parazoanthus were collected by SCUBA diving, keeping a minimum distance of 2 m between polyps to avoid sampling of clones. All samples were immediately fixed in 80% ethanol and refrigerated at 4ºC. DNA was extracted from single polyps using EuroClone EuroGold tissue DNA mini kit. Mitochondrial COI fragment was amplified using species-specific designed primers COIpaxFwd (sequence 5′-3′: CGG TAT GAT AGG AAC AGC ), and COIpaxRev (sequence5′-3′:CGG GGT CAA AGA AGG TAG TG). A fragment of the nuclear DNA including 18S, ITS-1, 5,8Sa, ITS-2, and 28S (hereafter ITS) was amplified using zoanthid-specific primers described in 54 . PCR was performed in a final volume of 25 µl per sample and included: 5 µl GoTaq Flexi Buffer 1x (Promega), 4 µl MgCl 25 mM, 0.5 µl dNTPs 10 mM, 0.5 µl of each primer (10 mM), and one unit of GO Taq G2 Flexi DNA polymerase (Promega) and filled with nuclease free water to reach the volume. Amplifications were conducted on a GeneAmp 2,700 thermal cycler (Applied Biosystems) under the following conditions: a hold at 94° C for 3′ followed by 30 cycles of denaturation at 94 °C for 45″, annealing at a primer specific temperature (59 °C for COI and 50° C for ITS) for 1′ and extension at 72 °C for 2′, finishing with a final extension at 72° for 7′.
PCR products were checked in 1.5% agarose gel stained with Gelred (BIOTIUM) 1% after a 20′ 120 V electrophoresis. They were then sent to Macrogen Europe Inc. for purification and sequencing.
Field and experimental protocols were approved by the University of Bologna, Italy and were performed in accordance with its relevant guidelines and regulations. No permit was required for the collection of the species.

Data analysis
Sequence quality check and alignment. Each sequence was checked for quality in MEGA v.6 and good quality sequences were aligned using MAFFT (Multiple Alignment using Fast Fourier Transform) 55 through phyloch R package 56 (see below for more detail).
Genetic diversity and structure among Parazoanthus axinellae samples. The number of haplotypes (H) and haplotype and nucleotide diversity (Hd and π respectively), were calculated for each sample and marker in DNAsp 57 . Haplotype networks for COI and ITS markers were built using pegas R package 35 in Rstudio (Version 1.1.453) 58 . Haplotypes were reconstructed with haploNet function using an infinite site model based on Hamming distances of DNA sequences. As there was an overlap between genetically approximate haplotypes on the Fig. 2 (in particular for the ITS gene), alluvial diagrams were created using an online platform RAWGraphs 59 to better clarify Morphotype/Haplotype/Sampling locations correlations, and Inkscape (version 0.92.4) was used to integrate the plots and finalize graph compilation.
A Minimum Spanning Tree was computed as in 60 to infer links between the most similar haplotypes based on a previously computed distance matrix, and to visualize the number of mutations between them. This was done by performing a Multidimensional Scaling (MDS) analysis on Hamming distances computed between haplotypes, using the show.mutation = T option to show the number of mismatches between linked haplotypes. The results were presented as a two-dimensional MDS plot.
Genetic differentiation between samples (with morphs from the same locality treated as different populations) was estimated using Φ statistics (Φ ST based on haplotype frequencies and molecular divergence) and its significance determined using a permutation test (10,000 permutations) for each marker in Arlequin v. 3.5 61 . Significance values were corrected for multiple comparisons following FDR correction method 62 .
To study the demographic history of the two morphotypes Fu and Li 33 and Tajima 63 ′s neutrality tests were performed in DNAsp for each locus and morphotype across the study area. Subsequently, the demographic history was also assessed by performing a mismatch distribution analysis, in which the frequencies of pairwise nucleotide differences between individuals were compared with the expected values under a sudden expansion model 34 using Arlequin. The best fit was tested by evaluating both the sum of squared deviation (SSD) and the Harpending's raggedness (HRI) indexes with a total of 1,000 permutations.
Phylogeny and differentiation within Family Parazoanthidae. We retrieved COI and ITS sequences belonging to the family Parazoanthidae from Genbank. Full details on the species, accession numbers, approximate geographic origin, and original reference are given in Supplementary Table S7 for each marker.
We decided to perform all downstream analyses at the taxonomic level of family due to large uncertainties regarding the phylogeny and systematics within Parazoanthidae at lower taxonomic categories, which led to significant and recent taxonomic modifications for these cnidarians (e.g. 25,64 55 through phyloch R package 56 . P. axinellae + Genbank COI sequences were aligned with default parameters, whereas L-INS-i algorithm with 1,000 iterations (iterative refinement method incorporating local pairwise alignment information) was used to align P. axinellae + Genbank ITS sequences due to a presence of gaps between individuals showing high sequence divergence. Gblocks Version 0.91b 65 was then ran from R to exclude gaps between highly divergent sequences, with parameters finally set to: b1 = 0.8 (the minimum number of sequences for a conserved position), b2 = 0.9, (the minimum number of sequences for a flank position), b4 = 2, (the minimum length of a block, default 2) and b5 = "h", to remove gaps present in > 50% of individuals.
Likelihood ratio test was performed with modelTest function in phangorn 66,67 to decide on the model of nucleotide evolution that best fits COI and ITS markers, and the best nucleotide substitution model was identified based on minimum AICc values for each marker separately.
Then, Bayesian inference (BI) in MrBayes v. 3.2.6 68 and Maximum Likelihood (ML) trees with UPGMA algorithm and clustering dendrograms were used to resolve phylogenetic situation of P. axinellae species complex within Parazoanthidae family based on COI and ITS gene markers. Distance matrices were calculated using the suggested F81 nucleotide substitution model and Maximum Likelihood (ML) trees with 1,000 bootstrap iterations were estimated from the obtained distance matrices using the UPGMA algorithm. This was done in phangorn 66,67 by computing the likelihood of a given tree with the function pml(), and with the function optim. pml(), which was used to optimize tree topology and branch length for F81 model of nucleotide evolution. Tree plotting was done using plotBS function, with haplotypes presented on the tree nodes in roman numerals, and bootstrap support values (based on 1,000 iterations) shown on tree branches with a cutoff value of p = 50 (%). To further clarify "haplotype-sampling location" associations, clustering dendrograms were reconstructed for each gene marker in RAWGraphs 59 and combined with ML trees in Inkscape (version 0.92.4). Hierarchically, the clustering dendrograms are composed out of 2 levels for sequences downloaded from GenBank (Species code, Number of sequences) and 3 levels for individuals sampled in our study (Morphotype Color, Sampling Location, Number of sequences). Sequences originating from our samples were colored in yellow and orange for "Slender" and "Stocky" morphotypes, respectively, while all sequences downloaded from GenBank are in grey.

Data availability
The dataset and the R codes supporting the conclusions of this article are available in fasta format as additional files in the Supplementary information.