The study was undertaken over a period of 1 year across four field sites within Tolo Harbour, Hong Kong: Center Island – CI, Che Lei Pai – CLP, Port Island – PI, and Tung Ping Chau – TPC (Figure 1). These sites represent a gradient in coral cover, CI < CLP < PI < TPC (Duprey et al., 2016; Wong et al., 2017) and several environmental parameters. Center Island is the site closest to the harbor, with almost no hard corals, and is also the most polluted compared to the other sites. Tung Ping Chau, on the other hand is a marine reserve with over 60 species of hard corals. The study used a simple toolkit to allow the quantification of key ecosystem functions and overall multifunctionality.

Primary productivity is the rate at which energy is stored as organic matter (Fahey and Knapp, 2007). We measured primary productivity as net mass gain (%) of macroalgae (Ulva) on substrata protected from grazers. We deployed four replicates at each site in plastic bottles with holes, 4 m from the ocean floor (to standardize light availability). Ropes were attached to bricks and deployed on the substratum at 5–10 m intervals (Figure 2). After 48 h, we collected the macroalgae, returned them to the laboratory and assessed their mass loss/gain (ΔM) as follows:

ΔM = (Mf − Mi)/Mi

Where i and f refer to initial and final wet weights as per Rasher et al. (2013).

Herbivores play a vital role in maintaining a healthy coral-dominated community through intense feeding and grazing of unwanted macroalgae that indirectly interfere with the growth, reproduction and survivorship of corals (Burkepile and Hay, 2008). To estimate herbivory, we exposed samples of macroalgae (referred to as “algae pops”) to determine their susceptibility to grazing. The macroalgae were deployed on 90 cm lengths of three-stranded nylon rope, upon which the algae had been grown for 2 weeks prior to use (Figure 2). At each site four replicates were deployed and spaced 5–10 m apart on the sea floor. Following exposure for 48 h, loss of biomass (ΔC) was calculated as follows:

ΔC = {[(1 + ΔM)^∗Mi] − Mf}/Mi

Where i and f refer to initial and final wet weights as per Rasher et al. (2013) and ΔM is the % mass gain calculated from the primary productivity assay.

Fish feeding intensity is often used as a measure of predation in a shallow water ecosystem. A simple assay was used to assess predation – “the squid pop” (Duffy et al., 2015). This refers to a piece of dried squid deployed at the end of a tether tied to a stake as bait for fish. In order to standardize the squid pop assay, we cut the dried squid into disks that all had the same size. Following this, the disks were deployed from the seafloor by tethering them to individual metal stakes using fishing line. Squid pops were deployed at three sites (CI, CLP, PI; 25 replicates per site). Squid pops were not deployed in the fourth site – TPC due to inclement weather forecast and budget restrictions. One hour after deployment, we returned to observe whether or not fish consumed the bait. The degree to which the squid pops were consumed was used as an indicator of the level of predation (or the activity of fish feeding) in the marine ecosystem (Duffy et al., 2015; Figure 2).

Bait loss = 1 − [(Si − Sf)/Si]

Where i and f refer to the initial and final number of squid pops, respectively.

To measure organic matter decomposition and carbon sequestration potential, we used the tea bag assay as per Keuskamp et al. (2013). The tea bag assay was developed to measure organic matter degradability in a standardized, cost-effective manner across diverse ecosystems globally. Recent studies have demonstrated its usefulness and its applicability across aquatic and terrestrial ecosystems (Al-Maliki and Al-Masoudi, 2018; Mueller et al., 2018; Duddigan et al., 2020). Further, Seelen et al. (2019) demonstrated the applicability of the tea bag method to marine environments by incorporating a leaching factor (40% for green tea and 20% for rooibos tea).

We deployed 10 pairs of commercial rooibos tea bags and green tea bags (rooibos tea – ASIN: B00BAUF9KA, Lipton, Unilever, United Kingdom and green tea – B0173LLHNM, Lipton, Unilever, United Kingdom) in nylon mesh (200 μm) that were held within a porous plastic bottle and buried ∼10 cm beneath the sea floor for a period of 6 months in the four study sites (Figure 2). Upon retrieval, the tea bags were oven-dried at 60°C for at least 72 h and weighed after removing the mesh bags. Decomposition rate (k, d^–1) of rooibos tea and organic matter stabilization (S, %) of green tea were calculated using the following equations:

