Confirmation of the southern African distribution of the marine sponge Hymeniacidon perlevis (Montagu, 1814) in the context of its global dispersal

Museum material

Stephens (1915) material, Hymeniacidon caruncula (Bowerbank, 1858), NMSZ 1921.143.1443, from False Bay and Saldanha Bay, and Leucophloeus stylifera Stephens, 1915, Syntype, NMSZ 1921.143.1443, from Saldanha Bay, were acquired on loan for comparative studies. The specimens are kept in 70% ethanol in the Department of Natural Sciences, National Museums Collection Centre in Edinburgh, United Kingdom. The holotype of Hymeniacidon sublittoralis Samaai & Gibbons, 2005, SAM-4903, is preserved in 70% ethanol and deposited at the Iziko Museums of South Africa.

South African Halichondria caruncula recorded by Stephens (1915), Leucophloeus styliferus, and H. sublittoralis were examined and compared with the South African and Namibian specimens of H. perlevis.

Sample collection

Between 2015 and 2019, field surveys were conducted during low tide at 53 intertidal habitats along the South African coast from Port Nolloth (Benguela Current system, 26°44′7.05″S; 15°5′43.38″E) to Kosi Bay near the Mozambique border (Agulhas Current system, 26°55′46.21″S; 32°52′41.23″E) (Fig. 2; Table S2). Dr. Maya Pfaff collected a sponge sample from Grosse Bucht south of Lüderitz (Namibia) in 2019 as part of the Department of Forestry, Fishries and Environment (DFFE) (South Africa) and Ministry of Fisheries and Marine Resources (Namibia) joint rocky shore monitoring program (Kreiner et al., 2019a; Kreiner et al., 2019b; Kreiner et al., 2019c), and it was identified as H. perlevis. This sample was included in this study. We were not able to sample the rocky intertidal area between Bettys Bay and Knysna (Fig. 2). Sponges were collected from the intertidal rocky shores by removing a representative piece of the animal. Observations on appearance in life, habitat type and depth were recorded in situ (Fig. 3). Colour photographs were taken in situ (Fig. 3) and in the laboratory. Upon collection, specimens were stored in 96% ethanol and processed for histological examinations according to Samaai & Gibbons (2005). Spicule dimensions are given as the mean length (range) × mean width (range) of 20 spicule measurements.

Molecular analyses

DNA was extracted from 54 representative tissue samples across the three major biogeographic provinces using the E.Z.N.A Tissue DNA kit according to the manufacturer’s protocol (Omega Bio-Tek). A fragment of the mitochondrial cytochrome c oxidase subunit I (COI) was amplified using primers LCO–1490 (5′—GGT CAA CAA ATC ATA AAG ATA TTG G—3′) and HCO–2198 (5′–TAA ACT TCA GGG TGA CCA AAA AAT CA–3′) (Folmer et al., 1994).

Polymerase chain reactions (PCR) were performed in volumes of 25 µl containing 12.5 µl Taq, 0.5 µl of each primer (10 mM), 1 µl of BSA, 5 µl of DNA template and 5.5 µl H₂O. The cycling profile included an initial denaturation step (3 min at 94 °C), 40 cycles of denaturation (30 s at 94 °C), annealing (20 s at 45 °C) and extension (1 min at 72 °C), and a final extension step (10 min at 72 °C). The amplified DNA was purified with a PCR Clean-Up Kit according to the manufacturer’s protocol. The final DNA product was sequenced in both directions on an Applied Biosystems 3730xl DNA Analyzer (see Teske, Bader & Golla, 2015 for Standard protocols), and the obtained chromatogram was edited using MEGA11: Molecular Evolutionary Genetics Analysis version 11 (Tamura, Stecher & Kumar, 2021). All the sequences were deposited in GenBank (NCBI; Benson et al., 2018) under the accession numbers ON062377 –ON062402 (see Table S3). No amplification product for Halichondria caruncula (Bowerbank, 1858) sensu Stephens (1915), Hymeniacidon stylifera (Stephens, 1915) and Hymeniacidon sublittoralis Samaai & Gibbons, 2005 from South Africa could be obtained.

