**1. Introduction**

Scleractinian corals are known to host a variety of organisms belonging to different phyla [1–6]. Among coral ectosymbionts, the acoel flatworms of the genus *Waminoa* (Order Acoela, Family Convolutidae) have been studied only recently. Formerly included in the Platyhelminthes and now placed within the Acoelomorpha [7], *Waminoa* is considered an enigmatic group since much about its ecology and diversity remains unknown [8]. *Waminoa* spp. are found mainly on scleractinian corals, but also on octocorals, sea anemones, corallimorpharians, zoantharians, and echinoderms [8–15], and they show a circumtropical distribution, with the exception of the Caribbean [16]. *Waminoa* flatworms are characterized by the presence of intracellular dinoflagellate symbionts, which have also been found in the worm oocytes, suggesting that these symbionts are inherited via a vertical transmission and not obtained from the host coral [8–10].

**Citation:** Maggioni, G.; Huang, D.; Maggioni, D.; Jain, S.S.; Quek, R.Z.B.; Poquita-Du, R.C.; Montano, S.; Montalbetti, E.; Seveso, D. The Association of *Waminoa* with Reef Corals in Singapore and Its Impact on Putative Immune- and Stress-Response Genes. *Diversity* **2022**, *14*, 300. https:// doi.org/10.3390/d14040300

Academic Editor: Michael Wink

Received: 2 March 2022 Accepted: 12 April 2022 Published: 15 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

These flatworms subsist via different ways, by ingesting the host coral mucus [11,17] and/or by "standing" on the polyps of host anthozoans to obtain floating zooplankton [18]. For these reasons, *Waminoa* spp. could represent a threat to corals, although the worms do not appear to consume coral tissue directly [9]. In fact, acoelomorph flatworms may limit the host feeding on zooplankton by competing for the prey and by physically blocking the coral oral disc, possibly resulting in kleptoparasitism [18]. In this regard, *Galaxea fascicularis* polyps infested by worms showed a significant decrease of prey ingestion rates compared to polyps without worms and between 5 to 50% of total prey captured by the polyps was stolen by *Waminoa* individuals [18]. Since the coral mucus layer aids in heterotrophic feeding and represents a protective physiochemical barrier [19], the removal of the mucus by *Waminoa* spp. may reduce the coral's resistance to pathogens and environmental stressors [11,17]. In addition, being able to reach high densities and cover a significant portion of the coral [20], *Waminoa* spp. may also cause light shading affecting the coral's photophysiology, reducing the productivity of the coral holobiont [11,13]. Despite the possible negative effects of *Waminoa* worms on their hosts, specific studies to diagnose the health of the infested corals have never been performed so far.

As sessile organisms, corals rely on the modulation of their cellular and molecular mechanisms as the first defensive line against environmental stresses and invading pathogens, and these mechanisms represent useful biomarkers of corals' health status [21–24]. In this context, corals possess an innate immune system consisting of different self/non-self-recognition receptors, which activate specialized cellular and humoral signaling pathways leading to diverse downstream effector responses [25,26]. Among the complex network of immune processes, the complement pathway is a proteolytic cascade, by which pattern recognition receptors (PPRs), such as lectins, initiate intracellular signaling to enact complement component factors, such as C3-like proteins [27]. Lectins are recognition receptors that bind to glycans and play a role in non-self-recognition, cell-cell adhesion, bacterial cell wall recognition, and phagocytosis [28]. In particular, C-type lectins are a superfamily of Ca2+-dependent carbohydrate-recognition proteins involved in the activation of several innate immune responses [29–31]. Complement C3-like proteins, whose activation relies on lectins, are important in allorecognition and involved in the opsonization of the pathogens, chemotaxis, and activation of leukocytes [25,27]. Although complement-encoding genes, such as those of C-type lectins and C3, have been identified and characterized in corals [32–36], their involvement and modulation in response to biotic stressors remain poorly tested.

In addition to the immune response components, other diagnostic tools in corals able to reflect changes in cellular integrity and functionality caused by stress exposure are cellular proteins, such as Actin and Heat shock proteins. Indeed, Actin, which is a major cytoskeletal protein involved in cell motility, growth, and division, is thought to be a proxy of the growth rate in corals [37,38]. Heat shock proteins (Hsps) are molecular chaperones involved in cytoprotection and maintenance of protein homeostasis, and their expression is usually up-regulated when organisms face conditions that may affect their cellular protein structure [39]. For this reason, Hsps have been frequently adopted as cellular stress biomarkers in corals subjected to different environmental stressors [40–45]. In addition, Hsps may play a role in the coral immune system, since they can be activated in response to epizootic diseases or other biotic stresses [46–49].

