Multifunctionality of an Urbanized Coastal Marine Ecosystem

Coastal marine ecosystems provide numerous invaluable services and benefits to humankind. However, urbanization of coastal areas has homogenized and reduced the biodiversity of the surrounding marine environment and the sustainability of the multiple ecosystem services it provides. Studies have focused on single ecosystem functions despite human populations relying on several functions being delivered at once (known as multifunctionality). This study investigates five ecosystem functions (primary productivity, herbivory, predation, organic matter decomposition and carbon sequestration) and overall multifunctionality in four sites along a gradient of 16 environmental parameters. Ecosystem function varied significantly between sites that were farthest apart. In determining factors that drove ecosystem functioning, we found a positive relationship between coral cover and primary productivity but negative relationships between coral cover and levels of herbivory and predation intensity. Higher temperatures and greater concentrations of chlorophyll-a had a positive impact on predation and herbivory, respectively. Notably, we found a significant negative impact of total inorganic nitrogen and significant positive impact of total Kjeldahl nitrogen on carbon sequestration. Further, individual functions were compared with fish abundance (obtained from seawater eDNA), and benthic community composition (obtained from plate % coverage of autonomous reef monitoring structures). Increasing herbivorous fish abundance had a positive impact on Ulva mass loss. Overall, relative abundance of predatory, omnivorous and planktivorous fish exerted overriding influences on primary productivity and predation intensity, implying that fishing pressure and marine protected area status are important factors. Importantly, we found significant effects from environmental parameters indicating that reliably predicting the effects of future anthropogenic impacts will not be straightforward as multiple drivers are likely to have complex effects. Taken together, urbanized coastal ecosystems exhibit varying levels of multifunctionality depending on the extent of human impact, and the functional diversity of the benthic community present.