k = −Ln [(final mass − initial mass)/duration of deployment]

s = 1 − [(final mass − initial mass)/initial mass of tea] × 100.

In order to assess the ability of a given site to perform multiple ecosystem functions, we calculated a multifunctionality index following a geometric mean approach and a multiple threshold based approach (Byrnes et al., 2014). For each of the five functions (primary productivity, herbivory, predation, organic matter decomposition and carbon sequestration), we defined “maximum functioning” as the highest value among the recorded observations. Following this we calculated the standardized percent functioning for each function and took the geometric mean of these values to get the multifunctionality index for every site. Next, we selected a multifunctionality threshold (50% of the calculated maximum functioning values) to evaluate the effect of coral species richness on ecosystem function. The rationale for this is that Gamfeldt et al. (2008) made reasonable observations assuming a 50% functioning threshold to study the effects of species loss on multifunctionality. We recorded 4 values between 0 and 1 for each site, where 1 represented a function was maintained over the selected threshold and 0 represented a function that was maintained below this threshold. We took the sum of these values for each site, which resulted in numbers between 0 and 4, with 4 indicating that the site was able to maintain all ecosystem functions above the specified threshold.

Plate photos from 12 autonomous reef monitoring structures (ARMS), that were deployed in the same four study sites were analyzed and annotated by functional group. The ARMS are similar to mini apartment blocks that are made of 9 stacked PVC plates and are deployed on the sea floor. The structures have cave-like spaces for marine fauna to crawl in, hide and settle. Following deployment for 2 years, the ARMS were retrieved, and everything on the plates was collected and identified (Leray and Knowlton, 2014). During this step, high resolution photographs of individual plates were taken (Figure 3). Similar to settlement plates, ARMS plate photos can be used to assess the diversity, coverage and size of sessile taxa. Recently, David et al. (2019) validated ARMS as a promising monitoring tool for hard bottomed communities and to investigate environmental effects. An online software CoralNet, a repository and resource for benthic image analysis was used that helped generate 50 random points on each plate photo (Beijbom et al., 2015). The software uses computer algorithms that allows fully and semi automated annotation. Based on already identified taxa, the organisms on the plates were assigned to various taxonomic groups such as sponges, bryozoans, bivalves, etc. The assignments were made based on annotation categories used by the NOAA CRED program that is already available on CoralNet. In total, 204 plate photos were annotated at the functional level. Relationships between the abundance of various taxonomic groups on ARMS plates and individual ecosystem functions were assessed to understand the mechanisms that drive variability in functionality.

Environmental data was obtained from surface seawater samples collected from the four study sites by the Hong Kong Environmental Protection Department in 2016 (10 replicates per site, every month). The parameters used for the analysis were – turbidity, total phosphorous, total nitrogen, total Kjeldahl nitrogen (TKN), total inorganic nitrogen (TIN), temperature, suspended solids, secchi depth, salinity, pH, orthophosphate, nitrate-nitrogen, ammonia-nitrogen, dissolved oxygen, chlorophyll-a, and coral cover (Supplementary Appendix Table A.1).

Data were screened for normality using a Shapiro–Wilk W test and for homoscedasticity using Levene’s test. One-way ANOVA was used to evaluate the variability of primary productivity (% mass gain), herbivory (% mass consumed), predation (% bait loss), decomposition rate and carbon sequestration by habitat. Post hoc comparisons to test for significant differences between sites were conducted using Student’s t-test and Tukey’s test. Generalized linear regression was performed with second-order AICc (Hurvich and Tsai, 1989) model selection to obtain the model that could best explain the relationship between the parameters. All statistical analyses were performed in JMP 15.0. Visualizations were carried out using Datagraph.

