Site-Level Variation in Parrotﬁsh Grazing and Bioerosion as a Function of Species-Speciﬁc Feeding Metrics

: Parrotﬁsh provide important ecological functions on coral reefs, including the provision of new settlement space through grazing and the generation of sediment through bioerosion of reef substrate. Estimating these functions at an ecosystem level depends on accurately quantifying the functional impact of individuals, yet parrotﬁsh feeding metrics are only available for a limited range of sites, species and size classes. We quantiﬁed bite rates, proportion of bites leaving scars and scar sizes in situ for the dominant excavator ( Cetoscarus ocellatus, Chlorurus strongylocephalus, Ch. sordidus ) and scraper species ( Scarus rubroviolaceus, S. frenatus, S. niger, S. tricolor, S. scaber, S. psittacus ) in the central Indian Ocean. This includes the ﬁrst record of scar frequencies and sizes for the latter three species. Bite rates varied with species and life phase and decreased with body size. The proportion of bites leaving scars and scar sizes di ﬀ ered among species and increased with body size. Species-level allometric relationships between body size and each of these feeding metrics were used to parameterize annual individual grazing and bioerosion rates which increase non-linearly with body size. Large individuals of C. ocellatus, Ch. strongylocephalus and S. rubroviolaceus can graze 200–400 m 2 and erode > 500 kg of reef substrate annually. Smaller species graze 1–100 m 2 yr − 1 and erode 0.2–30 kg yr − 1 . We used these individual functional rates to quantify community grazing and bioerosion levels at 15 sites across the Maldives and the Chagos Archipelago. Although parrotﬁsh density was 2.6 times higher on Maldivian reefs, average grazing (3.9 ± 1.4 m 2 m − 2 reef yr − 1 ) and bioerosion levels (3.1 ± 1.2 kg m − 2 reef yr − 1 ) were about 15% lower than in the Chagos Archipelago (4.5 ± 2.3 and 3.7 ± 3.0, respectively), due to the dominance of small species and individuals in the Maldives (90% < 30 cm length). This demonstrates that large-bodied species and individuals contribute disproportionally to both grazing and bioerosion. Across all sites, grazing increased by 66 ± 5 m 2 ha − 1 and bioerosion by 109 ± 9 kg ha − 1 for every kg increase in parrotﬁsh biomass. However, for a given level of parrotﬁsh biomass, grazing and bioerosion levels were higher on Maldivian reefs than in the Chagos Archipelago. This suggests that small-bodied ﬁsh assemblages can maintain ecosystem functions, but only if key species are present in su ﬃ ciently high numbers. the Maldives (Gaafu Dhaalu), which are dominated by small parrotﬁsh ( < 10% of ﬁsh > 30 cm), and the Chagos Archipelago (Peros Banhos and Salomon), which are dominated by large parrotﬁsh (30-90% of ﬁsh > 30 cm). Equations and R 2 describe bootstrapped regressions models. Abbreviations: G–Grazing, E–Bioerosion, B–Biomass.