INTRODUCTION
Marine environments are significant contributors to both local and global economies, providing ecosystem services and jobs to local people on the one hand, and representing a great bank of undiscovered species, including many of potential importance for human health, on the other. However, of the many changes underway in the world today, one of the most striking is the decline and homogenization of marine biodiversity through overexploitation, climate change, habitat loss, nutrient pollution and species invasions (Meyer et al., 2016). This raises considerable concern, since human society depends on the ocean for its services and the loss of marine biodiversity has detrimental impacts on the environment. For example, Villnas et al. (2013) demonstrated that repetitive anthropogenic stressors result in the gradual degradation of individual functions in sedimentary ecosystems. Additionally, based on (a) experimental systems and (b) individual functions, Gamfeldt et al. (2015) concluded that loss of marine biodiversity likely decreases ecosystem function and impacts habitat resilience. Specific to coral reefs, Lefcheck et al. (2019) demonstrated that tropical fish diversity significantly impacts functioning and Brandl et al. (2019) observed the same relationship between biodiversity and ecosystem functioning. In the last decade, an increasing number of studies have arrived at similar conclusions about the biodiversity-ecosystem function relationship in natural systems (non-experimental, 'real-world' ecosystems;drylands -Maestre et al., 2012;Berdugo et al., 2017, grasslands -Soliveres et al., 2016, forests -Zhang and Wang, 2012. Our understanding of ecosystem function, however, is not only far from complete but is also biased toward experiments that have been conducted (1) over a small-scale, (2) with a few species (3) across a few lower trophic levels Nelson et al., 2009;Berdugo et al., 2017) and (4) are biased toward single functions (Stachowicz et al., 2007;Cardinale et al., 2011;Byrnes et al., 2014). In the last decade, efforts have been made to change this trend, by expanding ecosystem function studies across spatial scales (Thompson et al., 2018;Gonzalez et al., 2020). A recent meta-analysis by Duffy et al. (2017) synthesizes these biodiversity effects from multiple ecosystem types. Specifically, Ptacnik et al. (2008); Mora et al. (2011), Zimmerman and Cardinale (2013); Duffy et al. (2015) Duffy et al. (2016), and García-Comas et al. (2016) present measured effects of eelgrass, fish species and plankton diversity on reef ecosystem functioning after statistically controlling for environmental covariables.
Ecosystem function components, when considered synergistically, or additively can have different and likely stronger effects on biodiversity (Hector and Bagchi, 2007;Hautier et al., 2017;Meyer et al., 2018). For instance, Maestre et al. (2012) used a natural experiment to test the biodiversity-ecosystem function relationship (BEF; Duffy, 2009;Gamfeldt and Roger, 2017) between plant species richness and multiple functions in semi-arid ecosystems. This proved to be timely because previous work on dryland systems had only assessed ecosystem functions individually. This led to the concept of a habitat's "multifunctionality" which was defined as the "simultaneous provision of multiple functions" (Hector and Bagchi, 2007;Gamfeldt et al., 2008). Recently, in a study by Byrnes et al. (2014), the authors suggested that aspects of the BEF relationship, which continue to remain unresolved in several marine ecosystems, could be tested stringently when multifunctionality is taken into account and quantified. Taking this forward, Lefcheck et al. (2015) highlighted the importance of taking multiple ecosystem functions into consideration to understand the effect of biodiversity on integrated functioning, an aspect that is unclear when only individual functions are analyzed. By presenting evidence of biodiversity's effect on multifunctionality, the authors demonstrated that communities with higher species richness maintain a higher number of ecosystem functions than those with lower species richness. In this manner, the concept of multifunctionality became useful to understand ecosystem functioning in fundamental and applied ecology.
The measurement of multifunctionality, however, has proven extremely challenging as it takes into consideration only a small subset of all possible individual ecosystem functions. To our knowledge, there is no standardized way of defining or carefully accounting for specific functions that will capture an ecosystem's "true" multifunctionality (Manning et al., 2018). For instance, Alsterberg et al. (2014) illustrated that there were no effects of nutrient enrichment, toxicants, sedimentation and warming on marine ecosystem multifunctionality. However, a few years later, Alsterberg et al. (2017) demonstrated that habitat diversity has a direct effect on marine ecosystem functionality. As a result, data on ecosystem functioning of marine ecosystems are especially limited and fragmented, despite the fact that such information could help inform effective management of crucial environments (Manning et al., 2018).
The last decade has seen a surge in the development of methods to quantify ecosystem multifunctionality (Maestre et al., 2012 -the averaging approach;Gamfeldt et al., 2008;Byrnes et al., 2014 -single and multiple threshold approaches). Manning et al. (2018) introduced two additional components -ecosystem function (EF) multifunctionality and ecosystem services (ES) multifunctionality to help environmental managers and economists quantify the value of ecosystems using monetary units and life satisfaction. Recently, Hoölting et al. (2019) proposed a new way of measuring multifunctionality across spatial scales. Although these methods have been reasonably successful in assessing an ecosystem's performance (Hensel and Silliman, 2013 reviews a few of the methods listed), the literature cited have not yet completely addressed how biodiversity when coupled with environmental factors influences ecosystem function, especially for urbanized coastal marine ecosystems. Therefore, it is crucial to understand the resilience of coastal marine ecosystems in response to global change and environmental stressors. This is especially relevant to a coastal megacity such as Hong Kong, where eutrophication (nutrient-driven marine pollution) has contributed to a marked reduction in water quality  and loss of critical habitats such as hard corals, thus reducing the complexity, diversity and function of benthic ecosystems (Scott, 1990;Fabricius and McCorry, 2006). For example, Tolo Harbour, in northeast Hong Kong was once a pristine, coral-fringed bay home to various coral communities and vibrant fisheries (Morton, 1989). However, several sites in the harbor have experienced a dramatic loss of coral cover from >60% to <10%. In 2015, we utilized stable isotope analysis to understand nitrogen source dynamics in both wastewater effluents and receiving seawaters in the Tolo Channel and found sewage effluent to be the dominant source of nitrogen pollution (Archana et al., 2016(Archana et al., , 2018. Consistently, Wong et al. (2017) measured δ 15 N of hard corals from several sites in this region, which also revealed strong human-derived nutrient signals. Today, we see a punctuated gradient in water quality with nutrient concentrations decreasing with distance from coastal populations. Yet we have little understanding of how these eutrophic conditions affect the ecosystem function of such benthic marine communities and their multifunctionality.
Here, we analyzed ecosystem processes -primary productivity, herbivory, predation, organic matter decomposition, carbon sequestration and overall ecosystem multifunctionality along a gradient of environmental parameters in an urbanized coastal marine environment. We used environmental data (nutrient concentrations, dissolved oxygen, chlorophyll A, pH, salinity, turbidity, secchi depth, temperature, and suspended solids) to test the hypothesis that, ecosystem multifunctionality is significantly correlated with geochemical parameters. Next, we used abundance data of functional groups from annotated settlement plate photos of autonomous reef monitoring structures (ARMS) deployed in the same sites to test the hypothesis that increased abundance of benthic taxa is significantly correlated with ecosystem functioning. We used fish abundance data from eDNA of water samples collected from the same sites to test the hypothesis that increased fish abundance significantly alters ecosystem functioning.

Study Area
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
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: Where i and f refer to initial and final wet weights as per Rasher et al. (2013).

Herbivory/Grazing Intensity
Herbivores play a vital role in maintaining a healthy coraldominated 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: 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.