Net mass gained by protected Ulva (primary productivity), net mass loss by exposed Ulva (herbivory), loss of squid pop (bait) to predators (predation intensity), and carbon sequestration varied significantly among the four habitats (see Supplementary Materials 1–5 for one-way ANOVA results); One-way ANOVA; n = 47; F_(3,44) = 8.9 p = 0.0001; One-way ANOVA; n = 47; F_(3,44) = 301.2 p < 0.001; One-way ANOVA; n = 47; F_(3,44) = 21.1 p < 0.001; n = 47; F_(3,44) = 3.4 p = 0.03, respectively; Figure 4. Post hoc Tukey HSD tests revealed significant differences between the farthest sites (Tung Ping Chau and Center Island) for primary productivity (p = 0.0001), herbivory (p = 0.0001), and predation (p = 0.0001); and adjacent sites for carbon sequestration (Che Lei Pai and Port Island p = 0.04), primary productivity (Tung Ping Chau and Che Lei Pai p = 0.0007; Tung Ping Chau and Port Island p = 0.02), herbivory (all site pair comparisons p < 0.001) and predation (Tung Ping Chau and Che Lei Pai p < 0.0001; Center Island and Port Island p = 0.0001; Tung Ping Chau and Port Island p = 0.04).

To identify the key environmental factors driving the variability, step-wise linear regression with model selection using the second-order Akaike Information Criterion (AICc) was performed (see Supplementary Material for detailed step-wise linear regression results with AICc scores: Supplementary Tables 6A–E). Chlorophyll-a and coral cover were significant drivers (p < 0.05) of primary productivity, on the other hand inorganic nitrogen and dissolved oxygen were closely significant (p = 0.05; best-fit model r² = 0.44). Coral cover was also critical and the most significant (p < 0.05) for determining herbivory, however, secchi depth, an indicator of light (p = 0.05) could not be eliminated (best-fit model r² = 0.96). Temperature and coral cover (p < 0.05) were important in predicting predation intensity (best-fit model r² = 0.66). In determining factors that drove the variability of organic matter decomposition, no significant environmental parameters were observed. However, inorganic nitrogen significantly impacted the carbon sequestration potential in the sites (p < 0.05; best-fit model r² = 0.44).

ARMS plate photo analysis has been proven to be a powerful way to compare marine benthic communities (David et al., 2019). We observed 11 broad groups of organisms present on all ARMS deployed. This included several types of algae – chlorophytes, rhodophytes, and phaeophytes, tunicates, bivalves, and tube worms. To identify the key groups of benthic organisms driving the observed ecosystem function variability, step-wise linear regression with model selection using the second-order Akaike Information Criterion (AICc) was performed (see Supplementary Table 7). We observed no significant effect of the ARMS plates community composition on herbivory, primary productivity, carbon sequestration and organic matter decomposition. On the other hand, the relative abundance of green macroalgae (p = 0.00001), red upright macroalgae (p = 0.00251), soft tube worms (p = 1.46e^–5) and tunicates (p = 0.00044) had a significant effect on the intensity of predation (best-fit model r² = 0.95).

Illumina MiSeq analysis yielded an average of 196540 million reads. High quality reads were obtained upon using a Phred score cutoff of 90% on the paired end pre-processed reads (Supplementary Appendix Tables A.4, A.5). Overall, ∼75% of all reads were annotated to local fish species. Samples from Center Island (the innermost site) did not amplify during the second PCR reaction. Therefore, these samples were excluded from library preparation.

eDNA sequences of a total of 91 taxa were obtained, in 44 families and 14 orders. Over 98% were registered under the Hong Kong Register of Marine Species, a comprehensive database for marine species in Hong Kong waters (Astudillo et al., 2018). Commercially valuable fish species were abundant, such as golden threadfin bream Nemipterus virgatus, flathead gray mullet Mugil cephalus, and yellow croaker Larimichthys crocea. Pollution-tolerant species including Mugil cephalus, and Siganus canaliculatus were also obtained from the samples. When compared with other conventional fish surveillance methods such as trawling (Leung, 1997), line fishing and netting (Hksar Agricultural, Fisheries and Conservation Department, 2016), as well as underwater visual census surveys (Sadovy and Cornish, 2000), eDNA captured ∼75% of taxa previously recorded in the same study sites. Moreover, there were 11 species that were only recorded using the eDNA method, but not observed in the conventional methods (Tsang, 2018).

