Cleaning symbiosis in coral reefs of Jardines de la Reina National Park

Study area

The Jardines de la Reina archipelago (Fig. 1) currently has more than 950 km² under the National System of Protected Areas in Cuba, with the management category of National Park (Pina-Amargós et al., 2014) and consists of an extensive chain of cays associated with mangrove islands and stretches of coral reef facing the Caribbean Sea. The reef structure in this area includes a reef terrace between 5 and 15 m deep, followed by a rocky slope culminating in a sandy esplanade, with the formation of large spur and groove area and a reef drop-off (Navarro-Martínez et al., 2022).

The sampling effort was concentrated on the largest group of cays in the archipelago, Las Doce Leguas, during March 2019, December 2021, and August 2022. A total of 26 sites continuously distributed from the far-East to the far-West of the MPA were sampled. For analysis purposes, sites were grouped into three main regions according to different levels of protection effectiveness: low protection level (eastern region), medium protection level (western region) and high protection level (central region) (Fig. 1, Table 1). This artificial grouping is similar to that applied by Pina-Amargós, González-Sansón & Cabrera-Páez (2008) and Pina-Amargós et al. (2014) and facilitates the evaluation and comparison of the results. Studied reefs are located south of the line of cays and include two ecosystems: the crest reef and the reef cliff. The sampling depth did not exceed 25 m.

Sampling methodology

Between one and four band transects of 50 m × 4 m were conducted at each site, representing a total of 19,200 m² monitored (n = 96). In each transect, the number of cleaning stations, the cleaner species and clients associated, and the number of individuals of each species was quantified.

Cleaning stations were considered as such when the reef sites that met the typical characteristics and conditions (i.e. coral, sponge or rock promontories, crevices, large anemones) and where occurred the observation of pairs of client and cleaner species interacting or only by the presence of species of well-known role as dedicated cleaners (Limbaugh, 1961). The client and cleaner diversity in each transect is the number of species of each role detected on transects associated with cleaning stations. The total abundance of cleaners and clients is the number of individuals, regardless of species, detected acting these roles associated with cleaning stations.

A team of two SCUBA divers conducted cleaning station surveys per transect. Measurements were started approximately five minutes after defining the transect positioning to minimize behavioral disruptions. Cleaning stations and individuals were counted as long as they were observed in front and to the sides of the observer, and most species were identified in situ.

Video recordings were made to identify the species involved and the number of times each pair of species interacted. Between 40 and 60 min of video recording was collected for each region. According to the criteria of Losey (1972) and Quimbayo et al. (2018b), we considered that there were cleaning interactions when physical contact between client and cleaner was detected. Client species were classified in functional groups following the criteria of Micheli et al. (2014) and Navarro-Martínez et al. (2022) with modifications according to Choat, Clements & Robbins (2004), Bonaldo, Hoey & Bellwood (2014) and Adam et al. (2015) and threat category according to the International Union for Conservation of Nature (IUCN). Functional groups were based on the combination of two functional traits: mobility and feeding guild as previously described in Navarro-Martínez et al. (2022) (Table S1). Specifically, the mobility categories were pelagic roving, cryptic, and midwater. Feeding guilds were based on nine categories: macroinvertivores (invertebrate feeders of 25–50 cm maximum length), microinvertivores (small invertebrate feeders of < 30 cm length), sessile invertivores (sessile invertebrate feeders), predators (invertebrate feeders and piscivores), piscivores, omnivores, planktivores, browsers, and grazers.

Statistical analysis

To assess the existence of a possible effect of the protection gradient on symbiotic cleaning interactions, MPA regions were compared according to number of cleaning stations, diversity of cleaner species, diversity of client species, total abundance of cleaners and total abundance of clients. This was done by using a smaller sample size (N = 37) than the initial one, so that three regions had an approximately similar number of sample units. Variables were tested for homoscedasticity (Levene test) and to verify if they fitted to a normal distribution (Kolmogorov–Smirnov test). Non-parametric Kuskal-Wallis test was used to verify statistical differences (α = 0.05) among regions.

Mutualistic interaction matrices between client species and cleaners were generated for each sector and the whole archipelago. Bipartite interaction networks were constructed for each region and for the MPA. The networks were generated from adjacency matrices with values for the number of interactions between pairs of species. The networks, due to their bipartite structure, consist of two groups of species segregated in different columns. Each column is formed by nodes, representing the species, the link between nodes represents the interaction and the thickness of the link is proportional to the number of interactions.

