CC-BY 4.0
Introduction
The resilience of wild populations to environmental change is closely linked to their ability to adapt (Eizaguirre & Baltazar-Soares, 2014). Several taxa show such ability in response to climate related changes of their environment, e.g., in plants (Jordan et al., 2017), birds (Gienapp et al., 2008), mammals (Boutin & Lane, 2014), fish (Muñoz et al., 2015) and corals (Matz et al., 2020). Effective evaluation of this ability is successfully integrated in conservation strategies (Funk et al., 2019). For instance, the efficiency of coral reef restoration is logically maximized when the ability of species to adapt is known and high (Baums et al., 2019). Assessing the relative environmental vs. genetic contribution to fitness variation, allows identifying the mechanisms shaping the ability of wild populations to adapt because the genetic variation of fitness is a prerequisite for adaptive evolution. This assessment is typically done by using quantitative genetic approaches to estimate genetic and environmental variance components for fitness proxies (Kruuk, 2004). Experimental evidence in the wild showed that the Lifetime Reproductive Success (LRS) is a good proxy for the real fitness of individuals in wild populations (Brommer et al., 2004). LRS is the total number of offspring captured on site in a given monitored population recorded as breeders produced by individuals over their lifetime (Hendry et al., 2018). However, it remains rarely used in wild populations because its measurement requires long-term datasets (Bonnet et al., 2022). It was studied in birds and mammals (e.g. Kruuk et al., 2000; Mccleery et al., 2004; Teplitsky et al., 2009; Mcfarlane et al., 2014), and only rarely in other taxa such as insects, fish and plants (Hendry et al., 2018; Salles et al., 2020). In our study, the LRS definition was slightly different because it was based on recruits defined as adults (or subadults) sampled and born on site rather than offspring produced by each individual and sampled as breeders on the following years (Salles et al., 2020). It is important to note that in long term surveys of a given wild population defined by its geographic perimeter, these two options are often equivalent.
Lack of knowledge on environmental and genetic variation of LRS in wild marine populations is obviously problematic when considering the stakes of their adaptation to global change. This is particularly true for coral reefs that are severely impacted and threatened by climate change (Hughes et al., 2003; van Hooidonk et al., 2016). The degradation and loss of coral reef habitats affects many fish species (Jones et al., 2004). Coral reef fish populations are experiencing a massive shift in selective forces that challenges their connectivity, self-recruitment, and therefore their conservation (Munday et al., 2008; 2009). Previous work in the long-term monitored wild population of the orange clownfish (Amphiprion percula) at Kimbe Island, Papua New Guinea, showed that the habitat of the breeders made the greatest contribution to LRS variation among individuals (Salles et al., 2020). Overall, this work demonstrated that there was little genetic variation in LRS (less than 2%) and therefore little microevolutionary adaptive change (Salles et al., 2020). Finally, the sustainability of this population is therefore mostly under the control of the environment that the local breeders are experiencing.
Most studies investigating the genetic variation of fitness in monitored wild populations found that only a little proportion of this variation was attributable to genetic variation (Hendry et al., 2018, Moiron et al., 2022), although Bonnet et al., (2022) found substantial genetic variation. This is expected from populations at evolutionary equilibrium when no adaptive evolution is ongoing over the rather limited timeframe of the survey as compared to the timeframe of the evolution of species. Ecological mechanisms driving fitness variation are rarely identified. In this study, we disentangle different ecological sources of variation in LRS that were previously combined into a unique notion of habitat (e.g., host anemone species, depth) and considered the role of spatial autocorrelation. These ecological effects may affect the contribution of individuals to recruitment through their LRS in this wild coral reef fish population. Our aim is also to improve the understanding of the capacity of wild clownfish to maintain self-recruitment above 50% in the population (Salles et al., 2016). This is crucial because the habitat of the clownfish is changing and self-recruitment has been shown as a key contribution for the population maintenance in the long term (Jones et al., 1999). In this system, most of the variation in LRS was explained by a breeder’s habitat and the lagoon where they live (Salles et al., 2020). In the worldwide context of natural habitat degradation that also affects host anemones, it is crucial to separate the role of specific interactions from other ecological factors affecting clownfish survival. The relative contribution of the anemone species, its depth and the local density in anemones to the LRS of clownfish remains unknown to date. Knowledge of this relationship would improve our understanding of clownfish ecology and adaptation to their complex habitat. Moreover, the spatial distribution of anemones where the fish live is not homogeneous around Kimbe Island. Because of this spatial heterogeneity, anemones are not spatially independent so that disentangling the effects of ecological conditions (density, depth, etc.) would be inappropriate if this spatial autocorrelation is unaccounted for.
Spatial autocorrelation of LRS is found if similarity in LRS exists amongst geographically close individuals (Legendre, 1993). Although studies on this topic remain rare, they have mainly documented spatial variation in fitness proxies (Kalisz, 1986; Stratton & Bennington, 1998; Wilkin et al., 2009). For example, Wilkin et al., (2009) found spatial variation for the average mass of great tit fledglings – which is a fitness proxy – in association with soil concentration in calcium in a woodland. Traditionally used in macro-ecological or population genetic approaches, spatial autocorrelation recently gained interest for studies of fitness variation (Stopher et al., 2012; Marrot et al., 2015). This is because spatial autocorrelation causes pseudoreplication if not taken into account, which might strongly bias statistical estimates (Haining, 2003). A meta-analysis conducted on 24 studies using linear regressions found that on average, model coefficients were biased by 25% when spatial autocorrelation was unaccounted for (Dormann, 2007). For instance, heritability and selection can be overestimated if spatial dependency among sampling units is ignored (Stopher et al., 2012; Marrot et al., 2015; Gervais et al., 2022). Around Kimbe Island, similar complications may arise from the spatial distribution of anemones associated with heterogeneity in depth and density by causing spatial dependency among clownfish. This study has two aims: (i) to build a spatially explicit model estimating the effect of environmental features on the LRS of clownfish, and (ii) to compare this spatially explicit model to a spatially naïve model to quantify the potential bias induced by spatial dependency. To accomplish these aims, we used a geostatistical method that takes into account spatial autocorrelation at multiple scales to disentangle the relative contribution of environmental effects independently from spatial effects. Overall, we expected the effect of the host anemone species, local density, and depth to be overestimated when spatial autocorrelation is not taken into account.
Material and methods
Study system
A natural population of orange clownfish (A. percula) living on the reef surrounding Kimbe Island (Fig. 1a; 5°12’22.56’’ S, 150°22’35.58’’ E), West New Britain Province of Papua New Guinea, was surveyed every second year from 2003 to 2013. A. percula lives in a mutualistic association with one of two host sea anemone species, Heteractis magnifica and Stichodactyla gigantea. All anemones occupied by a clownfish group have been geographically located, tagged for long term monitoring and sampled every second years for a total of 310 anemones (176 H. magnifica and 134 S. gigantea). A size-based dominance hierarchy exists within each clownfish group of typically three to five individuals: the female is the largest, the male is the second largest, and the immatures rank progressively lower in the hierarchy as they decrease in size. If the single female adult of a group dies, then the male changes sex to female, and the largest immature from the anemone moves up in the hierarchy to a sexually mature male. Reproduction occurs year-round, with females laying a few hundreds of eggs in a clutch near the pedal disk of the host anemone each lunar month. The eggs hatch after c.a. 7 days of paternal care into larvae that spend c.a. 10 days in the pelagic environment (Roux et al., 2020) before settling into an anemone, either at their natal location (Kimbe Island) or elsewhere along a declining dispersal kernel (Planes et al., 2009). On the reefs around Kimbe Island, anemone positions were recorded with a GPS, and depth was measured using a dive computer. Anemones were surveyed during each biannual sampling period, and all clownfish individuals within the anemones were counted, captured using hand nets, measured using callipers (total length, TL), and fin-clipped underwater for genetic analysis. Each fish was processed on SCUBA in situ, with each individual caught, processed next to the anemone, and then released back on the same anemone. The biggest fish in each anemone was identified as the female, the second largest individual was assumed to be the male, and all other individuals were recorded as subadults (TL, >25 mm) or new recruits (TL, <25 mm). Every sampled fish was genotyped at 22 polymorphic microsatellite loci (Bonin et al., 2016) and the pedigree of the population has been constructed (see Salles et al., 2016 and Salles et al., 2020 for a complete description of the pedigree construction of this population). The pedigree revealed half of the juveniles successfully recruiting were progeny of local breeding pairs (Salles et al., 2016).
Figure 1 - Spatial distribution of the mean LRS. Spatial distribution of the mean LRS, averaged among all breeding fishes living on the same anemone (A), depth (B) and anemone species (C). Each point represents an anemone for a total of 302 anemones sampled at Kimbe Island.
LRS & Ecological drivers
The high level of philopatry in this wild population of fish represents a rare opportunity to explore variation in LRS in a marine fish with a pelagic larval stage, given the fact that most estimations of LRS in the wild have been conducted in birds and mammals. Based on the pedigree, the LRS of an individual has been defined as the total number of descendants produced during its lifetime on a biannual basis that successfully recruited at Kimbe Island. Note that, contrary to most studies using LRS, our LRS measure does not refer to the number of offspring produced by each individual and sampled as breeders on the following years. A recruit is defined as an adult (or subadult) sampled and born on the reef surrounding Kimbe Island. This choice was made to adapt LRS to the group hierarchy of clownfish breeders and to be a sensible measurement of the relative ability of clownfish individuals to participate to the population self-recruitment. Our focus was made on the habitat of breeders and not on the habitat conditions where the larvae settle, which could have an effect on the LRS of their parents. The focus on the habitat of breeders rather than larvae should however not affect our results because clownfish offspring recruited in similar proportions in the two species of anemones around Kimbe Island (Salles et al., 2016). Furthermore, settlement habitat conditions were partially considered in the local density in anemone species within a 200m-radius around the focal anemone we investigated.
We chose to disentangle the effect of four ecological drivers: the anemone species where the parents are living (hereafter called the “focal anemone”), the density in H. magnifica in a 200m-radius around the focal anemone, the density in S. gigantea in a 200m-radius around the focal anemone and the depth. The micro-habitat provided by the two anemone species is remarkably different in terms of shape, size and toxicity (Dunn, 1981). Specifically, H. magnifica has larger and less toxic tentacles than S. gigantea (Nedosyko et al., 2014). We chose to measure the density of anemones in an area with a 200m-radius around the focal anemone because the spatial distribution of recruits is not heterogeneous around their parents in this population, and the mean distance between self-recruits and their parents was equal to 254m (± 189m), with 47% of recruits located at less than 200m from their parents (Appendix 1). This mean distance between recruits and their parents was significantly different (Wilcox test P-value < 0.001) from the mean distance between each anemone occupied by a self-recruit and the other anemones in this population (388m ± 37m), suggesting a potential effect of the local environment around the focal anemone on the probability of self-recruiting (Jones et al., 2005; Buston et al., 2011). We chose to investigate depth as a potential driver because the two anemones species are not equally distributed across depth (Fig. 1B, 1C) with H. magnifica more likely found in deeper lagoon waters (4.10m ± 2.61) than S. gigantea, which is found closer to shore (1.34m ± 0.51). Moreover, because larvae may recruit on deeper anemones than their anemones of birth (Buston et al., 2011), we expected that the probability to produce self-recruits increases as the depth of the anemone hosting the parents decreases, independently from the species of the anemone.
Statistical analysis
All statistical analyses were carried out using the R software (v. 4.0.3) (R Core Team, 2020). In this population, the mean of individual LRS was 1.27 (± 2.32) and reached a maximum of 20 self-recruits. Because the statistical distribution of individual LRS included an excess of zero (58% of individuals did not produce self-recruits in the population), we chose to test for the effect of ecological drivers by using a zero-inflated Poisson model. In this type of model, the excess of zero is analysed by a separate process (the probability that fish will not produce self-recruits, because in a zero-inflated model, it is the probability to have a zero that is estimated, generally with a logit link function) from the count values (the number of self-recruits produced by breeders, generally with a log link function). This choice was corroborated by a smaller Akaike information criterion (AIC, Akaike, 1998) of a null zero-inflated model (AIC=4246) than the AIC of a null generalized linear model following a Poisson distribution (AIC=5423). We included in the model the first year of sampling of each fish to take into account the effect of lifespan on LRS. Although this variable is an incomplete proxy of lifespan, it is equally spatially distributed in the population, preventing it from affecting the spatial distribution of LRS.
The explanatory variables included in the model were the species of the parental anemone (S. gigantea vs. H. magnifica), the depth and the density of the two species of anemone. This model is not spatially structured and we compared the estimates of responses variables to a spatially explicit model (see following section).
Geostatistics
The visualization of the spatial distribution of the LRS suggested that this variable is not homogeneous across space (Fig. 1A). Before conducting geostatistical analyses, we formally tested if LRS was significantly spatially autocorrelated around Kimbe Island by using spatial correlograms. Our analysis revealed a significant and positive spatial autocorrelation at 230m and 460m and a significant negative spatial autocorrelation at 690m and 920m (Fig. 2A). This result indicated that LRS is not spatially independent and geostatistical modelling is required to disentangle the effect of the ecological drivers. Based on the abundant literature that compare the performance of different geostatistical methods (Dormann et al., 2007; Kissling & Carl, 2008; Beguería & Pueyo, 2009; Bini et al., 2009; Diniz-Filho et al., 2009; Beale et al., 2010; Marrot et al., 2015), we chose to use the Principal Coordinates of Neighbourhood Matrix (PCNM) approach (Dray et al., 2006) because it is the best method to remove spatial autocorrelation, independently from its direction and its strength (relaxing the assumption of isotropy and stationary). The general principle of the PCNM is based on the extraction of eigenvectors from a truncated distance matrix among spatial units. Each of these eigenvectors (called after “Egs”) describes a specific spatial scale (regarding the spatial distribution of sampling units) and can be subsequently added as a covariate in any statistical model to model the spatial structure of the response variable (Borcard & Legendre, 2002). The use of the PCNM approach can be decomposed in four steps: (i) compute a truncated pairwise Euclidean distance matrix between spatial units, (ii) perform principal coordinate analysis of the distance matrix and extract the corresponding eigenvectors (Egs) for a total number of Egs equal to the number of sampling units (302 clownfish with non-null LRS at Kimbe Island), (iii) select a subset of Egs to include as covariates in any model. Although the whole set of Egs describes the whole spatial structure of the dataset, it is not possible to include all Egs in any model without over-parameterizing it. This Egs selection part of the analysis is crucial and strategies to select Egs are diversified (Blanchet et al., 2008). We chose to conduct a backward-forward selection based on AIC (Ficetola & Padoa-schioppa, 2009) to select the best set of Egs that explains spatial variation in LRS. This selection led us to keep 29 Egs. This set of Egs represented the minimum number of Egs necessary to describe spatial variation in LRS. (iv) Include the set of selected Egs in any model to take into account spatial autocorrelation between units of sampling. We included the set of 29 Egs as covariates in our model, and we compared the sign and the absolute value of ecological variables with a “naïve” model (not including the Egs). This approach allowed us to quantify the biased induced by spatial autocorrelation when it was ignored. In both models, all variables (ecological variables and Egs) have been standardized (o mean, unit variance).
Figure 2 - Moran’s I correlogram for LRS at Kimbe Island. Moran’s I with associated standard deviations before (A) and after (B) correcting for spatial autocorrelation with the Egs selected by the PCNM (see method). For B, Moran’s indexes have been estimated on the residuals from a model linking the response variable (LRS) and the Egs only. Moran’s index was calculated at 4 distance classes, from 230m to 920m, with a distance lag of 230m. There is a significant and positive spatial autocorrelation at 230m and 460m and a significant negative spatial autocorrelation at 690m and 920m in LRS before correcting for spatial autocorrelation (A). Spatial autocorrelation was undetected in LRS after correcting with the PCs (B).
Results
Strong spatial autocorrelation
Spatial autocorrelation detected in LRS (Fig. 2A) was totally removed by the set of 29 Egs selected by the PCNM method (Fig. 2B). We added this set of Egs as covariate in the zero-inflated model and we compared it with a naïve model. The estimates and their associated standard errors and P-values are summarised for the non-spatial and the spatially explicit model in Table 1. The inclusion of Egs altered the magnitude, the sign and the standard errors of estimates both in terms of probability of not self-recruiting and the count number of recruits produced (Table 1). Detailed results for each parameter are described below.
Table 1 - Zero-inflated non-spatial vs. spatially explicit model estimates associated with their standard errors. The probability to not self-recruit is represented by the logarithms of the odds ratios. Species: S. gigantea represents the effect of living on S. gigantea compared to H. magnifica. D. 200m S. gigantea and D. 200m H. magnifica represent respectively the effect of the density in S. gigantea and H. magnifica in a 200m-radius around the focal anemone. In bold the variables associated with P-values < 0.05.
Probability to not self-recruit | Non-spatial model | Spatial model (including Egs) | ||||
Estimates | Standard error | P-value | Estimates | Standard error | P-value | |
Intercept | 0.381 | 0.010 | <0.001 | 0.281 | 0.121 | 0.020 |
First year sampling | 0.431 | 0.073 | <0.001 | 0.475 | 0.081 | <0.001 |
Depth | 0.026 | 0.082 | 0.750 | -0.055 | 0.120 | 0.644 |
Species: S. gigantea | -0.527 | 0.174 | <0.01 | -0.693 | 0.248 | <0.01 |
D. 200m S. gigantea | 0.060 | 0.070 | 0.400 | 0.074 | 0.095 | 0.434 |
D. 200m H. magnifica | -0.143 | 0.070 | 0.040 | 0.002 | 0.091 | 0.981 |
Number of recruits produced | ||||||
| Intercept | 0.732 | 0.051 | <0.001 | 0.830 | 0.056 | <0.001 |
| First year sampling | -0.312 | 0.039 | <0.001 | -0.246 | 0.041 | <0.001 |
| Depth | 0.029 | 0.040 | 0.471 | 0.001 | 0.050 | 0.989 |
| Species: S. gigantea | 0.383 | 0.075 | <0.001 | -0.047 | 0.101 | 0.642 |
| D. 200m S. gigantea | 0.125 | 0.031 | <0.001 | 0.084 | 0.041 | 0.038 |
| D. 200m H. magnifica | 0.126 | 0.031 | <0.001 | 0.071 | 0.040 | 0.075 |
Degree of fit (including all variables) | ||||||
| Adjusted R² | 0.16 | 0.34 | ||||
AIC | 4042 | 3859 | ||||
Non spatial and spatially explicit model comparison
The comparison between the non-spatial full model and a spatially explicit the model (including all variables and the set of 29 Egs) showed that the inclusion of Egs increased the degree of fit of the model with an improvement of adjusted R² of 0.18 (the non-spatial and the spatially explicit model explained respectively 16% and 34% of the variance in LRS) and a drop of AIC of 183. This AIC drop means that the inclusion of the 29 Egs significantly improved the degree of fit of the model without over-parameterizing the model.
Year of sampling effect
In both non-spatial and spatially explicit models, the year of first sampling (first year sampling) increased the probability to not self-recruit and in the presence of self-recruitment decreased the number of recruits with the same magnitude (Table 1). This means that earliest sampled fish had the highest probability to produce at least one self-recruit and produced also more self-recruits than late sampled fish because mathematically, the monitoring period is more likely to overlap with the full lifespan of fish sampled earlier, which increases the probability to produce a self-recruit that we had the opportunity to sample.
Ecological effect on LRS components
The probability that fish will not produce a self-recruit (represented by the logarithm of the odds ratio; which estimate and 95% confidence interval are presented between brackets) was significantly affected by the focal anemone species (-0.693 [Standard Error: 0.248]) for breeders living on S. gigantea) in the spatially explicit model. Interestingly, this was also the case in the non-spatial model where the probability that fish will not produce a self-recruit was significantly affected by the focal anemone species (-0.527 [SE: 0.174]) for breeders living on S. gigantea, and the 200m-radius density in H. magnifica (-0.143 [SE: 0.070]).
We back transformed the results of the spatially explicit model from the latent (inferred through the model link function) to the observed scale for the sake of clarity. Fish living on S. gigantea had odds to not produce a self-recruit lowered by half (exp(-0.693) = 0.553), and therefore twice the odds to produce a self-recruit, in comparison to those living on H. magnifica. This means that, while our model estimated a probability to not self-recruiting of 57% for all fishes, this probability was lowered by 17% for fishes living on S. gigantea.
The other component of LRS, which was the number of recruits produced by breeders, was not affected by the anemone species and was positively (but marginally significantly) affected by the density of S. gigantea in a 200m radius around the parental anemone (0.084 [SE: 0.041]) in the spatially explicit model. It is interesting to note that it was significantly affected by the focal anemone species (0.383 [SE: 0.075]) for breeders living on S. gigantea) and the 200m-radius density in both anemone species (0.125 [SE: 0.031]) and 0.126 [SE: 0.031] for the density in S. gigantea and H. magnifica respectively) in the non-spatial model, again emphasizing some effect of spatial auto-correlation that is not taken into account in traditional analyses.
Discussion
Our spatially explicit analysis of long-term monitoring data from a wild clownfish population in which pedigree and individual geographic location are known revealed two main results: (i) the spatially explicit model outperformed the non-spatially-explicit model, so that ignoring spatial autocorrelation affected the estimates of our zero-inflated model, (ii) breeders living on S. gigantea contributed to a larger extent to the self-recruitment of the population than those living on H. magnifica, and this effect was independent from the local density, the depth and the spatial structure.
Why use a spatially explicit approach?
Beyond taking into account the spatial structure of clownfish LRS to identify its ecological drivers independently from spatial effects, our study aimed at evaluating the consequences of not accounting for spatial autocorrelation. We successfully removed spatial dependency at all scales in LRS by using the PCNM approach. This method presents numerous advantages: it allows to control for spatial autocorrelation at any scale, it does not need the assumption of stationary (spatial autocorrelation does not depend on direction) and isotropy (spatial autocorrelation is constant across the study site), and it can be incorporated with any statistical models commonly used by ecologists (Borcard & Legendre, 2002). The inclusion of Egs doubled the variance explained in the LRS at Kimbe Island (the non-spatial and the spatially explicit model explained respectively 16% and 34% of the variance in LRS, Table 1), suggesting that unaccounted spatially autocorrelated variables affect the LRS. Note however that although the variance explained by a model including only the ecological variables explained 16% of the variation in LRS, it is incorrect to attribute 18% of variation explained by the Egs because ecological variables and Egs explained shared variation in LRS. Moreover, a model including only the set of Egs explained 27% of the variation in LRS (not shown in Table 1).
The comparison of the non-spatial and the spatially explicit models including the Egs revealed that, when unaccounted for, spatial autocorrelation biased the absolute values of estimates and their confidence intervals. This result was not surprising and confirmed a large body of evidence, both on simulated and real data, showing the effect of spatial autocorrelation on estimates (e.g. Lennon, 2000; Dormann, 2007; Beale et al., 2010; Le Rest et al., 2013). As pseudoreplication, spatial autocorrelation decreases the precision (increasing the type 1 error) around estimates and alters (in a non-predictable way) their absolute values when unaccounted for if the response and the explicative variables are spatially autocorrelated (Beale et al., 2010). Among the ecological variables, three of them had their confidence intervals overlapping 0 (for both LRS components) when spatial autocorrelation was taken into account. Beyond statistical considerations, these results highlight how spatial autocorrelation may result in biological misinterpretations when unaccounted for. For example, while a non-spatial analysis would conclude that wild clownfish surrounded by a high density in anemone (within a 200m-radius) produce more recruits than those living in empty areas, our spatially explicit model did not detect any effect of the density in anemones.
Spatial autocorrelation of LRS
Our results showed that there was a spatial structure of self-recruitment around Kimbe Island. This is in line with a previous study conducted in the same clownfish population where it was suggested that the location of the lagoon could influence LRS (Salles et al., 2020). Our approach allowed us to go further and detected a strong spatial autocorrelation of LRS. As mentioned earlier, it would be inappropriate to try and to provide a precise estimate of its strength. This effect occurred over a relatively small spatial scale (less than 1km). This result is original because spatial autocorrelation is a property of ecosystems and can be found at any scale, from micrometres to the continental scale (Legendre, 1993), but it was mainly investigated by ecologists interested in macro-ecological patterns at very large scales (Dormann, 2007) and population geneticists estimating the limits of dispersal within and between populations (Sokal et al., 1989; Vekemans & Hardy, 2004). It was in fact rarely investigated in studies exploring individual fitness variation. Only a few examples from wild plant and animal populations indicate that spatial autocorrelation can strongly impact the estimation of natural selection (Marrot et al., 2015; 2022), estimates of genetic variation for fitness-related traits (Stopher et al., 2012; Gervais et al., 2022) and now ecological drivers of LRS with this study. The biological mechanisms that underlie this pattern of spatial autocorrelation remain unknown to date. Yet, identifying them would improve our understanding of clownfish LRS and self-recruitment.
One limitation of our spatial analysis comes from the Euclidean distances that we used. Euclidean distances likely may not accurately represent connectivity in the ocean, where islands and currents can impede routes of travelling larvae (Treml et al., 2008). Oceanic distances accounting for actual currents around the island would be a good alternative to Euclidean distances (White et al., 2010; Truelove et al., 2017) because currents are generally predicted to have a strong effect on dispersal pattern (Cowen et al., 2000; but see Van Wynsberge et al., 2017). They would likely help to explain the variation of LRS more accurately. Our results have to be interpreted with caution since currents and the role of the island of Kimbe were not considered in our spatially-explicit analysis of LRS. One should note that a previous study showed that this currents were likely not influencing self-recruitment at short distances of about one kilometer around Kimbe Island (Buston et al., 2011) although at the broader scale of Kimbe Bay, a direct relationship was found between geographic distance and genetic differentiation (Pinsky et al., 2017). Furthermore, although not explicitly modelled, the potential role of the island was at least partially taken into account by the PCNM method because it estimates patterns of spatial autocorrelation at any scale (Borcard & Legendre, 2002).
The host anemone species affects clownfish LRS
Our spatially explicit approach allowed us to identify some ecological drivers of wild A. percula clownfish LRS without potentially confounding their effect with spatial effects. Our results confirmed the findings from a previous study (Salles et al., 2020) conducted on the same clownfish population: an important part of variation in components of fitness linked to self-recruitment was explained by the ecological conditions of the habitat of breeders. Our approach allowed us to go further and disentangle the relative contribution of different ecological factors. Among them, we demonstrated that the anemone species is the only factor affecting the LRS. Specifically, we found a higher probability to produce at least one recruit for breeders living in S. gigantea than those living in H. magnifica. The spatially explicit approach showed that this effect was likely independent from the depth, the local density in anemones around the focal anemone, and the spatial structure of the population, suggesting a direct causal relationship between the anemone species and LRS.
Experiments will be necessary to identify the causal mechanism linking the anemone species to the variation in LRS that we observed in wild populations. Several hypotheses can be assessed to explain our finding. First, one may push forward the lower clownfish fecundity caused by hormonal stress resulting from host anemone bleaching (Beldade et al., 2017) as a potential explanation of our results if anemone bleaching affected differentially the two anemone species. However, anemone bleaching cannot explain our results because evidence suggests that bleaching affects all anemone species with the same magnitude (Hobbs et al., 2013) and no anemone bleaching event was reported over the duration of the survey. Second, one candidate hypothesis to explain the higher LRS of clownfish living on S. gigantea is that this anemone species is more toxic than H. magnifica (Nedosyko et al., 2014), which could offer more protection to clownfish, and ultimately affect positively its fitness because of a higher survival. Third, along with this difference in protection allocated to clownfish, female size may vary (Buston, 2004) as it is mediated by the group size which has been shown to be driven by the size and shape of the anemone (Chausson et al., 2018). A direct link can be drawn to a recent study on A. percula (Barbasch et al., 2020), using both experimental design and field work data, which revealed a strong and positive relationship between anemone area, clownfish size, and reproductive success, driven by food availability. Whether S. gigantea provides clownfish with a higher food availability remains however to be tested. Observations in the field nevertheless suggest that clownfish spend more time feeding out and above S. gigantea anemones than clownfish hosted by H. magnifica (S. Planes personal observation).
The host anemone species, whether through their toxicity, size or available food, or other unknown features, directly influences wild clownfish self-recruitment, and thereby their ability to locally replenish the population around Kimbe Island. Any changes to the ecology of the anemone species would therefore impact the sustainability of the clownfish population. The limits of the perspective on clownfish sustainability provided by our finding is restricted to the frontiers of the local population. Broader consequences cannot be revealed by our investigation of the local component of the clownfish lifetime reproductive success. Whether an integrative perspective considering long distance dispersal would reveal a trade-off between a breeder’s ability to produce self-recruits and successful settlers on foreign anemones remains unknown to date.
Conclusion
Ecological conditions experienced by wild clownfish are impacting components of fitness linked to self-recruitment such as the LRS. Our study showed that taking into account spatial autocorrelation is efficient for disentangling ecological sources of variation from confounding spatial effects and prevents biological misinterpretation. Geostatistical models such as PCNM represent a powerful tool to identify the spatial autocorrelation of LRS, even at small spatial scales. In a context of worldwide threat on coral reefs (Hughes et al., 2003; van Hooidonk et al., 2016), it is crucial to investigate the sources of individual variation in LRS to predict the resilience of wild coral reef fish populations. As part of a marine reserve (Almany et al., 2007), the lagoons around Kimbe Island where wild clownfish are monitored are protected from fishery or direct habitat degradation and this clownfish population has proven to persist in absence of migrants (Salles et al., 2015). However, our study revealed that the subsistence of the population could be altered if the proportion of S. gigantea was to change independently from density, depth and other environmental aspects and a similar concern could be raised if the spatial distribution of anemone species was to change.
Acknowledgements
Research visas were approved by the Papua New Guinea (PNG) government and research protocols were endorsed by the Board of Mahonia Na Dari Research and Conservation Centre, Kimbe, PNG. We thank the large number of volunteers who assisted in the field and collected tissue samples. Mahonia and FeBrina provided essential logistic support. We are grateful to the traditional owners of the reefs near Kimbe Island for allowing us access to their reefs.
Preprint version 3 of this article has been peer-reviewed and recommended by Peer Community In Evolutionary Biology (PCI Evol Biol) (https://doi.org/10.24072/pci.evolbiol.100698; de Villemereuil, 2024).
Funding
PM post-doc salary was funded by a grant from the Laboratoire d’Excellence CORAIL. This research was supported by Laboratoire d’Excellence CORAIL, Expenditure Review Committee, Coral Reef Initiatives for the Pacific, the Global Environment Facility Coral Reef, Targeted Research Connectivity Working Group, National Science Foundation, the Australian Research Council Centre of Excellence Coral Reef Studies, The Nature Conservancy, Total Foundation, James Cook University, King Abdullah University of Science and Technology.
Conflict of interest disclosure
The authors declare that they comply with the PCI rule of having no financial conflicts of interest in relation to the content of the article. Benoit Pujol is a recommender for several Peer Communities (PCI Evol Biol, PCI Ecology and PCI Registered Reports).
Data, scripts, code, and supplementary information availability
Data are available online: https://doi.org/10.5281/zenodo.13806778 (Marrot et al., 2024).
Scripts and code are available online: https://doi.org/10.5281/zenodo.13806778 (Marrot et al., 2024).
Supplementary information is available online: https://doi.org/10.5281/zenodo.13806778; (Marrot et al., 2024).
Appendices
Appendix 1 - Statistical distribution of the mean distance between each anemone occupied by a self-recruit and the other anemones (A) and the distance between self-recruits and their parents (B) at Kimbe Island. Red dotted lines represent the mean of each statistical distribution (A: 388m ± 37m; B: 254m ± 189m). Based on a Wilcoxon test, these two statistical distributions were significantly different (P < 0.001). Note that 47% of self-recruits were located at less than 200m from their parents.