Relative abundance of all fish varied significantly between sites (One-way ANOVA n = 11; F_(3,8) = 103.5 p < 0.0001). Post hoc Tukey HSD revealed significant differences between the farthest sites Che Lei Pai and Tung Ping Chau and also all adjacent site pairs except Port Island and Tung Ping Chau (Figure 5). It is worth noting that OTU relative abundance is only correlated but does not reflect the relative abundance of the fish species under discussion. This is due to limitations with sequencing methods that yields platform specific data (reads), therefore making relative abundance analysis challenging (Harrison et al., 2020). Nevertheless, acknowledging these constraints imposed by compositional data is vital when analyzing sequencing data.

We classified the observed fish species according to their dietary habits to herbivorous, omnivorous, planktivorous and predatory fish using FishBase (Figure 6). Specifically, relative abundance of herbivorous, omnivorous and planktivorous fish varied significantly between sites, however, the relative abundance of predatory fish showed no significant variability.

We performed step-wise linear regression with model selection using the second-order Akaike Information Criterion (AICc), to characterize the relationship between fish abundance, fish diversity, dietary preference and ecosystem function (see Supplementary Tables 8A–E). Mass loss significantly depended on the relative abundance of herbivorous fish (p = 0.00061; best fit model r² = 0.22). Mass gain significantly depended on the relative abundance of the following 12 genera of predatory fish (p = 0.00369; best fit model r² = 0.46): Platycephalus indicus, Protonibea diacanthus, Trichiurus japonicus, Larimichthys polyactis, Polynemus dubius, Anguilla japonica, Nemipterus virgatus, Collichthys lucidus, Silurus meridionalis, Istigobius campbelli, Cephalopholis boenak, and Larimichthys crocea. The relative abundance of all planktivorous fish (p = 0.01198), omnivorous fish (p = 0.00179), and certain genera of predatory fish (A. japonica, I. campbelli, N. virgatus, C. lucidus, L. crocea, P. indicus, P. dubius, T. japonicus, S. meridionalis, C. boenak, L. polyactis, P. diacanthus; p = 0.00431) were significant drivers of bait loss variability (best fit model r² = 0.52). The relative abundance of 12 genera of predatory fish – C. boenak, C. lucidus, I. campbelli, L. polyactis, N. virgatus, P. indicus, P. dubius, S. meridionalis, T. japonicus, A. japonica, L. crocea, and P. diacanthus significantly drove the variability of organic matter decomposition (p = 0.00126; best fit model r² = 0.24). Finally, N. virgatus, P. indicus, A. japonica, L. polyactis, P. diacanthus, C. boenak, L. crocea, T. japonicus, P. dubius, S. meridionalis, C. lucidus, and I. campbelli significantly determined carbon sequestration variability (p = 0.00224; best fit model r² = 0.18).

In this study, overall multifunctionality was calculated using two approaches – average approach and threshold approach. Multifunctionality using the average approach was calculated without bait loss, as there was no data for Tung Ping Chau. Overall multifunctionality index was 0.70 for Center Island, 0.72 for Tung Ping Chau, 0.92 for Che Lei Pai and 0.51 for Port Island. The number of functions that could be maintained above a 50% threshold was 4, 5, 3, 4, respectively for Center Island, Che Lei Pai, Port Island and Tung Ping Chau (Table 1).