The coral reefs of Singapore, an island megacity that has been experiencing intense urban development over the past 60 years, represent highly disturbed and urbanized coastal environments [50,51]. Coral communities here have been affected for decades by multiple chronic anthropogenic pressures, resulting in high levels of sedimentation and turbidity [52–55], and multiple bleaching events caused by climate change [56–58]. However, no work has examined the presence and the impact of *Waminoa* worms on Singapore's coral communities.

In this study, the association between *Waminoa* and scleractinian hosts was studied for the first time in Singapore reefs through an ecological survey and a molecular analysis.

In particular, we determined the prevalence of the association, the host range, the host preference of *Waminoa*, and the flatworm abundance on the coral surface. In addition, we tested a possible effect of *Waminoa* on components of the coral immune system and cellular stress response. For this purpose, the coral *Lobophyllia radians* was selected being a species highly colonized by the *Waminoa* in Singapore, and the gene expression of *C-type lectin*, *C3*, *Hsp70* and *Actin* were examined and used as biomarkers to provide new insights on the nature of the association.

#### **2. Materials and Methods**

#### *2.1. Ecological Analysis*

#### 2.1.1. Study Area and Sampling Design

Fringing reefs of two islands south of mainland Singapore, Pulau Hantu (1◦2274 N, 103◦7465 E) and Kusu Island (1◦2255 N, 103◦8585 E), were selected as study sites (Figure 1). Both reefs are characterized by a shore-adjacent reef flat leading seaward to the reef crest and down the reef slope to ~8–10 m maximum depth because of high levels of suspended sediments that cause extreme light attenuation [55,59].

**Figure 1.** Map of the study area showing the two investigated sites. Pulau Hantu is a sheltered site situated 8 km south of mainland Singapore, while Kusu Island is located 6.4 km south of Singapore's city center and 13.4 km east of Pulau Hantu.

In both sites, extensive surveys were conducted by SCUBA diving between August and October 2019 to detect the occurrence of *Waminoa* individuals on different coral genera (Figure 2). Images of *Waminoa* specimens were taken by a digital camera and analyzed. We distinguished a single morphotype of *Waminoa* in the surveyed area (as described in [8,15], Figure 2) and we treated it as *Waminoa* sp.

**Figure 2.** *Waminoa* sp. individuals in association with different coral genera in Singapore reefs. (**A**) Close up of a single *Waminoa* sp. individual characterized by an obcordate general shape, body with brown coloration, white peripheral outer edge, white random dots, and white comparatively larger one white internal spot. Flatworms colonizing corals belonging to the genera *Pachyseris* (**B**), *Fungia* (**C**), *Lobophyllia* (**D**,**E**), *Goniastrea* (**F**), *Pectinia* (**G**), *Merulina* (**H**), and *Favites* (**I**). Scale bars: 1 mm for (**A**) ~1 cm for (**C**,**E**,**F**,**G**) ~2 cm for (**B**,**D**,**H**,**I**).

Six 50 × 2 m belt transects (100 m<sup>2</sup> each), spaced 10 to 20 m apart and placed parallel to the coast at a constant depth between 5 and 8 m (depending on the tides) along the reef slope, were randomly laid at each site. Within belt transects all the coral colonies were identified to genus level, according to Huang et al. [60] and Wong et al. [61], and the number of colonies of each genus found colonized by *Waminoa* sp. was recorded. Moreover, in coral colonies hosting *Waminoa* individuals, their abundance was estimated by determining the percentage of coral surface covered by flatworms and was indicated with four coverage categories, as suggested in [62]: low (coral surface covered 1–10% by worms), moderate (11–25%), severe (26–50%), and extreme (>50%).

In each belt transect, the point intercept transect (PIT) method was also performed (by recording data every 10 cm [63]) to determine the composition and structure of the benthic community, as well as the cover percentage of each coral genus. Data were collected using the following benthic categories: algae, dead coral, coral, coral rubble, rock, sand, and other (sponges, soft corals, tunicates, zoantharians, and unknown). Furthermore, for the macro-category "coral", the genus was also recorded.

#### 2.1.2. Data Analysis

For each transect, the prevalence of the association was calculated as the ratio between the number of corals colonized by *Waminoa* sp. and the total number of colonies. By averaging the corresponding prevalence values measured on the six random belt transects for each site, both an overall and a series of taxon (coral genus)-specific prevalence values were determined. Data normality was verified using the Shapiro–Wilk test. A one-way ANOVA was used to test significant differences in the overall association prevalence between the two sites analyzed. The same analysis, followed by Tukey's honestly significant difference (HSD) post hoc tests for multiple pairwise comparisons of means, was performed to assess significant differences in the prevalence of the *Waminoa* sp.–coral association among the different host genera. Coral genera showing a prevalence < 5% were not included in the analysis.

The host preference of *Waminoa* sp. in terms of coral genus was tested through the Van der Ploeg and Scavia Selectivity coefficient (Ei), following [64]. This coefficient is defined for a group *i* as:

