Morphological and molecular evidence of cryptic speciation in sympatric colour morphotypes of Mycale (Carmia) cecilia (Porifera: Poecilosclerida) from the Mexican Pacific

Collection of material

Eight individuals of M. (C.) cecilia were collected in the morning on the same day (26 July 2019) from the Urias Estuary in Mazatlán, Sinaloa, Mexico (23°11’07.3”N, 106°25’18.4’’W). The sponges were found side by side, between one and two meters deep in an artificial jetty. We chose eight individuals: four with red body colour and red to reddish-orange specks and four with green body colour and green to greenish-yellow specks (Fig. 1). Hereafter, these are referred to as red (MCR) and green (MCG) morphotypes, respectively. The samples were transported to the laboratory in an aerated seawater aquarium. For morphological analysis and DNA extraction, a small fragment of each individual was fixed in 70% and 96% ethanol, respectively. For transcriptome analysis, a small portion of 50-100 µg of each fresh individual was dissected using a stereoscope to avoid macrosymbionts; the tissue was immediately homogenized in 1 mL of TRIzol^TM Reagent (Thermo Scientific; Waltham, MA) and stored at -20°C until RNA isolation. All specimens were deposited in the Mexican Pacific Sponge Collection held by the Institute of Marine Sciences and Limnology, National Autonomous University of Mexico (red morphotype, LEB-ICML-UNAM-3181 MCR1 to MCR4 and green morphotype, LEB-ICML-UNAM-3182 MCG5 to MCG8).

Transcriptome analysis

Three individuals of each morphotype were used for RNA isolation for standard RNA-seq, library preparation and sequencing. Raw reads of three individuals of each morphotype were pooled for transcriptome assembly and biological contamination detection based on Porifera best-hit. These procedures were performed and described in Castillo-Páez et al. (2021)Castillo-Páez A.Y., Llera-Herrera R., Cruz-Barraza J.A. 2021. De novo transcriptome assembly for two colour types of the marine sponge Mycale (Carmia) cecilia. Mol. Biol. Rep. 48: 1-4. https://doi.org/10.1007/s11033-021-06274-4. The raw reads of the transcriptomes for the pre-competent larvae, post-larvae and adults of Mycale (Carmia) phyllophila (NCBI SRA BioProject: PRJNA269144; accessed on 10 August 2020) were assembled and used as comparative material. The Porifera open reading frames (ORFs) of each transcriptome of M. (C.) cecilia and the M. (C.) phyllophila transcripts were used to identify orthologous by bidirectional reciprocal BLAST hit to obtain pairs of sequences with the best hit (E-value of 1e^-5). Histograms were used to present the numbers of ORFs and the percentage of amino acid identity for each transcriptome comparison (i. M. (C.) cecilia green vs. M. (C.) cecilia red; ii. M. (C.) cecilia green vs. M. (C.) phyllophila; and iii. M. (C.) cecilia red vs. M. (C.) phyllophila).