Coral reefs and coral communities are vital components of the Earth’s biosphere dominating several coastal habitats (Knowlton et al., 2010). In recent years, however, these habitats have been destroyed and homogenized owing to unsustainable coastal development, increased levels of pollution, overfishing and disease. This decline in coral cover is projected to continue even if the goals of the December 2015 Paris Agreement are implemented (Hoegh-Guldberg et al., 2019). Changes in habitat (coral reef complexity, coral cover, coral species richness, etc.) can have important implications for the functions that coral reefs provide such as productivity, sediment generation and so on (Perry et al., 2018). Our study sites in Hong Kong are along a gradient of coral cover (Center Island has almost no hard coral species to Tung Ping Chau which has ∼65 hard coral species) and serve as natural laboratories to discuss how changes to coral cover can affect multifunctionality in an ecosystem. Specifically, in testing for the effects of habitat on ecosystem function, we observed a positive relationship between coral cover and primary productivity. This is consistent with several literature showing that homogenized reef communities can jeopardize ecosystem functioning and productivity (Alvarez-Filip et al., 2015; Hughes et al., 2018; Richardson et al., 2018).

Functional diversity has been proven to be an indicator of ecosystem function (Petchey et al., 2004) and has been correlated to productivity (Cadotte et al., 2009). When coral communities shift, due to stressors, some visible changes often occur. One such example is the shift from a coral-dominated community to an algae-dominated one (Brooker et al., 2016). We observed this visually during scuba dives in Center Island and validated it through our analysis on the relative abundance of food sources as captured by the ARMS plates. Center Island recorded high relative abundance of algae, but no hard-coral species. On the other hand, Tung Ping Chau recorded relatively high abundance of both corals and algae (Figure 7).

When looking for functional diversity drivers of overall multifunctionality, the relative abundance of benthic filter-feeders emerged significant, outweighing effects from abiotic factors. This is consistent with literature citing bivalves as pivotal players in ecosystem functioning owing to their contribution to nutrient regeneration and productivity (Norkko et al., 2013).

Nitrogen is the biological nutrient limiting factor in marine ecosystems, as it is generally less available for ocean plants and animals. In this study, there were significant inorganic nitrogen mediated effects on primary productivity and carbon sequestration in the study sites. Consistent with Zhang et al. (2015) and Mueller et al. (2018) increased inorganic nitrogen significantly reduced organic carbon preservation in coastal marine sediments. The results were also was consistent with other studies such as Bristow et al. (2017) which observed that nitrogen exerts a critical control on primary productivity. On the other hand, we observed the opposite relationship between total Kjeldahl nitrogen (TKN: a measure of organic nitrogen) and carbon sequestration, where carbon storage decreased with decreasing TKN concentration. TKN is directly related to the source of bottom deposition. Wastewater discharges often contain relatively high concentrations of organic nitrogen and Hong Kong discharges over 3 million cubic meters of treated wastewater into its surrounding marine environment every day. Archana et al. (2016) recorded that wastewater was indeed the dominant source of the nitrogen pool using stable isotopes. Therefore, characterizing the effect of bottom sediment TKN is valuable in this highly urbanized coastal system.

We observed no nutrient-driven effects on primary productivity, herbivory, predation and organic matter decomposition. This is consistent with Alsterberg et al. (2014), which summarized findings from some studies that demonstrated no effects of nutrients on ecosystem multifunctionality in a marine habitat. We hypothesize that nutrient-driven effects were masked by effects from other significant abiotic factors such as secchi depth (an indicator of light), chlorophyll-a, dissolved oxygen and temperature.

The complexity of the coral-algae-herbivore relationship is well known (Holbrook et al., 2016). A popular hypothesis (that is also somewhat controversial) is that herbivore fishes help corals thrive in a reef benthos by managing the distribution and abundance of algae (Hixon, 2015). However, overfishing of herbivorous fishes can degrade this association (Heenan et al., 2016). The data recorded from the herbivory assay showed significant variability between the four sites. We observed relatively more algal consumption in Tung Ping Chau and Che Lei Pai when compared to what was observed in Center Island and Port Island. We hypothesize that this could be due to several reasons – firstly, we propose that the lack of consumption in Center Island and Port Island could be due to the availability and preference for other food sources, such as Sargassum, which was also present at some sites. This was confirmed by the ARMS plate photo analysis, which revealed several types of algae (Figure 7) in varying abundances at all the four study sites. Secondly, we observed herbivorous rabbitfish and pufferfish at the sites that recorded more herbivory. This was corroborated by fish abundance from eDNA analysis (Figure 6), which suggested the presence of herbivorous and omnivorous fish in increasing abundances from Che Lei Pai to Tung Ping Chau. However, it was unusual that we observed the presence of herbivorous and omnivorous fish in the eDNA analysis at Port Island, yet did not observe any algal mass loss in the herbivory assay. As the coral-algal-herbivore interaction does not occur in isolation of the rest of the ecosystem, we suspect the variability in our dataset to be due to a combination of abiotic factors, other organisms and excessive fishing.