Alignment, phylogenetic analysis, and Haplotype networks

The raw sequence data of the forward and reverse sequences were trimmed by removing uncertain sites. The forward and reverse sequences were then aligned using ClustalW (Thompson, Higgins & Gibson, 1994) in MEGA 11 (Tamura, Stecher & Kumar, 2021). Sequences were blasted in GenBank (Sayers et al., 2019) and the maximum score and E-values (Altschul et al., 1990) were used to select closely related specimens. The COI data set was checked for the potential occurrence of nuclear pseudogenes using the genetic code for invertebrate mitochondria, to detect frame-shift mutations, which would indicate that these sequences originate from a non–functional gene region, were identified. Sequences were compared to published data of Hymeniacidon sponges (see Table S3 in Turner, 2020), and thus sequences were jointly analysed with the Hymeniacidon data set used by Turner (2020) in MEGA 11 (Tamura, Stecher & Kumar, 2021). This included public sequences previously identified as H. sinapium and H. heliophila, which according to Turner (2020) are part of the single global species complex, H. perlevis (see Appendix S2–S6). The dataset compiled by Turner (2020), together with the South African sequences, was then used in the phylogenetic analysis and haplotype network. To see if the South African samples formed a distinct monophyletic clade in comparison to congeneric samples from other locations, we constructed a phylogenetic tree in MEGA 11 using Maximum Likelihood (ML) with the Tamura-3 parameter (T92), which was selected by the inbuilt model generator. Evolutionary distances were computed employing the Tamura-3 parameter (Tamura, 1992), and support for individual nodes was based on 1000 nonparametric bootstrap estimates (Felsenstein, 1985). The T92 distances were also used to compare levels of genetic differentiation between the sequences generated in this study and the published Hymeniacidon sequences (see Table S3 in Turner, 2020).

DnaSP 5.091 was used to evaluate haplotype (h) and nucleotide (π) diversities for individuals collected at the same location (Lourenço et al., 2017). Populations were divided into six groups according to world-wide presence (Turner, 2020; de Voogd et al., 2021). These were East Asia, North America (Pacific), North America (Atlantic), South America (Atlantic), Europe and South Africa. The South African populations were divided into three groups according to the national bioregional classification (see Sink et al., 2019). These include the Namaqua, Southern Benguela and Agulhas ecoregions. To determine how genetic variation is divided between groups, among locations within groups and within locations, the six and four groups described above were designated a priori following Lourenço et al. (2017).

A median-joining haplotype network was produced using the minimum spanning method (Bandelt, Forster & Röhl, 1999) as implemented in Popart (Leigh & Bryant, 2015). This analysis requires all included sequences to be the same length, so some sequences were trimmed whilst others were excluded. Alignments in the global dataset were 574 bp at CO1 (n = 115 sequences); the ITS alignment was 539 bp (n = 512). Alignments were longer when newly collected South African data were analyzed alone: 582 bp at CO1 (n = 29); 798 bp and ITS (n = 11). Sequence alignments were produced in Codon Code v.9 (CodonCode Corporation).

Ecological niche modelling of H. perlevis distribution

To gain a better perspective on the realized distribution of H. perlevis globally, ensemble species distribution modelling was applied (Appendix S1).

Occurrence data. For this purpose, occurrence/encounter data from multiple sources were compiled, including records from the World Porifera Database (WPD), GenBank, literature and South African observations. The location data were checked by Toufiek Samaai and Thomas Turner and only valid occurrence records were included in this study (Table S1). Bias, dubious and unverified data were excluded.

Because the data are occurrence/encounter data, background absence (pseudo-absence) data are required to apply standard correlative distributional models. Pseudo-absences were generated at random within the studies gridded spatial domain, with the thin layer of coastal area being globally generated. Given the environmental layers resolution (5 × 5 nautical mile grid cells), the first two adjacent grid cells were considered. The number of pseudo-absences generated were determined using the recommendations of Barbet-Massin et al. (2012), implemented in the R package SSDM (Schmitt et al., 2017), which was used in this study for ensemble species distribution modeling. Ensemble modeling was performed both with and without spatial thinning. Given that the encounter data are not globally uniform, spatial thinning, which is already implemented in the SSDM package, was used to deal with spatial bias (to reduce spatial bias due to non-random sampling while keeping most of the information).