Introduction
Parrotfishes (family Labridae, subfamily Scarinae) are an abundant, diverse and ecologically important group of coral reef fishes. Although commonly characterized as herbivores, they feed on cyanobacteria and other autotrophic microorganisms on (epilithic) or within (endolithic) calcareous substratum by cropping, scraping or biting the reef surface with their beak-like jaws [1,2]. Through this feeding strategy, parrotfish provide critical ecosystem functions on coral reefs, including the provision of settlement space for corals and other benthic organisms by clearing benthic substrate, and acting as primary agents for bioerosion and sediment production (reviewed by [3]). However, parrotfish are not a uniform group in their feeding behavior and functional impacts. Species with an "excavating" feeding mode typically erode more substrate than "scrapers" due to larger bite volumes, and larger individuals have a disproportionately greater grazing and bioerosion effect than smaller individuals (e.g., [4][5][6]).
With the rapid degradation of coral reef ecosystems, the quantification and preservation of reef ecosystem functions or processes, defined as the movement of material or energy, has gained increasing attention [7,8]. Based on parrotfish abundance and functional traits, studies have estimated and explored drivers of site-specific functions. For instance, Hoey and Bellwood [9] quantified cross-shelf variation in parrotfish densities and ecosystem functions on the Great Barrier Reef (GBR) and found that the differences in scraping, bioerosion and coral predation were driven by only a few species. Robinson et al. [10] looked at potential drivers of herbivory across the Indo-Pacific and found that cropping rates are influenced by substrate availability and structural complexity, while scraping rates are mainly controlled by fishing pressure. In contrast, Bellwood et al. [11] found that the processes of bioerosion and coral predation are highly sensitive to human activity, while scraping and sediment removal are not. Independent of direct human pressures, the presence of sea birds and associated nutrient inputs enhanced parrotfish biomass on adjacent reefs and thereby increased rates of bioerosion and grazing >3-fold compared to islands without birds [12]. Similarly, increases in parrotfish abundance and a changing population size structure after the 2015-2016 mass bleaching event led to increased bioerosion-derived sediment generation in the Maldives over several years post-bleaching [13]. In general, the quantification of ecosystem functions has important implications for reef management, and species-specific feeding metrics, normalized to grazable area of reef, have been shown to be better predictors of benthic change than fish biomass alone [14].
In order to infer ecosystem processes from fish assemblages, the empirical measurement of underlying individual functional rates is essential [7]. This means that to estimate community-level grazing and bioerosion, the impact of each individual parrotfish has to be quantified as accurately as possible to avoid inflation of errors during upscaling. Several metrics are required in this respect: (i) bite rate; (ii) the proportion of bites leaving scars; and (iii) scar area/volume. While bite rate data is available for approximately 55% of parrotfish species, data on scar sizes and the proportion of bites leaving scars are only available for 32% and 20% of species, respectively (available metrics reviewed in [15]). These ratios are much lower in the Indo-Pacific compared to the Caribbean, as species numbers are higher (69 compared to 13 species). Reasons for the scarcity of detailed feeding metrics for most species and locations include: the high expenditure of time to collect in situ observations; the large size range of many species that needs capturing; the variability in community composition across locations and habitats; low abundances or schooling behavior of some species; and the fact that for these reasons many studies have focused only on two to three dominant species. Difficulties in accurately measuring bite scar sizes additionally lead to a bias in data collection towards large species and individuals. These issues aside, previous studies in the Indo-Pacific region have shown that feeding metrics vary with species, life phase and body size [10,[16][17][18][19][20], and that functional impacts may differ regionally [3,21] and seasonally [22,23]. Furthermore, high wave exposure [24,25] and high sediment load [26,27] can negatively affect feeding rates. Given the scarcity of species-specific estimates of feeding metrics, the assessment of grazing and bioerosion impacts of parrotfish assemblages at a given site most often relies on genus-or group-level estimates or data obtained in different regions. This suggests Diversity 2020, 12, 379 3 of 16 that the quantification of ecosystem functions is inaccurate for many globally important reef regions. Additional species-and size-specific data informing functional rates are urgently required to better understand how and why the key ecosystem functions that parrotfish provide vary among locations.
In this study we: (i) Quantify the relative effects of "species", "life phase" and "body size" as drivers of individual parrotfish feeding metrics (i.e., bite rate, proportion of bites leaving scars, scar area/volume) and parameterize species-level relationships between body size and each of these feeding metrics; (ii) use these allometric feeding metric equations to compute species-and functional group-specific annual grazing and bioerosion rates; and (iii) combine these individual functional rates with visual census data to calculate and compare community grazing and bioerosion levels across three atolls in the Maldives and Chagos Archipelago.

Study Sites and Parrotfish Species
Feeding metrics were collected from a range of shallow water reef habitats (2-8 m depth) at five atoll interior reef sites in the southern Maldives (Kandahalagalaa, Kafigalaa, Maahutigalaa, Kodehutigalaa and Kadumaigalaa in Gaafu Dhaalu atoll) and at twelve fore reef and lagoonal sites in the Chagos Archipelago (across the atolls of Salomon, Peros Banhos, Great Chagos Bank, Egmont and Diego Garcia) ( Figure 1). All sites were affected by extensive coral mortality during the 2015-2016 bleaching event, decreasing coral cover from >25% in 2015/2016 to <10% in 2018/2019 [13,28]. Local parrotfish populations benefitted from the resulting change in food availability and increased in abundance and body size at both locations [13,29]. As a result, most of the substrate at both locations was occupied by closely cropped filamentous turf algae and calcareous algae at the time of study. Site-specific differences in sediment load and wave exposure were not quantified, but were generally low at all sites of feeding observations and therefore unlikely to influence parrotfish feeding behavior.
Diversity 2020, 12, x FOR PEER REVIEW 3 of 16 globally important reef regions. Additional species-and size-specific data informing functional rates are urgently required to better understand how and why the key ecosystem functions that parrotfish provide vary among locations. In this study we: i) Quantify the relative effects of "species", "life phase" and "body size" as drivers of individual parrotfish feeding metrics (i.e., bite rate, proportion of bites leaving scars, scar area/volume) and parameterize species-level relationships between body size and each of these feeding metrics; ii) use these allometric feeding metric equations to compute species-and functional group-specific annual grazing and bioerosion rates; and iii) combine these individual functional rates with visual census data to calculate and compare community grazing and bioerosion levels across three atolls in the Maldives and Chagos Archipelago.

Study Sites and Parrotfish Species
Feeding metrics were collected from a range of shallow water reef habitats (2-8 m depth) at five atoll interior reef sites in the southern Maldives (Kandahalagalaa, Kafigalaa, Maahutigalaa, Kodehutigalaa and Kadumaigalaa in Gaafu Dhaalu atoll) and at twelve fore reef and lagoonal sites in the Chagos Archipelago (across the atolls of Salomon, Peros Banhos, Great Chagos Bank, Egmont and Diego Garcia) ( Figure 1). All sites were affected by extensive coral mortality during the 2015-2016 bleaching event, decreasing coral cover from >25% in 2015/2016 to <10% in 2018/2019 [13,28]. Local parrotfish populations benefitted from the resulting change in food availability and increased in abundance and body size at both locations [13,29]. As a result, most of the substrate at both locations was occupied by closely cropped filamentous turf algae and calcareous algae at the time of study. Site-specific differences in sediment load and wave exposure were not quantified, but were generally low at all sites of feeding observations and therefore unlikely to influence parrotfish feeding behavior. Fish observations were performed during 09:00-16:00 h, the peak feeding period for parrotfishes (e.g., [5,30]). Feeding metrics were quantified with SCUBA or snorkeling using focal individual Fish observations were performed during 09:00-16:00 h, the peak feeding period for parrotfishes (e.g., [5,30]). Feeding metrics were quantified with SCUBA or snorkeling using focal individual observations: An adult parrotfish ( > 10 cm total length (TL)) was haphazardly selected and followed for a short period of acclimation (ca. 1 min), during which species, life phase (initial or terminal) and size (to the closest 2.5 cm size class: 10, 12.5, 15 cm, etc.) was noted. Studied species were the most common parrotfishes on the visited central Indian Ocean reefs, with observed size ranges given in brackets: Cetoscarus ocellatus  15-20 cm) are included in the raw data tables but not presented in this paper, as the replication over life phases and sizes was not sufficient to enable detailed analyses. Species were classified into functional groups of excavators (Cetoscarus sp. and Chlorurus spp.) and scrapers (Scarus spp.) [4,31]. Individuals <10 cm were not included in feeding observations as they are considered to make negligible contributions to grazing and bioerosion levels [5,16].