There continues to be some controversy on whether algal distribution and growth in the marine benthos is controlled by top-down herbivory/predation or bottom-up nutrient cycling (Lapointe et al., 2004; Poore et al., 2012). While bottom-up control by resources is relatively well understood, it is necessary to characterize the geography of top-down control by predators and herbivores through space and time (Meyer et al., 2016). One way to do this is to measure prey vulnerability – however, across geographies and larger scales, it is desirable to have a standardized food type and a simpler tool to measure the relative feeding intensity of generalist predators. In line with this, we utilized the squid pop technique as per Duffy et al. (2015). Feeding intensity (bait loss from squid pops) was higher in the sites with higher hard coral species richness. A recent study that employed the squid pop technique found increased feeding intensity to be directly correlated with fish abundance. This pattern could therefore reflect a higher level of predation in Port Island, a site where visual observations have also recorded more fish compared to inshore sites. Consistently, fish abundance from the same sites revealed a similar pattern, thereby corroborating our findings. However, we did not obtain predation intensity for Tung Ping Chau. Therefore, we did not obtain a holistic picture for the difference across habitats in relative predation intensity. We believe these limitations are compensated by the assay’s ease of applicability and standardization, making it possible to replicate the experiment and obtain data fairly quickly.

To determine the overall ecosystem multifunctionality, we followed two standard approaches – the average approach and the threshold approach. Both these approaches include all available measures of ecosystem functions in a given habitat. Consistently across both approaches, Che Lei Pai recorded a relatively higher multifunctionality index and multifunctionality threshold while Port Island recorded a relatively lower index and threshold compared to the other sites. The site with the least disturbance (Tung Ping Chau) did not record the highest multifunctionality index or threshold, neither did the site with the most disturbance (Center Island) record the lowest multifunctionality index or threshold. The site with intermediate disturbance (Che Lei Pai) recorded the highest multifunctionality index and threshold. This was contrary to our hypothesis (multifunctionality increases from inshore – more human impacts to open ocean - less human impacts). The intermediate disturbance hypothesis was presented originally to describe species-mediated effects on corals and trees (Connell, 1978). It predicts that sites with maximum disturbance and least disturbance will have low diversity and that intermediate disturbance will maximize diversity. In line with this, we propose that intermediate disturbance maximizes diversity and correspondingly multifunctionality in Che Lei Pai when compared to the other sites. However, it must be noted that both approaches to characterize ecosystem multifunctionality are not without limitations and multivariate methods for quantification is recommended.

In testing for drivers of multifunctionality, we observed abiotic-mediated (turbidity) and functional diversity-mediated (benthic filter feeder abundance) effects. However, maintaining high levels of multifunctionality in an ecosystem is not likely to be driven solely by these two factors. The maximization of functions favored by the habitat is likely to be largely responsible as per Hensel and Silliman (2013). The findings from our study revealed that when human impact decreases, the consequences to overall ecosystem functioning depend not only on which single functions are taken into consideration. This implies that ecosystem multifunctionality is also driven by the presence or absence of other metazoans and microbial organisms present that drive trophic level interactions and functions such as herbivory and predation. This is consistent with Gamfeldt et al. (2008) who proposed that even in eutrophication impacted coastal marine environments, overall multifunctionality is more susceptible to species loss than are individual ecosystem functions. So, taking these factors into consideration and shifting focus to a variety of functions provided by a diversity of species is likely to provide a holistic picture of ecosystem multifunctionality and will subsequently help in conservation management.