We calculated connectance, nestedness, H₂′ specialization (Blüthgen, Menzel & Bluthgen, 2006), niche overlap and complementarity for all networks. Connectance is calculated as proportion of realized interactions given a total number of possible links (Quimbayo et al., 2018a). Nestedness was assessed using NODF metric and it describes asymmetrical networks with some species cleaning the most of clients whereas others clean just a subset of those clients (Guimarães et al., 2006). H₂′ measures network-level specialization ranging between 0 (no specialization) and 1 (perfect specialization (Dormann et al., 2009). Niche overlap indicates similarities in interaction pattern between species from the same trophic levels whereas functional complementarity measures the ecological niche complementarity between species of the same trophic level as the total branch length of a dendrogram based on differences in interacting species assemblages (Devoto et al., 2012). These indices are generally used to describe the characteristics and structure and infer the robustness and the specialization of interaction networks (Quimbayo et al., 2018a). Processing was performed using R programming environment and bipartite package (R Core Team, 2023; Dormann, Gruber & Fründ, 2008).

Interaction ecology

During our surveys, we recorded 45 species engaging in cleaning mutualism interactions on the reefs of the Jardines de la Reina National Park. Of these, 39 species were identified as clients, including the taxonomic entities Pterois spp. and Kyphosus spp. whose in situ identification to species level proved difficult. Client species belong to 15 families (Table S1) and those with higher number of species were Labridae (eight species), Serranidae (seven species), and Pomacentridae (four species). Six species belonging to two families of fish: labridae (two species), and Gobiidae (two species); and two families of shrimps: Palaemonidae (one species), and Stenopodidae (one species) were identified as cleaners. In situ identification of Elacatinus genie and E. evelynae without disturbing cleaning behaviors was difficult, so they were treated as the Elacatinus spp. (Table 2).

During the study, 69 different combinations of cleaner –client pair species were detected. Cleaner fish, Bodianus rufus (Labridae) interacted with a total of 23 different client species, the Elacatinus spp. with 20 species, Thalassoma bifasciatum with 17, and the shrimps Stenopus hispidus and Ancylomenes pedersoni with six and three species, respectively. In the case of the clients, the species that interacted with the highest number of cleaner species were the serranids Epinephelus striatus and Mycteroperca tigris with five and four species respectively.

In terms of number of interactions, a total of 588 interactions were recorded. The most frequent cleaner species was B. rufus (Labridae) with 39% of interactions, followed by juveniles of T. bifasciatum (Labridae) with 35% of interactions, Elacatinus (Gobiidae) with 21% and the decapods A. pedersoni (Palaemonidae) and S. hispidus (Stenopodidae) whose interactions together represent 2% of the total (Fig. 2). According to Vaughan et al. (2017) three of the detected species are considered dedicated cleaners; these are A. pedersoni, Elacatinus spp., while the remaining cleaners (T. bifasciatum, B. rufus, and S. hispidus) express this behavior facultatively (Table 2).

The most frequently detected client species establishing interactions were Clepticus parrae (Labridae) with 30% of the total, followed by Caranx ruber (Carangidae) with 12%, and Acanthurus coeruleus (Acanthuridae), and Epinephelus striatus (Serranidae) both with 7%. These species constituted more than 50% of the total interactions detected in the study. C. parrae and T. bifasciatum, with 21% of the total interactions, constituted the cleaner- client pair of species that interacted the most; followed by C. ruber and B. rufus whose interactions accounted for 7% of the total. These were followed by the pairs C. parrae and B. rufus, A. coeruleus and B. rufus, and E. striatus and Elacatinus spp. with 7%, 5%, and 4% of the total interactions respectively.

We identified 15 different client functional groups. Roving macroinvertivores presented the highest number of species associated with cleaning stations (with seven species), followed by roving grazers (six) and roving piscivores (five) (Table S1). Regarding the frequency of interactions established by each group, the highest percentage of interactions was presented by the group of midwater planktivores with 30% of the total interactions, followed by roving herbivores with 25% and then, midwater piscivores with 11.5%.

Analysis of the number of interactions of cleaner species with the functional groups of clients revealed that B. rufus interacted with clients from 15 functional groups, and more than 50% of its interactions were established with roving grazers and midwater piscivores and planktivores. T. bifasciatum interacted with species from 11 functional groups; however, almost 90% of the time it cleaned midwater planktivores and roving grazers. The Elacatinus spp. clients belonged to 10 different functional groups and, although the number of interactions is more homogeneously distributed than in the rest of the species, in about 90% of the observations cleaning gobies interacted with piscivores, predators and roving grazers. Shrimp species S. hispidus and A. pedersoni interacted with species from seven and three functional groups, respectively. In both cases, roving piscivores predominated.

Effect of protection gradient on cleaning symbiosis

A total of 230 cleaning stations (N = 96) were observed throughout the study, with a mean density value of 2.3 ± 0.5 stations/200 m². This variable presented a similar pattern throughout the entire area of the PNJR and did not show statistically significant differences among MPA regions (p = 0.79, H = 0.46, Fig. 3A). However, highly protected region exhibit the highest medium density value of cleaning stations with 2.7 ± 1.3 stations/200 m², whereas east and west regions show similar densities of stations.

In the case of cleaner abundances, the differences found among regions were not statistically significant (p = 0.31, H = 2.43, Fig. 3B). Medium protection region exhibited the lowest mean abundances, whereas high and low protection region presented very similar values of abundances with 9.6 ± 1.7 and 9.5 ± 1.6 individuals respectively. Client abundance was also statistically similar across three regions of MPA (p = 0.43, H = 1.67, Fig. 3C). However, this variable presented the same pattern observed in number of cleaning stations, high protection region exhibited the highest mean abundances with 10 ± 2.4 individual clients, whereas the lowest mean client abundance values was registered in the lowly protected region with 6.5 ± 2 individuals.

The diversity of cleaners (p = 0.95, H = 0.10, Fig. 3D) and clients (p = 0.62, H = 0.95, Fig. 3E) was also not affected by the level of protection. However, we noticed that the high protection region had the highest mean values for client species number, with a mean of 3 ± 1.5 species in each transect and records of up to 6 client species.

The networks differed in their degree of specialization and nestedness according to the level of protection of the regions. The network of the high protection region had the highest H2′value (0.55) and the lowest connectance value (0.33), indicating that it was the most specialized network. The network of the low protection region had the lowest H2′value (0.41) and the highest connectance value (0.53), indicating that it was the least specialized network. The JRNP network had intermediate values for both metrics (0.39 and 0.54, respectively). All networks had similar nestedness values, with the JRNP network having the highest NODF value (54.9) (Table 3). Functional redundancy and niche overlap are measures of how similar the species are in their roles and interactions within the network. Functional redundancy reflects the extent to which species can be substituted by others without affecting the network structure or function, while niche overlap measures the degree to which species share the same resources or partners in the network. The network of the high protection region had the lowest niche overlap value (0.23) and the highest functional complementarity value (262.4), suggesting that its species were more distinct and less replaceable than those of the other networks.

Interaction ecology

The 45 species of fish and crustaceans recorded in this research have been previously reported for the JRNP and fish represent 17% of the total species (Pina-Amargós, González-Sansón & Cabrera-Páez, 2008). Cleaner species detected have been previously reported to occupy that role (Limbaugh, 1961) (Table 3). Most studies in the Caribbean have noted the importance of Elacatinus and juvenile individuals of T. bifasciatum in cleaning stations (Longley & Hildebrand, 1941; Losey, 1972; Darcy, Maisel & Ogden, 1974; Dunkley, Cable & Perkins, 2018). However, in this study B. rufus and T. bifasciatum were the main cleaners.

According to Guimarães et al. (2006), the ratio between species richness of clients and cleaners is the typically observed in reef cleaning networks. The total number of species involved in the interactions (45) indicates that PNJR cleaning mutualism network is a species-rich network, whose number of species in both groups exceeds those described in others in the region with similar sampling efforts (Johnson & Ruben, 1988). However, both cleaner and client richness could increase since that studies in Caribbean region (Titus et al., 2019; Dunkley et al., 2019) report interactions of cleaners detected in this research with client species inhabiting PNJR reefs that were not observed in cleaning stations. For example, species such as Pomacanthus paru, present in the MPA reefs (Pina-Amargós, González-Sansón & Cabrera-Páez, 2008) and frequently reported as cleaner in their juvenile phase (Vaughan et al., 2017), were not observed expressing cleaning behaviors during the study, which indicates more sampling effort is required.

Effect of protection gradient on cleaning symbiosis

The absence of significant differences among regions in the variables related to cleaning symbiosis may be explained by the fact that the protection gradient does not have a strong impact on the fish species that participate in this mutualistic interaction. These results are similar to those reported by Silvano, Tibbetts & Grutter (2012), who found that cleaner fish density, as well as the number of individuals and species cleaned, did not differ significantly bet ween a fished and a no-take site on the Great Barrier Reef. Therefore, the abundance and diversity of cleaners and clients may be more influenced by other factors that are independent of the level of protection, such as habitat complexity, benthic community composition, ectoparasite infection rates, and interspecific interactions.

Habitat complexity is one of the main factors that determine the number and distribution of cleaning stations, as they provide suitable substrates and shelters for cleaners and clients (Whittney et al., 2021). The coral reefs sites of JRNP have similar habitat complexity and benthic community composition throughout the park, as reported by Pina-Amargós, González-Sansón & Cabrera-Páez (2008) and Hernández-Fernández et al. (2019). This may explain why the density of cleaning stations did not vary significantly among regions. However, the high protection region had a slightly higher density of cleaning stations, which may be related to the higher abundance of fish in this region (Navarro-Martínez et al., 2022). Higher fish abundance may imply higher ectoparasite infection rates, which increase the demand for cleaning services and the availability of food resources for cleaners (Marcogliese, 2002; Sikkel, Cheney & Côté, 2004). This may enhance the trophic niche and abundance of cleaners, especially those that are dedicated cleaners and prefer ectoparasites in their diet, such as Elacatinus cleaning gobies (Soares et al., 2008; Soares et al., 2010).

We also expected to find higher abundance and species richness of clients in the highly protected region. However, we did not observe significant differences in these variables across the protection gradient. One possible explanation is that the competitive interactions among clients and cleaners at the cleaning stations may have masked the effect of protection. Cleaning stations are not homogeneous in size and structure, but they depend on the species of cleaners that inhabit them (Huebner & Chadwick, 2012; Whittney et al., 2021). However, this variation does not seem to influence the number of cleaners (Whiteman & Côté, 2004) or clients that use the stations at any given time (either as cleaners or as recipients of cleaning services). Therefore, the number of clients at the stations may be regulated by the availability of fish from the surrounding community that may seek cleaning services. This is consistent with the finding that the abundance of clients and cleaners is the main predictor of the visitation rate to the cleaning stations (Floeter, Vázquez & Grutter, 2007). Another factor that may limit the number of clients at the stations is the presence of predators or competitors that may deter potential clients from approaching the cleaners. This is supported by our results that showed that piscivorous and predatory species were among the most frequent visitors to the cleaning stations.

Another potential reason why we did not find a clear effect of protection on the variables related to cleaning symbioses is that protection may only affect some specific species that are targeted by fishing activities (Pina-Amargós et al., 2014). These species may not be involved in the cleaning mutualism network, or may have a minor role in it. Furthermore, the cleaner-client networks are not tightly structured or obligatory (Quimbayo et al., 2018a), so the dependence between the two groups is low. Therefore, protection may not have a strong impact on the dynamics of cleaning symbioses in coral reefs. However, we must keep in mind that this work evaluated different degrees of protection in a marine protected area, where even the least protected area also has protection. It would be essential to evaluate and compare with marine areas without any protection and with a marked effect of fishing (Silvano, Tibbetts & Grutter, 2012).

Network indices

Indices calculated for the networks in each region describe them in two important ways: structure and functional redundancy. In the first group, connectance provides information on how connected the two trophic levels that make up the network are in relation to the maximum value of possible connections between the nodes of both levels (Dormann et al., 2009). Higher connectance values mean that the network is more complete. We found that the network in the low protection region had the highest connectance, even though we expected the opposite. This may be related to the species richness in each region. More species means more nodes in the network, and therefore less chance of observing all the potential interactions. Since we recorded more species in the high protection region (23 species), we would expect that region to have lower connectance than the others. Nestedness shows how similar the patterns of interactions are among different subsets of species in the network (Bascompte & Jordano, 2013). Higher nestedness values indicate that there is a core group of species at each level that interact frequently with each other (generalists), and another group of species that interact less often and only with the core group (specialists). This way, the generalists provide a stable source of resources, in this case, cleaning services, that allows for specialization. Our results for nestedness are consistent with those reported for other cleaning networks in the Caribbean (Quimbayo et al., 2018a).

We also calculated three indices that reflect the functional redundancy of the network: H2′specialization, niche overlap and functional complementarity. These indices can be analyzed together, as they all indicate how similar or different the species are in their interactions with each other. In this case, the values we obtained for the archipelago network suggest that, as usual for these kinds of networks, there is a high degree of generalization and functional redundancy, which agrees with the findings of Quimbayo et al. (2018a). These authors argue that cleaning mutualism networks exhibit a high degree of niche overlap and generalization, due to the non-specific nature of the interactions. In the JRNP mutualisms, this is due to the predominance of interactions established by the facultative cleaners T. bifasciatum and B. rufus. These species exploit other food sources and often clean up clients that are common or use the same habitat. In localities with only facultative cleaners, client use overlaps more because they clean when they are juveniles or when predation risk is low (Côté, 2000; Vaughan et al., 2017) and therefore eventually and sporadically interact with the entire client group.

Our results also suggest that the network in the high protection region is slightly less generalized and redundant, and more specialized and complementary. This may be related to the higher species richness and abundance of clients, which may create better conditions for cleaners and clients to show more natural and specific behaviors and preferences related to cleaning symbiosis that are not possible in environments with fewer clients. Moreover, Silvano, Tibbetts & Grutter (2012) showed that the community of clients changed from a no-take site to a protected area. Therefore, cleaning symbiotic network structure may also vary according to the level of protection, as different client species may form different types and frequencies of links with cleaners.