Feeding metrics
To quantify bite rates, each fish was observed for 3-5 min, the number of bites was recorded using a handheld tally counter, and counts were converted to bites min −1 (bpm). Observations were discontinued if the fish showed a detectable response to the diver. Replicates for each species (average n = 44) were spread as evenly as possible over both life phases and all observed body sizes.
To quantify the proportion of bites leaving scars (%scars), a fish was followed until it took a number of successive bites (foray) in an area that could be clearly observed. Immediately, this area was inspected closely to determine how many of the recorded bites left visible scars. Data were modelled as %scars in individual feeding forays for each species (average n = 37) except for S. frenatus, for which replication was low.
If bite scars were visible in the observed bite forays, the maximum dimensions of scars were measured using calipers. Scar areas (area, in mm 2 ) were calculated by multiplying length and width of each bite (average n = 34). Data for C. ocellatus and S. frenatus were not analyzed as replication was low.
It was difficult to accurately measure the depth of each bite scar due to variable substrate morphology. For the purpose of calculating scar volume and bioerosion, depth of bite scars was therefore assumed to be 1.5 mm for C. ocellatus and Ch. strongylocephalus, 1 mm for S. rubroviolaceus ≥ 40 cm TL, 0.2 mm for Ch. sordidus and 0.1 mm for all other species and S. rubroviolaceus < 40 cm (own observations and [4,30]). Scar volumes (vol, in cm 3 ) were calculated assuming a cuboid form (vol = area × depth), which may have resulted in a slight overestimation, although no significant differences were found when comparing this method to direct volume measurements via wax casts of bite scars in past studies [17,30,32].
Data on feeding metrics were pooled over life phases and modelled as a function of body size (TL) to obtain species-and functional group-specific best-fit regression lines (power functions y = a × (TL) b for bpm, area and vol; logarithmic function y = a + (b × ln(TL)) for %scars). Although bite rates were significantly lower for terminal compared to initial phase individuals for six out of nine species, data on bpm were also pooled over life phases because i) life phase only explained an additional 0.75% of the variability in the fitted model (Table S1), and ii) the effect of life phase could not be disentangled from body size effects, as the bigger size of terminal individuals also decreases bite rates ( Figure S1).

Individual Grazing and Bioerosion Rates
The equations describing the allometric relationships between body size and feeding metrics were used to derive average bpm, %scar, area and vol for individuals within 10 cm size classes (10-19, 20-29, 30-39, 40-49, 50-59, >60 cm) using midpoints of size classes as TL (Table S2). Species-specific annual grazing and bioerosion rates per individual were then calculated with the following formulas: Grazing rate (m 2 ind. −1 yr −1 ) = area (mm 2 )/10 6 × %scars × bpm (min −1 ) × 60 (min) × hours of daylight × proportion of daytime feeding × 365 days (1) As the study sites are located very close to the equator, 12 h of daylight was used for calculations. Proportion of daytime feeding was adapted from Bellwood [30]: 83.3% for large parrotfish (C. ocellatus, Ch. strongylocephalus, S. rubroviolaceus) and 87.7% for smaller species (Ch. sordidus and other Scarus spp.). For (2), an average substrate density of 1.52 ± 0.19 g cm −3 was used, which was determined as the average bulk density of 138 coral skeleton samples from 23 different genera collected in the Chagos Archipelago in 2019 (ID Lange, unpublished data). Average coral density values measured in the Maldives are very similar (1.50 ± 0.19 g cm −3 [33]).

Community Grazing and Bioerosion Levels
Annual parrotfish grazing and bioerosion levels at a specific site (m 2 reef grazed or kg carbonate eroded per m 2 reef per year) were calculated by combining individual grazing and bioerosion rates with data on species' densities and sizes. Fish abundances were determined via underwater visual census (UVC) surveys at the same atoll interior reef sites where feeding observations were conducted in Gaafu Dhaalu atoll, Maldives (n = 5) in January 2019 (30 × 5 m belt transects, n = 8 per site) and at ten fore-reef sites in the Chagos Archipelago across the atolls of Peros Banhos (n = 6) and Salomon (n = 4) in March 2019 (50 x 5 m belt transects, n = 4 per site) ( Figure 1). Each observed parrotfish ≥10 cm was identified to species level and TL was estimated to the closest 10 cm size class (Maldives) or the closest cm (Chagos Archipelago). Overall species composition and abundance of parrotfish communities in 2019 were consistent with previous surveys at the same locations [13,34].
Biomass was calculated for each observed individual using TL or the mid-point of the corresponding size class with species-or genus-specific length-weight relationships obtained from FishBase (Biomass (kg) = a × (TL) b / 1000 [35]). Annual grazing and bioerosion rates were calculated for each individual parrotfish using formulas (1) and (2) and the equation parameters describing the allometric relationships between TL and each of the feeding metrics (Table S2). Four species that were occasionally observed in the censuses (together accounting for <7% of total individuals), for which feeding metrics were not quantified in this study, were assigned rates determined for S. frenatus (S. caudofasciatus, S. prasiognathos) or a combination of species-specific bpm extracted from raw data of Robinson et al. [10] and %scars, area and vol functions obtained in this study (Hipposcarus harid, S. viridifucatus). The biomass or grazing/bioerosion rates of individuals within each transect were summed up and transect-level values were averaged to yield site-level data.