$$\text{Ei} = \frac{\left[Wi - \left(\frac{1}{n}\right)\right]}{\left[Wi + \left(\frac{1}{n}\right)\right]}$$

where *Wi* represents the value of Chesson's α and *n* represents the number of habitat types [63], here represented by the coral genera found in the study area. Chesson's α value (*Wi*) is defined as:

$$Wi = \frac{r\dot{i}}{Pi} / \sum\_{\dot{i}} r\dot{i} / Pi\dot{i}$$

where *ri* represents the frequency of a habitat category (coral genus) in the environment, and Pi represents the frequency of the same habitat category in which the organism of interest (*Waminoa* sp.) is found [65]. Values of selectivity coefficient range between −1 and 1, with −1 meaning complete avoidance of a host coral genus, and 1 meaning exclusive preference for a specific coral genus [66].

A one-way ANOVA followed by a Tukey's HSD post hoc test was performed to assess differences in the coverage percentages of the flatworms on the surface of the host corals, between the four coverage categories.

#### *2.2. Molecular Analysis*

#### 2.2.1. Coral Collection

To assess the effect of *Waminoa* sp. on coral's gene expression, five colonies of *Lobophyllia radians* showing patches of *Waminoa* sp. individuals in a fixed position on their surface were selected, monitored for several days and sampled on Kusu Island. All five colonies were moderately (11–25%) covered by the worms. For each colony, fragments of approximately 2 cm<sup>2</sup> were collected using a hollow-point stainless steel spike [67]. Two coral fragments were sampled from each colony: one fragment located just underneath a patch of *Waminoa* sp. individuals (marked as "W"), and the other from a colony portion without worms, at least 5 cm away from the *Waminoa* individuals (marked as "W/0"). In addition, three healthy colonies of *L. radians* not colonized by the flatworms were randomly sampled as controls to test primer efficiencies.

All coral samples were taken at the same hour (around 9:00 a.m.), at the same shallow depth, and during high tide, to minimize any possible effects of abiotic variables, such as water temperature and light intensity, on the gene expression. During the sampling period, the sea surface temperature was continuously logged and its mean was 30.36 ± 0.39 ◦C, with very slight oscillation between 29.66 and 30.42 ◦C. In "W" fragments, *Waminoa* individuals were removed from their hosts by using the pipettes as previously described [8] and the morphological condition of the coral tissues, as well as the presence of any physical damages or lesions, were evaluated. All the coral portions were immediately placed in pre-labeled tubes and, at the end of the underwater sampling, the seawater

was decanted and replaced with RNAlater (Thermo Fisher Scientific, Singapore) with a tissue:RNAlater ratio of 1:10. The maximum time between collection and placement in RNAlater was about 30 min. The sample tubes were inverted to mix for 30 s and kept at 4 ◦C overnight to allow complete penetration of RNAlater into the coral tissues. The tubes were subsequently stored at −80 ◦C until RNA extraction.

#### 2.2.2. RNA Extraction and REVERSE Transcription (RT)

Total RNA was extracted from all the coral samples without homogenization to reduce RNA fragmentation using TRIzol Reagent (Life Technologies, Sigma-Aldrich, Singapore) following the manufacturer's protocol. RNA quality was checked by examining with gel electrophoresis for presence of clear bands of ribosomal RNAs and RNA concentration was estimated using Qubit (RNA Broad Range Assay Kit, Thermo Fisher Scientific). Complementary DNA (cDNA) was immediately prepared from 1 μg of total RNA for each sample, using a one-tube format of iScript RT Supermix (Bio-Rad, Hercules, CA, USA) for a reverse transcription quantitative polymerase chain reaction (RT-qPCR). Reaction setup was composed of iScript RT Supermix (4 μL), RNA template (varied depending on RNA sample concentration; 14.6–200 ng/μL), and nuclease-free water (variable), with a final volume of 20 μL as previously described [68]. Incubation of the reaction mix was performed in a Labcycler (Sensoquest, Göttingen, Germany) following the manufacturer's protocol: priming for 5 min at 25 ◦C, RT for 20 min at 45 ◦C, and RT inactivation for 1 min at 95 ◦C.

#### 2.2.3. Primer Design and Validation