Niche modelling. Multiple correlative statistical models are widely used to model the distributions of many taxa. The majority of these widely used models are already included in the SSDM package. For the ensemble modeling of the distribution of H. perlevis, eight correlative statistical models were considered: generalized linear model (GLM), generalized additive model (GAM), support vector machine (SVM), classification tree analysis (CTA), generalized boosted model (GBM), random forest (RF), multivariate adaptive regression spline (MARS), and artificial neural network (ANN). Appendix S1 contains a brief discussion of each of the eight distribution models. These models were trained on a randomly selected portion of the data (70%) and their prediction performance was evaluated using the hold-out set. Each model was fitted and evaluated four times to account for sources of variability due to random selection of training and evaluation sets as well as random selection of pseudo-absences. When analyzing classification models, multiple measures of performance can be used. For the purposes of this study, Kappa and area under the curve (AUC) were used. AUC usually has a value in the range of 0.5 to 1.

Models with AUCs of 0.5 are generally considered random classifiers, while those with values between 0.7 and 0.8 are considered fair classifiers, and those with values close to 1 are considered excellent classifiers (Kleinbaum & Klein, 2010). The inclusion of individual distribution models into the ensemble distribution modeling was based on whether the model had an AUC value greater than 0.7. To generate the ensemble species distribution map, all models with AUCs greater than 0.7 were pooled by weighting their predicted probability of occurrence by their AUC. The uncertainty map was also computed primarily to identify regions of high agreement and low agreement among the models considered, which correspond to low and high uncertainty regions, respectively. Uncertainty was computed as cell by cell variance in the predicted probability of occurrence by the models included in the ensemble. The response curves (partial effects) of each of the variables considered was generated by predicting the probability of occurrence for the variable of interest while keeping the remaining variables at their mean. This was done for each of the eight models considered (Appendix S1).