Statistical Analysis
All analysis were performed in R 3.5.1 [36]. The relative impacts of predictor variables "species" (factor with 9 levels, see above), "life phase" (2 levels: initial, terminal) and "body size" (continuous) on each feeding metric (bpm, %scar, area, vol) were assessed using generalized linear models (GLMs), thereby allowing for different error distributions of the response variables. Based on data type and the evaluation of model residual plots, quasipoisson error distribution with "log" link was used for bpm (count data, ratio, overdispersion), quasibinomial distribution with "logit" link for %scars (proportional data), and Gamma error distribution with "log" link for area and vol (continuous data). The adequacy of these models was evaluated by examining (i) the plot of residuals vs. fitted values to look for heteroscedasticity and (ii) the normality Q-Q plot to assess the normality of errors [37]. Interaction terms were included in the final model if significant (Table S1). The variable "functional group" (2 levels: excavator, scraper) was collinear with "species" and was therefore, excluded from GLMs and instead tested separately for impacts on each feeding metric using one-way ANOVAs. Collinearity between other predictors was inspected using variance inflation factors (VIF) and found to be negligible (VIF < 2). Factor "species" was subject to Tukey's HSD pairwise comparisons to determine which species differed significantly from others.
To estimate the uncertainty of individual grazing and bioerosion rates, we used the Monte-Carlo method of error propagation [38]. For each species and feeding metric, 1000 possible regression lines were computed within the error probability distribution of the fitted model using mvrnorm() [39], and their coefficients were extracted. Using random coefficients from each metric, the calculation of grazing and bioerosion rates for fish sizes 10 to 60 cm was repeated 1000 times and a power function was fitted to data from each iteration ( Figure S2). The resulting distribution of model coefficients was used to extract the mean and standard deviation for the final grazing/bioerosion regression and confidence interval.
Community-level density, biomass and annual grazing/bioerosion of parrotfish assemblages at sites in the Maldives and Chagos Archipelago were compared at atoll level using non-parametric Kruskal-Wallis and, if significant, post-hoc Dunn tests, as a non-parametric comparison of medians was more appropriate for the sample size and error distribution of data.
To test if parrotfish biomass is a good proxy for community-level ecosystem functions, grazing/bioerosion levels were modelled as a function of transect-level biomass using bootstrapping (R = 2000) as a non-parametric alternative to linear regression analysis [40]. Models were fitted for (i) all sites, (ii) Maldives-communities dominated by small parrotfish (<10% of fish >30 cm), and (iii) Chagos Archipelago-dominated by large parrotfish (30-90% of fish >30 cm) (following [10]). One outlier was removed for this analysis, as a school of 28 S. prasiognathos in a single transect at a Maldivian reef site caused unrepresentatively high biomass values.

Feeding metrics
Bite rates (bpm) differed among species (Figure 2A, Table S1) and were significantly higher in initial compared to terminal life phases for all species except Chlorurus sordidus, Scarus rubroviolaceus and S. scaber ( Figure S1B). Furthermore, bpm significantly decreased with body size for all species except Cetoscarus ocellatus, S. frenatus and S. scaber ( Figure 2A, Table S2). Comparing functional groups, bpm were generally higher for scrapers compared to excavators (F 1,416 = 130.5, p < 0.001, Figure 2B).
The proportion of bites leaving scars (%scars) differed significantly between functional groups (F 1,314 = 38.96, p < 0.001, Figure 2D) and among species ( Figure 2C, Table S1), the latter mainly driven by differences between S. scaber and S. psittacus compared to large-bodied parrotfish (Tukey HSD p < 0.05). Body size had a significant positive effect on %scars, with Chlorurus spp., S. rubroviolaceus and S. tricolor showing more pronounced increases in %scars with increasing fish size compared to other Scarus species (Table S2).
Scar areas (area) were mostly influenced by body size, as seen from overlapping confidence intervals comparing species and functional groups ( Figure 2E, F). However, functional group (F 1,365 = 44.37, p < 0.001) and species (Table S1) also had a significant influence on area, with Ch. strongylocephalus and S. rubroviolaceus leaving larger bite scars than all other species and Ch. sordidus leaving larger scars than S. niger and S. scaber (Tukey HSD p < 0.05).
Scar volumes (vol) were also influenced by body size and species (Table S1), with a clearer distinction of Ch. strongylocephalus and S. rubroviolaceus allometric relationships due to differences in bite depth ( Figure 2G). As a result, differences between functional groups are also more pronounced than for area (F 1,365 = 54.71, p < 0.001; Figure 2H).

Individual Grazing and Bioerosion Rates
Large individuals (>50 cm TL) of C. ocellatus, Ch. strongylocephalus, S. rubroviolaceus and S. frenatus are able to graze 200-400 m 2 of reef annually. However, small individuals of these species exert lower grazing pressure than similar-sized scrapers and Ch. sordidus, as bite rates of the latter are usually higher (Table 1, Figure 3A). Annual bioerosion rates can exceed 500 kg yr −1 for large individuals of C. ocellatus, Ch. strongylocephalus and S. rubroviolaceus. Rates for smaller Scarus species are similar to each other within the same size classes, except for S. psittacus, which showed very small bioerosion rates