A local search in the *L. radians* transcriptome previously sequenced (raw sequencing data available at NCBI Sequence Read Archive under accession number PRJNA512601) and assembled by [69] using BLASTn against the GenBank database was performed. Coral transcripts matching genes from the genomes of *Acropora digitifera*, *Orbicella faveolata*, or *Stylophora pistillata* were identified and orthologous genes on the *L. radians* transcriptome selected. The accuracy of the sequence of each gene of interest (GOI) was checked using the NCBI Nucleotide BLAST tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 15 November 2019). A megablast search was performed to ensure that the sequences producing significant alignments with the query sequence were corresponding to the GOIs. A BLASTn search was performed in an open-access coral genomic database at http://reefgenomics.org/blast/ [70,71] with *Goniastrea aspera* genome as subject and the sequence of each GOI as query. A suitable region within exons was selected and the query corresponding to the longer hit with the highest identities value was selected to design primers. Primers for each gene were designed using the online tool by NCBI that incorporated Primer3 and BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast). The set of primers given as output by the NCBI tool were analyzed and used in http:// reefgenomics.org/blast/to perform a BLASTn search against Symbiodiniaceae nucleotide databases (Supplementary materials Table S1) to verify primer specificity for coral DNA. The selected set of primer was then used for a megablast search using the NCBI Nucleotide BLAST tool to verify the absence of significant alignment with marine species found in the same environment as *L. radians*. The designed primers are shown in Table 1. The specificity of the selected primers for *L. radians* was tested by performing PCR using the GoTaq Green Master Mix in a LifeECO Thermal Cycler (Bioer Technology, Hangzhou, China). PCR reaction steps were (1) denaturation: 95 ◦C for 45 s, (2) annealing: 50–55 ◦C for 45 s, and (3) extension: 72 ◦C for 3 min, repeated for 35 cycles.


**Table 1.** List of genes of interest (GOIs) with primer designs. The melting temperature (Tm), the % of guanine and cytosine (GC), and the length of PCR product (bp) are also reported for each gene. Additional information of these sequences can be found in Supplementary Materials.

To test primer efficiencies, a series of twofold dilutions of *L. radians* cDNA starting from 5 ng/μ<sup>L</sup> were performed for the samples collected from the control colonies [38]. Each dilution was used in triplicate for each primer to assess the primer efficiency through an RT-qPCR using the CFX96 Real-Time PCR System (Bio-Rad). Calculations of efficiencies (E, the amplification factor per PCR cycle) needed to correct for amplification efficiencies per primer were undertaken using an MCMC.qpcr package in R developed by Matz et al. [72]. The function, PrimEff( ), calculates E and plots the regression slopes and E based on dilution series. GOIs with E values outside the 1.85–2.16 range had primers redesigned and re-validated. GOIs that failed amplification were excluded from downstream analyses [68].

#### 2.2.4. Gene Expression Quantification

RT-qPCRs were performed in a CFX96TM Real-Time PCR System (Bio-Rad) using SsoAdvanced inhibitor-tolerant SYBR Green Supermix following the manufacturer's protocol (polymerase activation and DNA denaturation: 3 min at 98 ◦C, denaturation: 15 s at 95 ◦C, and annealing/extension: 30 s at 60 ◦C) repeated for 40 cycles. The reaction mix was prepared as in Table S2. Each sample was tested in duplicates for each of the four genes. To control for variations in expressions of genes caused by differences in RNA concentration of each sample, the amount of cDNA template was standardized to ~10 ng of cDNA for every reaction mix.

#### 2.2.5. Data Analysis

Data obtained from RT-qPCRs expressed as "cycle of quantification values" (i.e., Cq values) were collated and sorted for subsequent analysis. RT-qPCR data were analyzed using generalized linear mixed models based on lognormal-Poisson error distribution, fitted using the MCMC.qpcr package Version 1.2.3 in R Studio Version 1.1.463 as previously reported [72]. Molecule count data with corrections for primer efficiencies were derived with amplification efficiencies (E) per gene and Cq for a single target molecule using the formula: Count = E(Cq1 − Cq). The Cq-to-counts conversion is the key transformation in this method, in which higher variation at the low gene expression values is properly accounted for by the relative quantification model. The transformation makes it possible to fit the resulting data to generalized linear mixed models to account for Poisson-distributed fluctuations when the number of the molecule count is low. Similar Bayesian approaches for analyzing qPCR data have been used in several other reports [68]. Results of the mcmc.qpcr() function were then plotted with HPDsummary() to visualize fold changes in gene expression in response to the presence of *Waminoa* worms. HPDsummary() also calculates all the pairwise differences between treatments and their statistical significance for each gene. Each gene profile was examined to determine the differential expression level between coral samples underneath the surface of the coral colonized by *Waminoa* sp. and samples of the same colonies at least 5 cm apart from the flatworms.

A multivariate analysis was performed using the statistical package PRIMER-E v.7 with the PERMANOVA+ add on [73,74] to investigate together the modulation of all biomarkers in response to *Waminoa* sp. colonization. In particular, data related to the levels of all the biomarkers were square root transformed to calculate a matrix based on the Bray–Curtis similarity. To test for differences in biomarker levels between corals with and without worms, a non-parametric permutational multivariate analysis of variance (PERMANOVA) was performed using 999 permutations with partial sum of squares and unrestricted permutation of raw data. Values were considered statistically significant at *p* < 0.05.