Predation
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" . 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 FIGURE 2 | Ecosystem function assays for (a) primary productivity using Ulva in closed plastic bottles (b) herbivory using Ulva tied to ropes (c) green and rooibos tea bags shown before decomposition and upon retrieval (d) after 6 months of deployment (e) and (f) predation using squid pops.
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 ; Figure 2).
Where i and f refer to the initial and final number of squid pops, respectively.

Organic Matter Decomposition and Carbon Sequestration
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, costeffective 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:

Multifunctionality
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.

ARMS Plate Photo Annotations
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.

Seawater Sampling
Seawater samples were collected in August 2016 from the surface in 1-L Nalgene bottles (3 replicates per site from the four study sites) and filtered through a 47 mm diameter PES filter (nominal pore size, 0.22 µm; Millipore Express). Filters were wrapped in aluminum foil and stored at −20 • C for DNA extraction.

DNA Extraction
DNA was extracted in triplicate from every sampling site, using the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer protocol and supplemented with a few steps: first, the filter was cut into 4 pieces and put into a UV-sterilized 1.5 mL tube. Next, 180 µL ATL was added to each tube. Following this, 20 µL Proteinase K was added to each tube individually and vortexed. The tubes were vortexed a few more times during the day and incubated at 55 • C overnight. To increase the yield, the DNA was eluted with 100 µL Buffer AE twice. Next, following incubation at room temperature for 5 min, the spin column was spun at 8000 RPM for 1 min. For the remaining steps, we followed the manufacturer's protocol. The extracted DNA was stored at −20 • C. To monitor contamination, a blank filter was used as negative control. DNA from four fish species (Nemipterus bathybius, Gambusia affinis, Xiphophorus helleri, and Macropodus opercularis) was extracted and used as the positive control.

Library Preparation
We chose a set of universal PCR primers for metabarcoding eDNA from fishes in subtropical habitats (MiFish-U). The primers amplify an approximately 163-185 bp region of the 12S rRNA gene, which is known to resolve taxonomy of most fishes to the family level (Miya et al., 2015). We used a two-step PCR amplification protocol. In the first round of PCR amplification, 12 samples with the universal fish primer pairs (MiFish-U) were used. The PCR amplification was done in 30 cycles and the total reaction volume was 12 µL, containing 6 µL 2 × KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, MA, United States), 0.36 µL of each primer (10 µM), 2.0 µL of template and 3.28 µL of sterile distilled H 2 O. The thermal cycle profile was as follows: initial 3 min denaturation step at 95 • C followed by denaturation at 98 • C for 30 s; annealing at 63 • C for 30 s; extension at 72 • C for 20 s; and a final extension at 72 • C for 5 min. Collected PCR products were purified using the MinElute Kit (Qiagen) with an elution volume of 10 µL. DNA fragments were recovered and ready to use in downstream analysis. The concentration of DNA was measured using a Thermo Scientific TM µDropTM Plate. The second round of PCR followed by library preparation and sequencing were performed at Genewiz (Suzhou, China) following protocols described in Miya et al. (2015). Indexed amplicon libraries were constructed, pooled and multiplexed sequenced using the Illumina MiSeq platform.

Sequencing Data Optimization
For base calling and preliminary analysis of the raw data, we used Bcl2fastq (v2.17.1.14). To optimize the raw sequencing data -(a) two sequences of each read pair were merged according to overlapping sequences, after which undetermined bases were removed (b) Primers and adapter sequences were removed (c) the 5 and 3 bases with Qphred score < 20 were trimmed (d) resulting sequences >200 bp in length passed this stage of optimization. Cutadapt (v1.9.6), Qiime (1.9.1), and Vsearch (1.9.6) were used (see Supplementary Information for raw data statistics).

Taxonomic Assignment
OTU refers to an operational definition of a classification unit (genus, species, grouping, etc.) used commonly for data analysis. All the sequences in a sample were classified to obtain information on species and genus. By classification, the sequences were grouped according to their similarity, and one group is an OTU. First, unique sequences were extracted from optimized sequences with read count information. Next, OTU clustering of unique sequences (read count > 1) was performed with similarity of 97%, and chimeric sequences were further removed to obtain representative OTU sequences. The RDP classifier was used to select and annotate the representative sequence for each OTU to obtain the community composition of each sample. Community composition was subsequently analyzed and summarized using Silva_128 12S rRNA database for species annotation. The relationship between fish community composition (Supplementary Appendix Tables A.4, A.5) and ecosystem function was characterized to understand the factors driving its variability.

Environmental Data
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

Data Analysis
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.