To confirm speciation trends between the two morphotypes based on phylogenetic inference of single-copy orthologous genes, we retrieved RNA-seq from the NCBI SRA BioProject (accessed on 4 August 2023) for Mycale grandis (PRJNA317402), Mycale laevis (PRJNA612016), Mycale acerata (PRJNA612010) and Mycale tridens (PRJNA612015) and Clatrhia prolifera (PRJNA901949) as an outgroup. The raw reads were transformed to paired-end .fastq files and individually assembled with Trinity software v2.15.1 (Haas et al. 2013Haas B.J., Papanicolaou A., Yassour M., et al. 2013. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat. Protoc. 8: 1494-1512. https://doi.org/10.1038/nprot.2013.084), using the -trinotate option without considering the strand-specific orientation of reads. The obtained transcriptome assemblies were assessed for basic metrics with the Trinitystats.pl script. The conceptual translations of proteins were processed with Transdecoder v. 5.7.0 included in Trinity software, retaining the HMMER and blastp hits against pFAM and UniProt-SwissProt, respectively. Only the proteins annotated as “complete” by the Transdecoder tool were retrieved for further analysis. We used the OrthoFinder tool (Emms and Kelly 2019Emms D.M., Kelly S. 2019. OrthoFinder: Phylogenetic orthology inference for comparative genomics. Genome Biol. 20: 1-14. https://doi.org/10.1186/s13059-019-1832-y) to assess the potential orthogroups among all the included species (FastME option) or a selection of those more closely related to M. (C.) cecilia morphotypes. A gene tree was obtained using the STAG approach (Emms and Kelly 2017Emms D.M., Kelly S. 2017. STRIDE: Species Tree Root Inference from Gene Duplication Events. Mol. Biol. Evol. 34: 3267-3278. https://doi.org/10.1093/molbev/msx259, Emms and Kelly 2018Emms D.M., Kelly S. 2018. STAG: Species Tree Inference from All Genes. BioRxiv. 267914. https://doi.org/10.1101/267914), in which the most likely phylogenomic relationship’s tree is obtained by collapsing individual trees obtained by all the single-copy orthologous genes. A detailed description of the retrieved RNA-seq NCBI SRA accessions and assembly metrics and the scripts and commands used to create the protein collection for each species are presented in Supplementary Material Table S1 and Table S2.

DNA extraction and ITS1 amplification and sequencing

Total DNA was extracted for preserved tissues (70% ethanol) using the Quick-DNA^TM Fecal/Soil Microbe Kits (Zymo Research, Orange, CA, USA) according to the manufacturer’s protocol. The Internal Transcribed Spacer 1 (ITS1) locus was amplified with primers designed by Adlard and Lester (1995)Adlard R., Lester R.J. 1995. Development of a diagnostic test for Mikrocytos roughleyi, the aetiological agent of Australian winter mortality of the commercial rock oyster, Saccostrea commercialis (Iredale and Roughley). J. Fish Dis. 18: 609-614. https://doi.org/10.1111/j.1365-2761.1995.tb00365.x RA2 5’-GTCCCTGCCCTTTGTACACA-3’ and ITS2.2 5’-CCTGGTTAGTTTCTTTTCCTCCGC-3’. PCR reactions were conducted in a final volume of 15 µL containing ~20 ng of template DNA, 1x NH4 buffer (Bioline), 2 mM MgCl₂, 0.2 mM dNTP, 0.4 µM of each primer, and 0.5 U of MyTaq DNA polymerase (Bioline). Amplification conditions were: 95°C for 5 min; then 35 cycles of 94°C for 30 s, 49°C for 30 s and 68°C for 1 min; and a final elongation step of 72°C for 10 min. Because of unspecific bands, PCR products were separated by electrophoresis in 1.5% agarose gels. DNA bands were recovered with the Wizard® SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA) and were processed by capillary sequencing through a third-party service at Macrogene Inc. (Korea).

Phylogenetic relationships and species delimitation based on traditional loci

Genetic loci traditionally used in phylogenetic studies were recovered from each transcript set [nuclear 18S and 28S rRNA and mitochondrial 16S rRNA and cytochrome c oxidase subunit 1 (COX1) genes] or Sanger sequences for ITS1 (MCR and MCG). BLASTx searches were performed using Mycale spp. sequences downloaded from GenBank as a reference for each genetic marker. Genetic distances were calculated with the uncorrected p-distance index for each marker, assuming uniform rates among sites and conducting pairwise deletion of gaps and missing data. The best-fitting DNA substitution model was inferred from the complete alignments using an initial neighbour-joining tree and the maximum likelihood statistical method using all sites, and a strong branch swap was used for gaps and missing data (5 threads). Phylogenetic reconstructions were made from each complete alignment with the maximum likelihood method with 300 bootstrap replicates using the best-fitting substitution model defined previously, and the tree inference nearest-neighbour-interchange with a very strong branch swap filter. These analyses were done in MEGA X (Kumar et al. 2018Kumar S., Stecher G., Li M., et al. 2018. MEGA X: Molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35: 1547-1549. https://doi.org/10.1093/molbev/msy096). Assemble Species by Automatic Partitioning (ASAP) software on the web server was run for species delimitation, applying the Kimura (K80 ts/tv 2.0) substitution model to compute the distances (Puillandre et al. 2021Puillandre N., Brouillet S., Achaz G. 2021. ASAP: assemble species by automatic partitioning. Mol. Ecol. Resour. 21: 609-620. https://doi.org/10.1111/1755-0998.13281).

Morphology and spicule size

Spicule preparations were made following the techniques described by Carballo and Cruz-Barraza (2008)Carballo J.L., Cruz-Barraza J.A. 2008. First record of Axinyssa Lendenfeld, 1897 (Demospongiae, Halichondrida) from the East Pacific Ocean, with the description of Axinyssa isabela sp. nov. Zootaxa. 1784: 58-68. https://doi.org/10.11646/zootaxa.1784.1.4. The different types of spicules were visualized under an Olympus® CX21 light microscope with an ocular micrometer rule. Twenty-five mycalostyles, 25 sigmas and 50 palmate anisochelae were randomly measured for each individual. Length, shaft width and head width were measured for each mycalostyle. Total and head lengths were taken for palmate anisochelae. For sigmas, only the length was recorded. A histogram was used to visualize the frequency and density of the spicule sizes. Statistical analyses were conducted to test significant differences in the spicular mean between individuals. After evaluating the homoscedasticity with the Bartlett test or Fligner-Killeen test, transformations (natural logarithm or square root) were implemented for normalization. We used analysis of variance with a significance level of 0.05. An a posteriori Tukey HSD test was conducted to identify the most different individuals. The multidimensional distribution of the samples based on morphometric data was explored using a principal component analysis (PCA). All analyses were performed in RStudio (version 1.2.5001, RStudio, Inc.).

Phylogenetic relationships and species delimitation inferred from traditional loci

Sequences from ITS1 were successfully obtained for all specimens examined. ITS1 sequences from the red individuals were 332 nt long, and those from the green individuals were 333 or 334 nt long. The number of differences (0 to 1) and genetic distances (0% to 0.3%) within each colour (red or green) morphotype group was extremely low when they occurred. However, when comparisons were made between individuals of different morphotypes of M. (C.) cecilia, the substitutions and differences were more evident (MCR vs. MCG in Table 1). Notably, the differences between M. (C.) cecilia morphotypes (MCR vs. MCG) were greater than those between MCG and M. (C.) phyllophila (Table 1). This locus was highly variable compared with the other analysed genes. The best-fitting substitution model was T92 + I. Individuals of M. (C.) cecilia with the same colour were grouped in a strongly supported monophyletic clade (bootstrap value 100%, Fig. 2). M. (C.) phyllophila clustered into the MCG clade, suggesting that this is a sister species to MCG. Four best partitions were found by ASAP for the ITS region sequences. The first partition was the best supported (p-value <0.05) and indicated that there were six subsets. In one of them there were only the samples of the green morphotype and in another only the samples of the red morphotype. M. phyllophila did not group with the green morphotype.

The partial consensus sequences for M. (C.) cecilia from the 28S rRNA D1-D2 domain had a length of 905 nt. The best-fitting model of substitution was TN93+G. The sequences were identical between MCG and M. (C.) phyllophila (Table 1). In addition, a low value was observed in pairwise comparisons between the morphotypes (MCR and MCG) of M. (C.) cecilia, and this value was the lowest when compared with another Mycale spp. (except M. (C.) phyllophila). The phylogenetic construction was well resolved. Both the morphotypes and M. phyllophyla were clustered, although MCG and M. phyllophyla may be siblings because they clustered in a highly supported branch (100% bootstrap, Fig. 3). Nine partitions were detected by ASAP analyses, but they were not informative (p-value > 0.1).

The other domain (D3-D5) was not informative, as no nucleotide differences were found between morphotypes or in other species, such as M. mirabilis vs. M. setosa, M. fistulifera vs. M. microsigmatosa and M. (C.) cecilia from Hawaii (KC961688.1) vs. M. adhaerens. Phylogenetic reconstruction was performed with 420 nt. The best-fitting model was K2+G+I. None of the clades was well supported and all bootstrap values were less than 75% (even 0). However, the two morphotypes (MCR and MCG) clustered in a moderately supported clade (65% bootstrap; Fig. S1). Partitions detected for species delimitation were not useful. The two morphotypes were clustered even with other Mycale species.

Partial consensus sequences of the 18S rRNA gene (1728 nt) for each morphotype were retrieved from the assembled transcriptome. The uncorrected p-distances ranged between 0% and 5.6% (Table 1). The lowest value was found between MCG and M. (C.) phyllophila, where sequences were identical. The genetic distance between MCR and MCG was slightly higher than that between other species, i.e. MCG vs. M. microsigmatosa (0.17) or M. sanguinea (0.23). We found 14 substitutions between M. (C.) cecilia morphotypes and none between M. (C.) phyllophila and MCG. The best-fitting model of substitution was K2+G+I. The phylogenetic reconstruction had low support values but was useful to show that our morphotypes were not monophyletic. MCG vs. M (C.) phyllophila clustered with a high bootstrap value (Fig. S2). The species delimitation analysis was not congruent; the two morphotypes were grouped even with other Mycale species.

We analysed two mitochondrial genes. COX1 was more informative than the 16S rRNA gene. For COX1, there were seven substitutions at 1261 nt between MCR and MCG. The sequences were identical between MCG and M. (C.) phyllophila (Table 1). The phylogenetic tree was made with HKY+G as the best-fitting substitution model. Both M. (C.) cecilia and M. (C.) phyllophila morphotypes clustered in a single well-supported clade (100%, Fig. 4). The ten best partitions were found by ASAP, but they were not informative (>0.1).

For 16S, the partial consensus sequence was 586 nt long for both morphotypes. The best-fitting model of substitution was HKY+G. The value of uncorrected p-distance ranged between 0% and 15%. The lowest values were between MCG or MCR vs. other Mycale species including M. (C.) phyllophila and M. (C.) cecilia (Table 1). The number of divergences ranged from 0 to 90 (Table 1). There was no difference between MCR and one of the M. (C.) cecilia sequences from Hawaii. Neither the other M. (C.) cecilia sequence from Hawaii nor M. (C.) phyllophila differed from MCG. There was just one substitution (p-distance = 0.17%) between M. (C.) cecilia colour morphotypes. The phylogenetic relationship was not well resolved. Both morphotypes were in a polytomy clade, perhaps because of extremely low divergence values caused by the low substitution rate that this gene shows (see Supplemental Material Fig. S3). The ten best partitions were found by ASAP, but they were not informative (>0.1).

Similarity between red and green morphotype transcriptomes

The sequencing and transcriptome data (NCBI BioProject PRJNA678333) of morphotypes of M. (C.) cecilia are described in detail by Castillo-Páez et al. (2021)Castillo-Páez A.Y., Llera-Herrera R., Cruz-Barraza J.A. 2021. De novo transcriptome assembly for two colour types of the marine sponge Mycale (Carmia) cecilia. Mol. Biol. Rep. 48: 1-4. https://doi.org/10.1007/s11033-021-06274-4. In summary, the RNA-seq of six samples yielded ~129 million raw reads with a length of 150 nt. Each morphotype was assembled separately because of the global assembly’s low overlap. Transcript totals of 461266 for MCG and 342237 for MCR were obtained. Of these, 153498 (MCG) and 122979 (MCR) were long ORFs. Only 48023 (MCG) and 37443 (MCR) ORFs had hits with Porifera. The raw RNA-seq reads for the three life cycle stages of M. (C.) phyllophila were assembled de novo, and 384743 transcripts were obtained.

Potentially orthologous genes (2977) were recovered by bidirectional best-hit with BLASTP after contrasting the sequences assembled from each morphotype of M. (C.) cecilia with the M. (C.) phyllophila transcripts and the Porifera ORFs (hereafter referred to as orthologous genes). The percentage of amino acid identity of the 2977 orthologous genes between MCR and MCG was lower than that between MCG and M. (C.) phyllophila (Fig. 5). MCR vs. MCG had an average identity of 95.14%: 102 genes were identical, while 245 genes were 99% identical between each data set. This is much lower than the figure for MCG vs. M. (C.) phyllophila, which showed an identity of 98.59%: 661 genes showed homology at 100% and 1210 at 99%. The percentage of amino acid identity values of MCR vs. M. (C.) phyllophila was quite similar to that of the M. (C.) cecilia morphotypes, with an average identity of 95.02%: 106 were identical and 252 were 99%.

Phylogenetic relationships inferred from transcriptomes

The number of transcripts after assembly for Mycale species and outgroup ranged from 85915 for Clathria prolifera to 299581 for M. tridens (see Supplementary Material Table S1). The number of complete proteins ranged from 17184 for Mycale acerata to 54915 for M. laevis. Only 12 single-copy orthologous genes were present in all species analysed. The phylogenetic reconstruction showed the same patterns as the phylogenetic reconstructions of traditional loci: M. (C.) phyllophila appears to be a sister to MCG (Fig. 6). The same pattern was observed when the number of complete genes (75 and 106) was increased by excluding samples that did not have them (see Supplementary Material Fig. S4).

External body shape and colour

Individuals displayed a similar external body shape in life, differing only in coloration. The organisms were mottled patchwork red to reddish-orange or green to greenish-yellow and covered by a translucent dermal membrane (Fig. 1). They were encrusted about 0.5 cm thick and fine-grained; the consistency was very soft and fragile with a smooth surface.

Spicule morphometry

Regarding the specimens’ spicules of both morphotypes (MCR and MCG), the length of mycalostyles ranged from 205 to 300 µm. The averages varied from 236.7 (individual MCR2) to 261.8 µm (individual MCG5; Table 2). In most individuals, there was a normal distribution in the length-frequency of mycalostyles, but it was most robust when individuals with the same colour were grouped (see Supplementary Material Fig. S5). There were strong significant differences (F = 6.73, p = 3.8e^-07) between the eight individuals analysed and when pooled together with individuals with the same body colour (red vs. green, F = 24.27, p = 1.76e^-06). An a posteriori test revealed that the MCG5 individual was the most different, showing significant differences from three red organisms and from one of the same body colour (green). Furthermore, the shaft width of mycalostyles varied between individuals. It fluctuated from 2.5 µm in both morphotypes to 7.5 µm in the green (MCG5) and 10 µm in the red (MCR2) individuals. The mean ranged from 4.52 (MCG6) to 5.94 µm (MCR2; Table 2). There were significant differences (F = 5.61, p = 6.6e^-06) in the shaft width of mycalostyles. Individual MCR2 was the most different. It varied with two individuals of the same body colour and with three of the green body colour. Individual MCG6 was different from one individual of each body colour. The other measure we analysed was the width of the head of mycalostyles. It ranged from 2.5 to 7 µm in red individuals and from 3 to 8.75 µm in green individuals. The means ranged from 4.40 (MCG7) to 5.97 µm (MCR2). Significant differences were found between individuals (F = 6.84, p = 2.9e^-07). The MCR2 individual was the most divergent in width of the head of mycalostyles, and it was different for all individuals except for one individual of the different body colour (green).

The length of sigmas ranged from 25 µm in the red and 27.5 µm in the green to 42.5 µm in the organisms of both colours. The average varied from 31.20 (MCR4) to 37.45 µm (MCG7; Table 2). Normal distribution in the frequency was found in most individuals (see Supplementary Material Fig. S6). Significant differences were found between individuals (F = 8.72, p = 2.6e-⁰⁹). The individual MCG7 was the most different in sigma length, which varied from all the individuals with red body colour and one with green body colour. We also found highly significant differences (F = 34.51, p = 1.76e-⁰⁸) when we collapsed all sigma lengths by colour into a data set and compared them.

The total length range of palmate anisochelae was 13.75 to 25 µm in the red morphotype and 12.5 to 25 µm in the green morphotype (Table 2). The average ranged from 17.72 (MCG7) to 19.51 µm (MCG5). All green individuals had less frequency in length between 16.5 and 19.75 µm (Fig. 7). This possibly corresponds to two categories: Category I, 12.5 to 16.25 µm; and Category II, 20 to 25 µm. The density curve for grouped green individuals (Fig. 7, MCG-All) was bimodal with a clear gap, but this was not found in grouped red individuals, which did not have a clear gap (MCR-All). Homoscedasticity was solely found when individuals were grouped by colour. Significant differences were found between red individuals (F = 5.75, p = 8.47e^-04), specifically between MCR1-MCR2 and MCR1-MCR3. Green individuals were also significantly different (F = 5.08, 2.04e^-03): MCG5-MCR6, MCG5-MCR7 and MCG7-MCR8. The head length of palmate anisochelae ranged from 7.5 to 12.5 µm in the red individuals and from 5 to 13.75 µm in the green individuals (Table 2). The average of the red morphotype was 10 µm, but the green ones ranged between 8.94 and 10.2 µm. The variation pattern in head length matched (or mirrored) the pattern of total length reported above. Only green individuals had less frequency in head length between 8.25 and 9.75 µm (Fig. 8), suggesting that there are two categories: Category I, 6.25 to 8 µm; and Category II, 10 to 13.75 µm. The density curve for grouped green individuals (Fig. 8, MCG-All) was bimodal. Homoscedasticity was found when individuals were grouped by colour. No significant differences were found between red (F = 0.99, p = 0.39) individuals and differences (F = 4.51, p = 4.45e^-03) were detected between only MCG8-MCG6 and MCG8-MCG7 green individuals.

No colour morphotype clusters were observed in PCA analyses. This suggests that the measured variables did not help differentiate between morphotypes. The length of the different types of spicules was correlated. However, these values were negatively correlated with the head or shaft width of mycalostyles, suggesting that it is sufficient to measure only the head or shaft width and the length of mycalostyles (Fig. 9). The size of the palmate anisochelae was the variable that contributed least to the morphometric ordination because it did not have a normal distribution, as mentioned above. In conclusion, no clusters by colour morphotype were found, and the variable that contributed the least in the PCA analysis (i.e. anisochelae length) is one of those that shows an evident difference in the histogram analysis.

Congruence in nuclear and mitochondrial traditional loci

Molecular markers have been established as powerful tools for revealing the biodiversity of “hidden” marine sponges, such as cryptic or complex species that cannot be resolved with only morphological characters (Uriz and Turon 2012Uriz M.J., Turon X. 2012. Sponge ecology in the molecular era. Adv. Mar. Biol. 61: 345-410. https://doi.org/10.1016/B978-0-12-387787-1.00006-4, Shaffer et al. 2019Shaffer M.R., Davy SK., Bell J.J. 2019. Hidden diversity in the genus Tethya: comparing molecular and morphological techniques for species identification. Heredity. 122: 354-369. https://doi.org/10.1038/s41437-018-0134-6). In some sponge studies, nuclear gene fragments have been more informative than mitochondrial ones (Bucklin et al. 2011Bucklin A., Steinke D., Blanco-Bercial L. 2011. DNA Barcoding of Marine Metazoa. Ann. Rev. Mar. Sci. 3: 471-508. https://doi.org/10.1146/annurev-marine-120308-080950, Uriz and Turon 2012Uriz M.J., Turon X. 2012. Sponge ecology in the molecular era. Adv. Mar. Biol. 61: 345-410. https://doi.org/10.1016/B978-0-12-387787-1.00006-4). Partial ITS1 nuclear locus showed low genetic variation levels within M. (C.) cecilia morphotypes and high variation between them. Surprisingly, the sequence of M. (C.) phyllophila organisms collected from China was more closely related to the MCG. The sequences of M. (C.) phyllophila were quite similar to the MCG sequences, with genetic distances and the number of substitutions being the lowest in the pairwise comparisons and clustered in a highly supported clade. However, in the species delimitation analysis, M. (C.) phyllophila and MCG appeared to belong to different species because they were included in different subsets. The same occurred with the two morphotypes of M. (C.) cecilia where the MCR sequences were grouped into a single subset and the MCG sequences into another subset. These results indicate that each morphotype belongs to a different species. Park et al. (2007)Park M-H.H., Sim C-J.J., Baek J., Min G-S.S. 2007. Identification of genes suitable for DNA barcoding of morphologically indistinguishable Korean Halichondriidae sponges. Mol. Cells. 23: 220-227. suggest that the ITS region (ITS1 and ITS2) could be more appropriate for molecular barcoding in sponges than mitochondrial genes. They indicated that the ITS1 locus exhibited the highest levels of intra- and interspecific variation in uncorrected p-distances compared with ITS2 and 5.8S. Intraspecific (1.52%) and interspecific (17.28%) maximum values of genetic divergence for Halichondriidae were an order of magnitude higher (Park et al. 2007Park M-H.H., Sim C-J.J., Baek J., Min G-S.S. 2007. Identification of genes suitable for DNA barcoding of morphologically indistinguishable Korean Halichondriidae sponges. Mol. Cells. 23: 220-227.). The same pattern was observed for Suberites diversicolor. The average p-distance within the lineages (lineage A 0.25% or lineage B 0.44%) was lower than that between lineages (7.26%) (Becking et al. 2013Becking L.E., Erpenbeck D., Peijnenburg K.T.C.A., de Voogd N.J. 2013. Phylogeography of the Sponge Suberites diversicolor in Indonesia: Insights into the Evolution of Marine Lake Populations. PLoS ONE 8: e75996. https://doi.org/10.1371/journal.pone.0075996), and similar values were found in the morphotypes of M. (C.) cecilia (0.3% intra- and 8.78% inter-lineage). ITS1 could be useful for taxonomic identification and genetic divergence for the other body colour morphotypes registered for M. (C.) cecilia and could be extended to other species of Mycale.

We obtained two domains of the 28S rRNA of the transcriptome. The D1-D2 domain was useful and matched with the phylogenetic reconstructions of another nuclear marker (ITS1). The genetic distance between morphotypes was low but not the lowest. MCG and M. (C.) phyllophila clustered in a strongly supported clade and were included in the same subset, suggesting that M. (C.) phyllophila could be a sister or the same species to the green morphotype. This extremely low divergence was recorded in another sponge group where only one substitution between Pachymatisma johnstonia and Pachymatisma normani was found for this domain (Cárdenas et al. 2010Cárdenas P., Rapp H.T., Schander C., Tendal O.S. 2010. Molecular taxonomy and phylogeny of the Geodiidae (Porifera, Demospongiae, Astrophorida) - Combining phylogenetic and Linnaean classification. Zool. Sci. 39: 89-106. https://doi.org/10.1111/j.1463-6409.2009.00402.x). Borchiellini et al. (2000)Borchiellini C., Chombard C., Lafay B. NB-E. 2000. Molecular systematics of sponges (Porifera). Hydrobiologia. 420: 15-27. https://doi.org/10.1007/978-94-017-2184-4_2 suggest that the D1-D2 domain appears more informative and provides better resolution for internal topologies in sufficiently closely related species. The genus Hemimycale was the only phylogenetic reconstruction with monophyly and low bootstrap value (Uriz et al. 2017Uriz M.J., Garate L., Agell G. 2017. Molecular phylogenies confirm the presence of two cryptic Hemimycale species in the mediterranean and reveal the polyphyly of the genera Crella and Hemimycale (Demospongiae: Poecilosclerida). PeerJ 5:e2958. https://doi.org/10.7717/peerj.2958). Cárdenas et al. (2010)Cárdenas P., Rapp H.T., Schander C., Tendal O.S. 2010. Molecular taxonomy and phylogeny of the Geodiidae (Porifera, Demospongiae, Astrophorida) - Combining phylogenetic and Linnaean classification. Zool. Sci. 39: 89-106. https://doi.org/10.1111/j.1463-6409.2009.00402.x used D1-D2 and D3-D5 of the 28S rRNA and COX1 to evaluate the phylogeny of Geodiidae. They observed patterns similar to those found in M. (C.) cecilia. First, there was extremely low divergence in the D1-D2 domain and COX1, and the D3-D5 phylogeny was inconsistent with the D1-D2 domain and COX1. The authors suggest that the D3-D5 domain could evolve quite differently from COX1 and the D1-D2 domain. Furthermore, they suggest that the latter could be an excellent nuclear marker for barcode rather than COX1. This is exactly what happened among the morphotypes of M. (C.) Cecilia. We found D1-D2 had more substitutions than COX1, and there was phylogeny of the D3-D5 incongruence.

The COX1 and I3-M11 partitions are commonly used for species delimitation and population structure research. COX1 is also a marker of slow evolution in sponges and sometimes fails to resolve phylogenetic relationships and speciation (Belinky et al. 2012Belinky F., Szitenberg A., Goldfarb I., et al. 2012. ALG11 - A new variable DNA marker for sponge phylogeny: Comparison of phylogenetic performances with the 18S rDNA and the COI gene. Mol. Phylogenet. Evol. 63: 702-713. https://doi.org/10.1016/j.ympev.2012.02.008; see Uriz and Turon 2012Uriz M.J., Turon X. 2012. Sponge ecology in the molecular era. Adv. Mar. Biol. 61: 345-410. https://doi.org/10.1016/B978-0-12-387787-1.00006-4). Though both morphotypes of M (C.) cecilia showed only seven substitutions in 1261 nt, and five of them were in the first 717 nt, the phylogenetic relationship (Fig. 4) between morphotypes was well resolved. In 99% of the phylogenetic reconstructions, MCG was clustered with M. (C.) phyllophila, rather than MCR. The genetic divergence between the morphotypes of M. (C.) cecilia was low (0.5) when compared with other Mycale species, except for the green morphotype and M. (C.) phyllophila, which was zero. A lower value than this was obtained between lineage A and B for S. diversicolor, where the divergence was 0.38 (Becking et al. 2013Becking L.E., Erpenbeck D., Peijnenburg K.T.C.A., de Voogd N.J. 2013. Phylogeography of the Sponge Suberites diversicolor in Indonesia: Insights into the Evolution of Marine Lake Populations. PLoS ONE 8: e75996. https://doi.org/10.1371/journal.pone.0075996).

The nuclear 18S rRNA gene and the mitochondrial 16S rRNA were non-informative in our phylogenetic reconstructions and speciation delimitation. However, the trees have similar patterns to the other phylogenetic reconstructions. All subsets of different partitions overlapped well-defined species, indicating the low resolution for species delimitation. Borchiellini et al. (2000)Borchiellini C., Chombard C., Lafay B. NB-E. 2000. Molecular systematics of sponges (Porifera). Hydrobiologia. 420: 15-27. https://doi.org/10.1007/978-94-017-2184-4_2 suggest that 18S rRNA loci are not informative for resolving relationships between genera.

Our molecular findings are opposite to those found in the congener M. (Mycale) laevis, where four morphotypes (orange or white, semi-cryptic or massive) showed no genetic differences, suggesting that they all constitute a single species (Loh et al. 2012Loh T.L., López-Legentil S., Song B., Pawlik J.R. 2012. Phenotypic variability in the Caribbean Orange Icing sponge Mycale laevis (Demospongiae: Poecilosclerida). Hydrobiologia. 687: 205-217. https://doi.org/10.1007/s10750-011-0806-1). This finding is relevant because it suggests that body colour is unrelated to genetic similarity. There could be more Mycale species with very low levels of divergence and they are probably unpredictable when evaluated by few genes. However, the molecular markers used may not have been sufficiently informative. It would be interesting to assay this species again using the ITS1, 28S rRNA D1-D2 and COX1.

Reports of cryptic sympatric speciation in sponges have increased in recent years in different groups of marine sponges. Cryptic species are widespread across Porifera (Knapp et al. 2015Knapp I.S., Forsman Z.H., Williams G.J., et al. 2015. Cryptic species obscure introduction pathway of the blue Caribbean sponge (Haliclona (Soestella) caerulea), (order: Haplosclerida) to Palmyra Atoll, Central Pacific. PeerJ. 2015. https://doi.org/10.7717/peerj.1170). Both mitochondrial and nuclear markers revealed at least three different species for the only known Scopalina lophyropoda, as in the case of our M. (C.) cecilia. Scopalina blanensis can even be found in contact with the true S. lophyropoda. These species of Scopalina also have little morphological differentiation including coloration: bright brownish orange in S. lophyropoda and salmon-orange in S. blanensis (Blanquer and Uriz 2008Blanquer A., Uriz M-J. 2008. “A posteriori” searching for phenotypic characters to describe new cryptic species of sponges revealed by molecular markers (Dictyonellidae: Scopalina). Invertebr. Syst. 22: 489-502. https://doi.org/10.1071/IS07004). Individuals of Haliclona (Soestella) caerulea species were morphologically indistinguishable, but in mitochondrial analysis, haplotypes can be up to 13.9% divergent, and even in the concatenated mitochondrial and nuclear tree they revealed two major clades (Knapp et al. 2015Knapp I.S., Forsman Z.H., Williams G.J., et al. 2015. Cryptic species obscure introduction pathway of the blue Caribbean sponge (Haliclona (Soestella) caerulea), (order: Haplosclerida) to Palmyra Atoll, Central Pacific. PeerJ. 2015. https://doi.org/10.7717/peerj.1170).

Morphotypes similarity at the transcriptome level

Our study is probably the first assessment at the transcriptomic level for species delimitation in marine sponges. Where we found the same results in individual genes, we found much more divergence than expected between the two M. (C.) cecilia morphotypes. The global assembly (all reads of M. (C.) cecilia morphotypes) had a low overlap between body colour morphotypes, suggesting a high degree of divergence between them. Further, the average amino acid identity of the orthologous genes between MCR and MCG (95.14%) was much lower than that between MCG and M. (C.) phyllophila (98.59%); thus, many of these orthologous genes were more identical between MCG and M. (C.) phyllophila. The same was true in the phylogenetic reconstruction of the complete genes of the transcriptomes with different Mycale species. The phylogenetic reconstructions of 12, 75 and 106 complete single-copy orthologous genes show the same pattern among them and with the trees of traditional genes. Namely, both M. (C.) cecilia morphotypes and M. (C.) phyllophila clustered but the green morphotype appears to be a sister to the M. (C.) phyllophila in place of the red morphotype.

The number of orthologous genes found among the three species (M. (C.) cecilia morphotypes and M. (C.) phyllophila) was quite similar to that found in recent research. In the phylogenomic inference of the Spongillidae, 2710 orthogroups were detected among 11 sponge species. These concatenated orthologous genes were enough to reveal conclusively the inter-relationships of freshwater sponges (Kenny and Itskovich 2020Kenny N.J., Itskovich V.B. 2020. Phylogenomic inference of the interrelationships of Lake Baikal sponges. Syst. Biodivers. 0: 1-9.).

The results obtained from the transcriptome comparison of MCG vs. M. (C.) phyllophila were unexpected. Firstly, the coloration of M. (C.) phyllophila is red; so, we would expect it to be more similar to MCR, but the opposite was found. Second, we did not expect a high level of similarity in the identity of the amino acid of orthologous genes between MCG and M. (C.) phyllophila (98%). We expected the similarity to have a similar value to other species comparisons (95%). This great similarity is interesting because M. (C.) phyllophila was sampled in China and MCG in Mexico. So far, it is unclear why two species that inhabit areas almost 8000 miles apart are much more similar than those that coexist.

Consistent morphological differences

Following the original description of M. (C.) cecilia, both Mazatlán morphotypes (red and green) had mycalostyles in a single size category, without apparent differences between them. These mycalostyles were slightly smaller and narrower than those of the Panama individuals but similar to those from Hawaii (see the Introduction and Table 2). Both morphotypes had sigmas in a single size category without length differences. The sigmas were also identical in size to specimens from Panama and Hawaii. Two green individuals, one for mycalostyle length and one for sigma length, were different even with individuals of the same body colour, suggesting that other spicular measurements overlapped between morphotypes. Further, we found differences between morphotypes when the analysis was done in body colour pools, i.e. all the measures without considering the individuals, only the body colour. We did not find morphotype clusters in PCA analysis, suggesting that spicule morphometry was unsuitable for delimiting morphs, even when some of them were correlated by length. The main spicule difference between morphotypes, and perhaps the most significant diagnostic key, along with coloration, was the presence/absence of two discrete categories in the total and head length of palmate anisochelae. According to de Laubenfels (1936)de Laubenfels M.W. 1936. A comparison of the shallow-water sponges near The Pacific and of The Panama Canal with those at The Caribbean end. Proc U. S. Natl. Mus. 83: 441-466. https://doi.org/10.5479/si.00963801.83-2993.441, the type material of M. (C.) cecilia was green in life and showed two size categories of palmate anisochelae, which seems similar to our green specimens. The red ones only showed a single anisochela category. de Laubenfels (1950)de Laubenfels M.W. 1950. The sponges of Kaneohe Bay, Oahu. Pac. Sci. 4: 3-36. reported Hawaii individuals with different colourations (pink and lavender or yellow and orange), but without differences in anisochela lengths. Some authors have argued that spicule size categories can be useful as a diagnostic character (Hajdu and Rützler 1998Hajdu E., Riitzler K. 1998. Sponges, genus Mycale (Poecilosclerida: Demospongiae: Porifera), from a Caribbean mangrove and comments on subgeneric classification. Proc. Biol. Soc. Wash. 111: 737-773.). We suggest then that the colour differences, along with the presence/absence of anisochela size categories, could provide sufficient morphological arguments to separate the two morphotypes into different species, but perhaps they are not the only ones.

We must consider whether other body colour morphotypes (i.e. orange, yellow and blue) recorded along the M. (C.) cecilia distribution (Tropical Eastern Pacific and Hawaii) could be different species, especially red and orange morphotypes because they have widespread distribution (Carballo and Cruz-Barraza 2010Carballo J.L., Cruz-Barraza J.A. 2010. A revision of the genus Mycale (Poecilosclerida: Mycalidae) from the Mexican Pacific Ocean. Contrib. Zool. 79: 165-191. https://doi.org/10.1163/18759866-07904003). Commonly, coloration has been considered to be without taxonomic value in sponge systematics because the colour could vary as a function of the environmental conditions or when the animals are moved from one place to another (Park et al. 2007Park M-H.H., Sim C-J.J., Baek J., Min G-S.S. 2007. Identification of genes suitable for DNA barcoding of morphologically indistinguishable Korean Halichondriidae sponges. Mol. Cells. 23: 220-227., Blanquer and Uriz 2008Blanquer A., Uriz M-J. 2008. “A posteriori” searching for phenotypic characters to describe new cryptic species of sponges revealed by molecular markers (Dictyonellidae: Scopalina). Invertebr. Syst. 22: 489-502. https://doi.org/10.1071/IS07004). Nevertheless, in some species identification, the body colour could be used as an additional character, especially when colour differences are slight but constant (Blanquer and Uriz 2008Blanquer A., Uriz M-J. 2008. “A posteriori” searching for phenotypic characters to describe new cryptic species of sponges revealed by molecular markers (Dictyonellidae: Scopalina). Invertebr. Syst. 22: 489-502. https://doi.org/10.1071/IS07004, Esteves et al. 2018Esteves E.L., de Paula T.S., Lerner C., et al. 2018. Morphological and molecular systematics of the “Monanchora arbuscula complex” (Poecilosclerida: Crambeidae), with the description of five new species and a biogeographic discussion of the genus in the Tropical Western Atlantic. Invertebr. Syst. 32: 457-503. https://doi.org/10.1071/IS16088).

According to our results, it is possible to attribute the green individuals to the M. (C.) cecilia species. These specimens shared the colouration and the presence of two anisochela categories described by de Laubenfels. In contrast, we found that the red individuals could correspond to another species of Mycale, perhaps a new species. To confirm whether the MCR is a new species, the other colour morphotypes described in the M. (C.) cecilia distribution and other species belonging to the genus Mycale need to be examined. It is possible that these (red and green) and other morphotypes described for M. (C.) cecilia in the same geographic location could be involved in a species complex because they are closely related and morphologically very similar, but the boundaries between them are unclear. We shall continue working to unravel the phenotypic plasticity bases of M. (C.) cecilia cryptic species at different levels, including genetic markers and gene expression profiles, to improve our knowledge of this species complex.