The intertidal area was delimited by extracting the coastal cells covering a range from −2 to 1 m from the General Bathymetric Chart of the Oceans (GEBCO) gridded bathymetric data set with a spatial resolution of 30 arc-seconds (http://www.gebco.net/).

Variable importance was computed on the holdout set. The amount of correlation changes between predicted values before and after permuting (reshuffling) a variable was used to measure its importance (expressed in percentage). ${\displaystyle I_r=1-max\left(cor\left(P_f,P_v\right),0\right)}$

where I_r is index of importance of a variable, Cor is correlation coefficient, P_f is prediction from the full model, P_v is prediction after permuting/reshuffling the variable v. Partial effect of each of the predictor were computed by predicting the response variable for the variable of interest while holding the other predictors at their mean.

Environmental variables. BIO-ORACLE was used to download global and readily available environmental layers. The environmental layers used in this study considered the minimum, maximum, mean, and range of surface temperature, surface salinity, and surface current velocity. Given the limitations of the environmental layers resolution (5 × 5 nm grid) and the fact that we are dealing with intertidal/coastal invertebrates, only grid cells within 10 km of the coastline were retained. The environmental variables were checked for multi-collinearity using the variance inflation factor (VIF) before the distribution of H. perlevis was modelled.

A VIF value of one indicates an absence of multi-collinearity, but larger values typically indicate the presence of a problem. Variables with VIF values > 5 are generally considered to be linearly related, and in the context of regression, their variance of the estimated parameters will be large, and its parameter will be poorly estimated (Hay-Jahans, 2011). Variables having VIF greater than 5 were thus excluded from this analysis. The final set of variables retained were: mean surface temperature, range surface temperature, mean surface salinity, range surface salinity, and mean surface current velocity. Appendix S1 shows the layers of environmental variables used in the distribution modelling.

The following model was fitted to model occurrence of the sponges. ${\displaystyle \begin{array}{cc}\hfill \text{Model}\hfill & \hfill \text{formula}\hfill \\ \hfill occurrence\hfill & \hfill spong{e}_{occ}\sim T_{mean}+T_{range}+S_{mean}+S_{range}+V_{mean}\hfill \end{array}}$

where sponge_occ is the occurrence is of Hymeniacidon perlevis; T_mean and T_range are the mean and range of coastal surface temperature respectively; S_mean and S_range are the mean and range of coastal surface salinity respectively; V_mean is the mean coastal surface current velocity.

All the analysis, visualization and report generation were done in R (R Core Team, 2021). Multiple R packages were utilized for data processing, visualization, analysis, and summary of results including (Alathea, 2015; Allaire et al., 2021; Henry & Wickham, 2020; Robinson, Hayes & Couch, 2022; Wickham, Chang & Henry, 2018; Wickham et al., 2021; Xie, 2021).

Material and acquisition

All recently collected voucher samples are housed at the Iziko Museum, Cape Town, South Africa under museum numbers SAMC-A091444–SAMC-A091463; MB-A094583–MB-A094599; MB-A094600–MB-A094614 (Table S3). Toufiek Samaai was granted permission to collect specimens during his various field excursions by the Department of Forestry, Fisheries, and Environment under Research Permits RES2014/DEA–RES2019/DEA.

Distribution of Hymeniacidon perlevis along the temperate southern African intertidal region

Hymeniacidon perlevis was found at 23 of the 53 locations sampled (Fig. 2; Table S2). The confirmed species range spans from Grosse Bucht, south of Lüderitz, Namibia on the southern African Atlantic coast to the Dwesa/Cwebe Marine Protected Area in the Agulhas region of South Africa’s Indian Ocean (east) coast. Kreiner et al. (2019a), Kreiner et al. (2019b) and Kreiner et al. (2019c) reported the species from Grosse Bucht, Diaz Point, Patrysberg, Mile 4, Badewane, and Möwe Bay in Namibia, but no samples were available to confirm the identifications, with the exception of the specimen collected from Grosse Bucht. The species represents the most conspicuous and common sponge on the intertidal shores of the temperate bioregion around the tip of southern Africa. The species was found mainly associated with intertidal rocky shores at locations both associated with high anthropogenic impact such as Saldanha bay, Table bay and False bay, all areas with commercial harbours and international shipping traffic, and remote, natural locations including the Tsitsikamma National Park, South Africa’s oldest Marine Protected Area.

We only found H. perlevis covered by sediment in one location (Strand), with surface projections extending beyond the sandy layer (Fig. 3G). Hymeniacidon perlevis was also found in high nutrient concentration areas throughout the Cape Peninsula (De Villiers, 2017; T Samaai, 2017, pers. obs.), in the kelp forest near Bettys Bay at a depth of 15 m, and on the west coast of South Africa by Stephens (1915) at a depth of 25 m. Although the species was not sampled between Bettys Bay and Knysna, it is present at various locations between these areas (Prof. G Branch, pers. comm., 2021).

Genetic analysis

Partial COI sequences were obtained for 29 specimens of H. perlevis (Table S3) with an alignment length of 691 base pairs. This data set reduced to two haplotypes differing by a single base pair (Fig. S1). The haplotype network of H. perlevis from South Africa is shown in Fig. 4. One haplotype was present only in the Namaqua region, while the other was shared across regions.

BLAST-n results revealed these sequences to be 99–100% identical to sequences previously identified as H. perlevis and H. sinapium from other regions (for example Portugal and California; note that H. sinapium is now considered a junior synonym of H. perlevis). This was further corroborated by phylogenetic reconstruction, with sequences of H. perlevis from South Africa grouping with those from other regions and H. sinapium in the maximum likelihood tree (Fig. S2), with weak maximum likelihood support. The H. perlevis sequences further showed about 97% similarity to H. flavia Sim & Lee, 2003 (EF217333.1, Korea), which is the most closely related outgroup species (Fig. S2) (Park et al., 2007). A haplotype network of this global sample is shown in Fig. 5. We employed a COI haplotype network to visualize the southern African specimens in the context of a global sample of publicly available data. One of the two haplotypes from southern Africa was among the most common haplotypes globally, also present in samples from Florida, Brazil, Portugal, Turkey, Korea, and California (Fig. 5).

In the global (539 bp) ITS alignment, four of the South African samples are identical to the most common haplotype, found in a sample from Northern California, five samples from Korea, and a large number of samples from Japan (Fig. 5). The other seven southern African samples differ from the most common haplotype by a single base pair (0.19% sequence divergence). Consistent with a previous analysis, the genetic variation within the sample of H. flavia is similar to the genetic variation in H. perlevis, despite the entire H. flavia sample being from Japan and Korea alone. When the South African ITS samples are analyzed alone, a longer alignment is possible (798 bp), and three variable sites are present (Fig. 4). One of these differentiates samples from Agulhas from the other regions.

Environmental niche modelling

Hymeniacidon perlevis has a very wide geographic distribution ranging from the Northern Hemisphere to Argentina, as well as South Africa to the middle latitudes of New Zealand. Visual exploration of occurrence data for H. perlevis are shown in Fig. 1. As can be seen in Fig. 1, most of the occurrence records are from coastal area around the UK, New Zealand and South Africa, with the remaining records from different parts of the globe, such as harbours, mudflats and inlets (Table S1).

Figure 6 depicts the predicted distribution, as well as the probability of occurrence, of H. perlevis from each of the eight models and the ensemble (Appendix S1). The ensemble prediction on the raw and thinned data is shown in detail in Fig. 6. The projected likelihood of occurrence was not significantly different for models based on thinning or raw occurrence data when seen visually. The ensemble produced with the best models resulted in an accurate overall description of H. perlevis distribution, including its expanding front (Fig. 6).

Along southern Africa, the niche model predicted a distribution further north into Namibia and in South America from southern Chile into northern Peru (Fig. 6). In addition, the prediction indicated that suitable habitat could potentially be found along southern Australia and the south Island of New Zealand. While the probability of H. perlevis being present in the Mediterranean and British Columbia shores was high, towards Northeast America the predicted likelihood decreased. No suitable habitat was detected along the tropical West African coast, the Indo-Pacific region, the Arabian Peninsula or India. Mean surface temperature was the most important predictor of the distribution of H. perlevis globally (Appendix S1), followed by the range of surface temperature. Most of the models performed reasonably well with AUC mostly above 0.8.

Taxonomy—species description

Systematic information with detailed morphological and spicule descriptions are provided below. The classification follows Morrow & Cardenas (2015).

Type species. Hymeniacidon caruncula Bowerbank, 1859: 286 (by subsequent designation; Bowerbank, 1864: 191) (this is considered a junior synonym of Hymeniacidon perlevis (Montagu, 1818: 86)).

Material examined. Table S4

Other material examined. Hymeniacidon perlevis voucher specimens. SAM-H4904 (Ts 305), Jacobs Bay, near Saldanha Bay (32 31′S, 17 30′E), depth 3–5 m, collected by T Samaai, 20 October 1997. Ts 329, Ts 331, Ts 337, Ts 338, Ts 343c, Elands Bay (32 20′S, 18 20′E), depth 3–6 m, collected by T Samaai, 15 November 1997. Ts 359, Ts 370, Ts 381, Ts 391, Groenrivier (30 29′S, 17 20′E), depth 3 m. Collected by T Samaai, 20 December 1997.

Hymeniacidon caruncula (Bowerbank, 1858) sensu Stephens (1915). NMSZ 1921.143.1443. Two fragments removed for loan. Station 479, False Bay shore 6th May 1904; Station 482, Saldanha Bay shore, 19 May 1904.

Leucophloeus stylifera Stephens, 1915. Syntype NMSZ 1921.143.1443. Two fragment removed for loan. Station 482, Saldanha Bay shore, 19th May 1904; Station 483, Entrance to Saldanha Bay, 45 m, 21 May 1904.

Hymeniacidon sublittoralis Samaai & Gibbons, 2005. Holotype. SAM-4903 (Ts 212), Vulcan Rock (34 04′S, 18 18′E), depth 27 m, collected by P Coetzee, 24 April 1996.

Description. A thin or thickly encrusting (Fig. 7A) to cushion-like sponge that varies greatly in form (Fig. 2). Diameter ranging from 5 cm long × 3 cm long × 4 cm thick to 14. five cm long × 8 cm wide × 6 cm thick; with processes of 1–4 mm high, 1–1.5 mm wide. In regions where there is considerably more wave exposure, this species is encrusting with a smoother surface. In sheltered or somewhat exposed locations, H. perlevis has upright processes that emerge from a basal mat. Surface variable, may be smooth and tuberculate, thrown into irregular folds, or covered with digitate processes. Oscules 0.5–1.5 mm in diameter, dispersed, level with the surface, or raised on low digitate processes. Firm, soft, fleshy texture that is compact and compressible. Live specimens of H. perlevis have distinct colour patterns being different shades of orange depending on geographical location (Fig. 2). Though this species can be blood red in other regions, in the current study, intertidal encrusting sponges of that colour were species of the genus Tedania (Fig. 2) (see also Samaai & Gibbons, 2005; Ngwakum et al., 2021).

Skeleton. The choanosomal skeleton, especially in the deeper regions, composed of a confused, disordered mass of styles, not organized into tracts (Fig. 7C). Towards the surface, tracts become ill-defined and with ascending fibres, ∼200 µm wide, with no separation between the primary and secondary tracts. The ascending tracts do not branch at the surface to form spicule brushes and tend to be vertically arranged. Numerous loose interstitial spicules. Large canals are present. The ectosomal skeleton consists of a dense tangential layer of spicules, ∼200–500 µm thick (Fig. 7C). Spongin scarce.

Spicules. Megascleres (Fig. 7B; Table 1; Fig. 8): styles smooth, straight, or slightly curved, thickest centrally, 250 (155–337) × 7 (7) µm, n = 20. Microscleres: absent.

Habitat and distribution. Found on the rocky intertidal areas, shallow subtidal reefs, mudflats on hard objects, harbours and inlets. Depth range 0–25 m (Table S1).

Status. Species may be native to southern Africa.

DNA barcodes. 691bp fragment of the universal mitochondrial cytochrome c oxidase subunit 1 gene, primer pair: LCO1490 and HCO2198 (Folmer et al., 1994). GenBank accession numbers ON062377 –ON062402& MT491492 –MT491502 (see Table S3). 539 bp fragment of the ITS gene, Genbank accession numbers MT501787 –MT501797 (see Table S3).

Remarks. The morphological features of H. perlevis and H. caruncula sensu Stevens (1915) were compared. Hymeniacidon caruncula sensu Stevens (1915) is similar to H. perlevis. In terms of spicule form, spicule size, external morphology, and colouring (Table S5), H. sublittoralis and H. stylifera differ from H. perlevis. Hymeniacidon sublittoralis is a thick, massive, erect, amorphous sponge, with numerous papillate processes that vary greatly in length. Surface smooth with various ridge-like structures, finely hispid and colour in situ yellow. Styles are large and thick with heads slightly subtylote, 394 (255–601) × 14 (14) µm (Samaai & Gibbons, 2005).

Hymeniacidon stylifera is easily distinguished from H. perlevis on spiculation; H. stylifera differs from H. perlevis by the larger size of the style megascleres, having a smooth surface and and have a very firm texture (Table S5).

The spicule size range of South African H. perlevis overlaps with specimens from Wales, Korea and Ireland (Fig. 8). A large spicule size range is found for the South African west coast and Wales specimens (see Fig. 8). California and Canada specimens have a smaller spicule range, with lengths very similar to New Zealand and Argentina (see Fig. 8). No spicule lengths are available for specimens from Portugal, Spain, France, Brazil, Florida, China and the Black Sea. The South African specimens conform to all previous descriptions in terms of morphology and spicule complements, as well as in habit (Bergquist, 1970; Ackers et al., 2007; Erpenbeck & Van Soest, 2002; Harbo et al., 2021).