RESULTS
We present our analyses on the level of individual functions and joint ecosystem functioning (or multifunctionality) in four habitats that have varying environmental parameters, fish abundance and benthic community composition. We analyzed the effect of urbanization, as reflected in the environmental variables on multiple ecosystem functions -primary productivity, predation, herbivory, decomposition and carbon storage. Next, we use two methods to assess multifunctionality -the average approach and the 0.5 threshold approach. Finally, we evaluate the effects of fish abundance and benthic community composition on individual ecosystem functions and overall multifunctionality.

Individual Ecosystem Functions
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 FIGURE 4 | Ecosystem function data for primary productivity, herbivory, predation, organic matter decomposition and carbon sequestration across the four study sites (Center Island, Che Lei Pai, Port Island, Tung Ping Chau). Primary productivity was measured as the % gain in mass of macroalgae over a deployment period of 48 h; herbivory was measured as the % loss in mass of macroalgae over a deployment period of 48 h; predation was measured as the % loss in bait (dried squid) deployed over a period of 1 h -note: data missing from TPC; organic matter decomposition and carbon sequestration were calculated using the tea bag assay (refer to section "Materials and Methods"). Data is represented as the average value ± standard deviation. Connecting letters indicate significant differences between sites revealed by post hoc Tukey HSD tests (α = 0.05).
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 2 = 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 2 = 0.96). Temperature and coral cover (p < 0.05) were important in predicting predation intensity (best-fit model r 2 = 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 2 = 0.44).

Analysis of Plate Imagery: Abundance of Sessile Fauna on ARMS Plates
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 algaechlorophytes, 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 2 = 0.95).

eDNA Data on Fish Abundance
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.

Overall Ecosystem Multifunctionality
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).

DISCUSSION
In the last decade, research into the concept of multifunctionality (ability of ecosystems to provide multiple ecosystem functions) has increased owing to tremendous anthropogenic pressure on natural resources. This has resulted in the need to design and manage ecosystems so that they can efficiently provide multiple ecosystem functions simultaneously (Manning et al., 2018). However, not all ecosystem functions exhibit the same response to environmental drivers (Bradford et al., 2014). Given the complexity and diversity of ecosystem functions, the present study discusses five vital functions (primary productivity, herbivory, predation, organic matter decomposition and carbon sequestration) in an urbanized coastal marine environment. Data from each individual ecosystem function has been used to derive the overall multifunctionality for the human-impacted habitat. Linkages between multifunctionality, individual functions, environmental factors, fish abundance and benthic community composition datasets have been presented here.

Habitat as a Driver of Ecosystem Function
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 as a Driver of Ecosystem Function
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 FIGURE 6 | Results of relative abundance of predatory, herbivorous, planktivorous and omnivorous fish obtained from eDNA analysis of seawater samples collected from the same study sites. We recorded no fish in Center Island. Data is represented as the average value ± standard deviation. 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 filterfeeders 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 .

Abiotic Factors as Drivers of Ecosystem Function
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 Coral-Algae-Herbivore Triangle
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 reasonsfirstly, 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 coralalgal-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.

Relative Feeding Intensity of Generalist Predators
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.

Intermediate Disturbance Affects Overall Ecosystem Multifunctionality
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.

CONCLUSION AND IMPLICATIONS FOR CONSERVATION MANAGEMENT
The results from this study have important implications for understanding processes regulating the structure of human-impacted coastal marine habitats. By characterizing important ecosystem functions (primary productivity, herbivory, predation, decomposition and carbon sequestration) and overall multifunctionality, the findings significantly expand on what we already know about urbanization impacts on key ecosystem processes. The highlight of this study is the finding that urbanization impacts multifunctionality in coastal marine systems. However, these results warrant a more detailed investigation on the links between multifunctionality and community composition (microbial and metazoan). Further research into the role of decomposition in overall ecosystem multifunctionality and associated carbon storage mechanisms can be used to help in global carbon inventory management. Moreover, research into functional level responses, microbial interactions and additional ecosystem services (Manning et al., 2018) such as shoreline stabilization, biogeochemical cycling and carbon sequestration will help shed light on the holistic properties of the ecosystem.

DATA AVAILABILITY STATEMENT
All data generated for this study has been made publicly available in the Appendix and Supplementary Material of this article.

AUTHOR CONTRIBUTIONS
AA and DB conceptualized the research. AA conducted the fieldwork and did the analysis. DB reviewed the manuscript and analysis. Both authors contributed to the article and approved the submitted version.

FUNDING
We are grateful to the International Coral Reef Society for the graduate student fellowship awarded to AA that helped fund this work. We are grateful to the Environment Conservation Fund (ECF, Grant 67-2016) and Collaborative Research Fund (CRF, C7013-19GF) that helped fund the work on the autonomous reef monitoring structures (ARMS) and eDNA. This paper is contribution 67 to the Smithsonian's MarineGEO-TMON program.