Individual Grazing and Bioerosion Rates
Large individuals (>50 cm TL) of C. ocellatus, Ch. strongylocephalus, S. rubroviolaceus and S. frenatus are able to graze 200-400 m 2 of reef annually. However, small individuals of these species exert lower grazing pressure than similar-sized scrapers and Ch. sordidus, as bite rates of the latter are usually higher (Table 1, Figure 3A). Annual bioerosion rates can exceed 500 kg yr −1 for large individuals of C. ocellatus, Ch. strongylocephalus and S. rubroviolaceus. Rates for smaller Scarus species are similar to each other within the same size classes, except for S. psittacus, which showed very small bioerosion rates due to low percentages of bites leaving scars. Ch. sordidus erodes more substrate than similar sized Scarus species due to larger and deeper bite scars (Table 1, Figure 3C). As large individuals of S. rubroviolaceus leave large bites, and because the grazing and bioerosion rates of Ch. sordidus are more similar to scrapers than larger excavators, the categorization into functional groups hides much of the variability in grazing and bioerosion rates among species (Figure 3B,D). Table 1. Individual annual grazing (m 2 ind. −1 yr −1 ) and bioerosion rates (m 2 ind. −1 yr −1 ) for a range of parrotfish species and size classes. Empty cells indicate size classes above the maximum total length (TL) of species (based on own observations and FishBase). Coefficients a and b describe allometric regressions between body size and individual functional rates (grazing/bioerosion = a × (TL) b ).
Grazing (m 2 ind. − due to low percentages of bites leaving scars. Ch. sordidus erodes more substrate than similar sized Scarus species due to larger and deeper bite scars (Table 1, Figure 3C). As large individuals of S. rubroviolaceus leave large bites, and because the grazing and bioerosion rates of Ch. sordidus are more similar to scrapers than larger excavators, the categorization into functional groups hides much of the variability in grazing and bioerosion rates among species ( Figure  3B, D). Table 1. Individual annual grazing (m 2 ind. -1 yr -1 ) and bioerosion rates (m 2 ind. -1 yr -1 ) for a range of parrotfish species and size classes. Empty cells indicate size classes above the maximum total length (TL) of species (based on own observations and FishBase). Coefficients a and b describe allometric regressions between body size and individual functional rates (grazing/bioerosion = a × (TL) b ).
These observations indicate that grazing and bioerosion levels are not driven solely by parrotfish densities, as these were 2-3 times higher in Gaafu Dhaalu compared to Peros Banhos and Salomon ( Figure 4C). Instead, function levels mirror biomass values, which were slightly lower in the Maldives compared to the Chagos Archipelago ( Figure 4D). The main driver of this pattern is the size distribution of parrotfish communities ( Figure 4E), which in Gaafu Dhaalu is heavily skewed towards small individuals (only 1-9% of fish >30 cm TL), while all reefs in the Chagos Archipelago have much higher abundances of bigger fish (32-93% of fish >30 cm TL). Biomass increases non-linearly with fish size, and the presence of large excavators (Peros Banhos) and scrapers (Salomon) in the Chagos Archipelago increased biomass values at these sites disproportionally ( Figure 4F). Given that individual grazing and bioerosion rates also increase non-linearly with fish size, the different size structure of parrotfish assemblages explains most of the among-site differences in reef ecosystem functions.
These observations indicate that grazing and bioerosion levels are not driven solely by parrotfish densities, as these were 2-3 times higher in Gaafu Dhaalu compared to Peros Banhos and Salomon ( Figure 4C). Instead, function levels mirror biomass values, which were slightly lower in the Maldives compared to the Chagos Archipelago ( Figure 4D). The main driver of this pattern is the size distribution of parrotfish communities ( Figure 4E), which in Gaafu Dhaalu is heavily skewed towards small individuals (only 1-9% of fish >30 cm TL), while all reefs in the Chagos Archipelago have much higher abundances of bigger fish (32-93% of fish >30 cm TL). Biomass increases non-linearly with fish size, and the presence of large excavators (Peros Banhos) and scrapers (Salomon) in the Chagos Archipelago increased biomass values at these sites disproportionally ( Figure 4F). Given that individual grazing and bioerosion rates also increase non-linearly with fish size, the different size structure of parrotfish assemblages explains most of the among-site differences in reef ecosystem functions.  Besides size structure of parrotfish assemblages, species composition may influence ecosystem grazing and bioerosion rates. In all three atolls, Ch. sordidus was the most abundant parrotfish (GD: 803, PB: 165, SA: 345 ind. ha −1 ), followed by S. niger (312 ind. ha −1 ) and S. psittacus (118 ind. ha −1 ) in Gaafu Dhaalu, by Ch. strongylocephalus (163 ind. ha −1 ) in Peros Banhos and by S. rubroviolaceus (128 ind. ha −1 ) in Salomon. Despite being classified as an excavator, Ch. sordidus contributed relatively little to total bioerosion, but around 30% to total grazing function in all three atolls ( Figure 5). Further important contributors to grazing at Gaafu Dhaalu were S. niger and S. prasiognathos, the latter however only due to a school of 28 individuals observed at one reef site, which increased mean grazing function at Gaafu Dhaalu (see outlier in Figure 4A). At Peros Banhos, high abundances of large Ch. strongylocephalus, and at Salomon, high abundances of S. rubroviolaceus contributed significantly to total grazing and bioerosion levels ( Figure 5).
Diversity 2020, 12, x FOR PEER REVIEW 10 of 16 Besides size structure of parrotfish assemblages, species composition may influence ecosystem grazing and bioerosion rates. In all three atolls, Ch. sordidus was the most abundant parrotfish (GD: 803, PB: 165, SA: 345 ind. ha −1 ), followed by S. niger (312 ind. ha −1 ) and S. psittacus (118 ind. ha −1 ) in Gaafu Dhaalu, by Ch. strongylocephalus (163 ind. ha −1 ) in Peros Banhos and by S. rubroviolaceus (128 ind. ha −1 ) in Salomon. Despite being classified as an excavator, Ch. sordidus contributed relatively little to total bioerosion, but around 30% to total grazing function in all three atolls ( Figure 5). Further important contributors to grazing at Gaafu Dhaalu were S. niger and S. prasiognathos, the latter however only due to a school of 28 individuals observed at one reef site, which increased mean grazing function at Gaafu Dhaalu (see outlier in Figure 4A). At Peros Banhos, high abundances of large Ch. strongylocephalus, and at Salomon, high abundances of S. rubroviolaceus contributed significantly to total grazing and bioerosion levels ( Figure 5). As biomass and individual grazing and bioerosion rates increase with increasing body size, it is logical to assume that parrotfish biomass at a given site is a good proxy for the magnitude of ecosystem functions. Indeed, both grazing and bioerosion levels could be predicted by transect-level biomass values with a high certainty (R 2 = 0.83 and 0.78, respectively). Across all sites, average (mean ± SE) grazing increased by 66 ± 5 m 2 ha −1 and bioerosion increased by 109 ± 9 kg ha −1 for every kg increase in biomass (Grazing (m 2 m −2 yr −1 ) = 0.0066 × Biomass (kg ha −1 ) + 1.194; Bioerosion (kg m −2 yr −1 ) = 0.0109 × Biomass (kg ha −1 ) − 0.319). However, for a given level of biomass, grazing and bioerosion levels are higher when parrotfish assemblages are dominated by large numbers of small individuals as observed in the Maldives ( Figure 6). For grazing, this seems to be true across the whole range of observed biomass values, whereas for bioerosion, the difference is only true for sites with biomass values >300 kg ha −1 . In the Maldives, this was only the case for 20% of transects. As biomass and individual grazing and bioerosion rates increase with increasing body size, it is logical to assume that parrotfish biomass at a given site is a good proxy for the magnitude of ecosystem functions. Indeed, both grazing and bioerosion levels could be predicted by transect-level biomass values with a high certainty (R 2 = 0.83 and 0.78, respectively). Across all sites, average (mean ± SE) grazing increased by 66 ± 5 m 2 ha −1 and bioerosion increased by 109 ± 9 kg ha −1 for every kg increase in biomass (Grazing (m 2 m −2 yr −1 ) = 0.0066 × Biomass (kg ha −1 ) + 1.194; Bioerosion (kg m −2 yr −1 ) = 0.0109 × Biomass (kg ha −1 ) − 0.319). However, for a given level of biomass, grazing and bioerosion levels are higher when parrotfish assemblages are dominated by large numbers of small individuals as observed in the Maldives ( Figure 6). For grazing, this seems to be true across the whole range of observed biomass values, whereas for bioerosion, the difference is only true for sites with biomass values >300 kg ha −1 . In the Maldives, this was only the case for 20% of transects.

Discussion
We quantified feeding metrics for a whole suite of parrotfishes in the central Indian Ocean, including for species with no prior data on scar frequencies and sizes. Individual grazing and bioerosion rates demonstrate pronounced differences in functional impact among species, even within functional groups. Ch. sordidus, S. frenatus, S. niger, and S. scaber show highest grazing rates for a given body size, and C. ocellatus, Ch. strongylocephalus and S. rubroviolaceus show highest bioerosion rates due to large bite volumes and body sizes. Functional rates of all species showed nonlinear relationships with body size, emphasizing the disproportional effect of large species and individuals. At our study sites, the dominance of large-bodied fish in the Chagos Archipelago was therefore able to compensate for 2-3 times higher parrotfish densities in the Maldives, causing similar community-level grazing and bioerosion across the studied atolls.

Individual Grazing and Bioerosion Rates
The strong relationships of parrotfish grazing and bioerosion rates with body size that we found for all studied species is consistent with previous studies [5,10,16,20,41]. This effect is mainly driven by non-linear increases of scar size and proportion of bites leaving scars with increasing body size. We therefore stress that proportion of bites leaving scars represents an important metric in parrotfish feeding studies, although it has only been sporadically quantified in the past (but see [17,18,20,41,42]). Even including this study, data on this metric is only available for 11 Indo-Pacific parrotfishes, compared to data on scar sizes and bite rates for 24 and 35 species, respectively. Species-specific comparisons with similar studies show that bioerosion rates of large Ch. strongylocephalus seem to be fairly consistent within the central Indian Ocean (400-460 kg ind. −1 yr −1 [19,20]), yet smaller than rates of the sister species Ch. microrhinos on the Great Barrier Reef (up to 1018 kg ind. −1 yr −1 [30]). This suggests potentially significant differences in functional rates of even closely related species or regional differences in behavior or substrate densities. Biorosion rates of Ch. sordidus and sister species Ch. spilurus however are comparable across the Red Sea, Indian Ocean and Great Barrier Reef, despite differences in habitats and environmental conditions (24-55 kg ind. yr −1 [19,20,30,43]). Rates for S. rubroviolaceus, S. frenatus and S. niger in this study are similar to those reported elsewhere in the Indo Pacific (230-380, 44 and 15 kg ind. −1 yr −1 , respectively [17,43]), but much larger than those reported by Yarlett et al. [20] for the Maldives (3.4 kg ind. −1 yr −1 for S. rubroviolaceus, 1.4-1.6 kg ind. −1 yr −1 for smaller Scarus spp.). This indicates that differences in

Discussion
We quantified feeding metrics for a whole suite of parrotfishes in the central Indian Ocean, including for species with no prior data on scar frequencies and sizes. Individual grazing and bioerosion rates demonstrate pronounced differences in functional impact among species, even within functional groups. Ch. sordidus, S. frenatus, S. niger, and S. scaber show highest grazing rates for a given body size, and C. ocellatus, Ch. strongylocephalus and S. rubroviolaceus show highest bioerosion rates due to large bite volumes and body sizes. Functional rates of all species showed non-linear relationships with body size, emphasizing the disproportional effect of large species and individuals. At our study sites, the dominance of large-bodied fish in the Chagos Archipelago was therefore able to compensate for 2-3 times higher parrotfish densities in the Maldives, causing similar community-level grazing and bioerosion across the studied atolls.

Individual Grazing and Bioerosion Rates
The strong relationships of parrotfish grazing and bioerosion rates with body size that we found for all studied species is consistent with previous studies [5,10,16,20,41]. This effect is mainly driven by non-linear increases of scar size and proportion of bites leaving scars with increasing body size. We therefore stress that proportion of bites leaving scars represents an important metric in parrotfish feeding studies, although it has only been sporadically quantified in the past (but see [17,18,20,41,42]). Even including this study, data on this metric is only available for 11 Indo-Pacific parrotfishes, compared to data on scar sizes and bite rates for 24 and 35 species, respectively. Species-specific comparisons with similar studies show that bioerosion rates of large Ch. strongylocephalus seem to be fairly consistent within the central Indian Ocean (400-460 kg ind. −1 yr −1 [19,20]), yet smaller than rates of the sister species Ch. microrhinos on the Great Barrier Reef (up to 1018 kg ind. −1 yr −1 [30]). This suggests potentially significant differences in functional rates of even closely related species or regional differences in behavior or substrate densities. Biorosion rates of Ch. sordidus and sister species Ch. spilurus however are comparable across the Red Sea, Indian Ocean and Great Barrier Reef, despite differences in habitats and environmental conditions (24-55 kg ind. yr −1 [19,20,30,43]). Rates for S. rubroviolaceus, S. frenatus and S. niger in this study are similar to those reported elsewhere in the Indo Pacific (230-380, 44 and 15 kg ind. −1 yr −1 , respectively [17,43]), but much larger than those reported by Yarlett et al. [20] for the Maldives (3.4 kg ind. −1 yr −1 for S. rubroviolaceus, 1.4-1.6 kg ind. −1 yr −1 for smaller Scarus spp.). This indicates that differences in methodology, especially challenges associated with measuring bite scar sizes, may impede direct comparisons among studies (for instance the use of few size classes and very small scar depths in [20]).
For S. scaber, S. tricolor and S. psittacus, some of the smaller but frequently found species on Indian Ocean reefs, this study provides the first annual grazing and bioerosion rates. They indicate that despite being of similar size and morphology, S. scaber contributes substantially to grazing, while S. tricolor left few and S. psittacus hardly any bite scars on the substrate. This may partly be a regional anomaly, as S. psittacus and initial phase S. tricolor were observed to intensively feed on the calcifying green alga Halimeda spp., a behavior that has not been reported previously. This might be a function of exploiting a newly abundant food source, as Halimeda has been proliferating in Maldivian reefs since the 2015-2016 bleaching event, reaching 21 ± 8% cover in 2019 [13]. However, other scraper species have also been found to contribute very little to community-level grazing or bioerosion, as they primarily forage on long, sediment-laden turf algae and rarely leave bite scars [21,42]. Therefore, these findings may support recent work on interspecific resource partitioning, as parrotfish species were found to target different substrates with distinct microautotroph communities [44].
The large variation among individual grazing and bioerosion rates within functional groups emphasizes that the categorization in excavators and scrapers, which was based on behavior, evolutionary division and jaw structure of species [4,31,45], does not directly translate to their functional impacts on reefs. Clearly, this is not a new realization and prior studies have differentiated between large and small excavators [9,11,46] or characterized large individuals of scraper species as occasional excavators [4,17,43]. However, this study demonstrates that large variations in functional rates also exist among small scrapers and highlights the species-specific differences in the rates of change with body size. The feeding metrics and grazing and bioerosion rates provided here thus add significantly to the available data on parrotfish functional impacts and support more accurate quantification of community-level grazing and erosion in the Indian Ocean and beyond.
The approach of merging feeding metric regressions into grazing and bioerosion formulas offers the advantage of estimating ecosystem functions directly from survey data. However, the natural variability in feeding behavior and the propagation of uncertainties from feeding metrics to grazing/bioerosion rates result in relatively large confidence intervals for some species, especially in the larger size classes (Figure 3). This potential variability has to be considered when interpreting community-level ecosystem functions and could be reduced by higher replication of feeding observations. As wave exposure may change feeding rates of herbivores, additional observations might be required at high energy reef sites. However, most parrotfish species show consistent feeding rates across exposure gradients [24,25]. In higher latitudes, seasonal variations in feeding rates should also be taken into account [22,23].

Community Grazing and Bioerosion Levels
After the 2015-2016 bleaching event, parrotfish abundances and body sizes in the Maldives and the Chagos Archipelago increased considerably [13,29], likely as a result of increased food availability in the form of photoautotrophic microbes in dead reef substrate [2]. The same pattern has been observed for Moorea after large-scale coral mortality following an outbreak of crown-of-thorns starfish [47]. These changes in parrotfish assemblages naturally have implications on ecosystem functions and can facilitate the return to coral dominated states after disturbance by keeping algae in a closely cropped state [47,48]. In support of this we observed no long, sediment laden turf in either location in the present study, and there was evidence of early coral recovery, especially in the Chagos Archipelago. High community-level bioerosion was also apparent as the reef structure was actively degrading and sand accumulated on the reef slopes.
Despite 2-3 times higher parrotfish densities in the Maldives than in the Chagos Archipelago, grazing and bioerosion levels were similar at both locations, because the higher proportion of large parrotfish in the Chagos Archipelago compensated for lower densities. Reasons for the skewed distribution towards small parrotfish in the Maldives are likely habitat characteristics of the studied reefs, which can influence parrotfish assemblages [9,49], and may also reflect the lack of top predators which predominantly target small parrotfish [50]. Higher exposure to human influence in the Maldives compared to the remote Chagos Archipelago may also be a factor [11], although the pressure of reef fishery in the Maldives is generally considered to be low [51].
Despite confirming the finding that large individuals and species contribute disproportionally to ecosystem functions (e.g. reviewed in [3,21]), we observed that for a given level of biomass, the dominance of small-bodied fishes in the Maldives enhanced community grazing and bioerosion levels compared to reefs in the Chagos Archipelago, supporting recent observations by Robinson et al. [10]. This suggests that in regions where the mean size of parrotfish is reduced due to targeted fishing of large reef fish [11,52], an increase in smaller individuals can potentially offset the associated decrease in ecosystem functions. This has been observed in the Caribbean, where very small parrotfish exerted grazing pressure strong enough to pre-empt macroalgal dominance following coral mortality [53]. However, this is only true if recruitment of juveniles remains high, as 2-3 times higher abundances were necessary to compensate the effects of smaller body sizes in this study. Furthermore, an analysis of Caribbean-wide patterns in herbivore and algal communities showed that parrotfish size rather than biomass drives algal cover across sites, emphasizing the importance of large individuals [54]. For the Great Barrier Reef, it has been shown that the disproportionally high bioerosion and coral predation impacts by Bolbometopon muricatum cannot be replaced by any other species [9,11]. To a lesser extent this may be also be true for other species.
Our study shows that, in the central Indian Ocean, key species for community-level grazing are Ch. sordidus and S. niger due to their high abundances, and Ch. strongylocephalus and C. ocellatus due to their large body size. The latter two species are also the most critical ones to provide bioerosion functions on reefs. Where present, S. prasiognathos and S. rubroviolaceus significantly contribute to grazing, and in the latter case also to bioerosion levels. These species often occur in big schools or small groups that move across the reef, so it is difficult to scale their impact on a representative unit area and emphasizes the need for sufficient replication of fish censuses within sites and over time.
Comparing community grazing and bioerosion levels in the Maldives and the Chagos Archipelago to other coral reef regions, average annual impacts are around twice those determined for reefs in the Red Sea and Arabian Sea (1.5-2.1 m 2 m −2 yr −1 and 0.4-1.6 kg m −2 yr −1 ), where parrotfish biomass is considerably lower (150-200 kg ha −1 [55]). Compared to the GBR, grazing and bioerosion levels in the Central Indian Ocean are similar to those on mid shelf reefs (5.2-6.3 m 2 m −2 yr −1 and 5.2-8.4 kg m −2 yr −1 ). On inner shelf reefs, grazing is higher (8.5-11.0 m 2 m −2 yr −1 ), as abundances of scrapers are very high. Bioerosion is highest on outer shelf reefs (23-32 kg m −2 yr −1 ), due to the occurrence of B. muricatum [9].
The variance in community grazing and bioerosion levels across our study sites and other reef regions demonstrates that these ecosystem functions are not only determined by total parrotfish biomass, but also by species composition and population size structure. Quantifying how these functions vary with species density and biomass has considerable potential not only for enhancing our ability to understand the functional impacts of fishing effort and ecosystem decline, but also, as a result, for improving management efforts. In particular, there should be a strong emphasis on assessing which parrotfish species and size classes drive most of the community grazing and bioerosion within regions, and ideally some recognition of the biomass levels needed to sustain these functions at site-appropriate levels.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/10/379/s1. Figure S1: Body sizes (A) and bite rates (B) of parrotfish species comparing initial (I) and terminal (T) life phases. Boxes are 25th and 75th percentiles with median line; whiskers and points extend to smallest and largest values. Results of Welch two sample t-tests comparing life phases are significant for all species in (A) and as indicated in (B): '***' p<0.001, '**' p<0.01, '*' p<0.05" "ns" not significant; Table S1: Analysis of deviance and F-tests for generalized linear models (GLM), testing the influence of factors "species", "size" and "life phase" on feeding metrics bite rate (bpm), proportion of bites leaving scars (%scars), scar area (area) and scar volume (vol). R 2 -pseudo states the proportion of variation explained by the full model (R 2 -pseudo = (Dev -Resid. Dev)/Resid. Dev) and the absolute contribution of each predictor variable; Table S2: Individual feeding metrics for a range of parrotfish species and size classes. Coefficients a and b describe regressions between body size (TL) and individual feeding metrics and sign. reports the significance level of the slope b. Figure S2: Monte-Carlo simulations to propagate uncertainties of feeding metric models into individual parrotfish grazing (m 2 ind. −1 yr −1 ) and bioerosion rates (kg ind. −1 yr −1 ). Grey lines represent 1000 possible regression lines calculated from random coefficients within the error distribution of each feeding metric model following formulas (1) and (2). Red lines represent mean models over all calculated coefficients. Blue lines represent the models using mean coefficients of feeding metric regressions