Client fish traits underlying variation in service quality in a marine cleaning mutualism

https://doi.org/10.1016/j.anbehav.2021.03.005 0003-3472/© 2021 The Authors. Published by Elsevie license (http://creativecommons.org/licenses/by/4.0/) Game-theoretical models help us understand how and when cooperation can evolve and persist. However, current models fall short of explaining the striking amount of variation in cooperation levels that we observe in nature, even within a system. For example, an animal's ability to choose partners with which to interact can explain the maintenance of cooperative interactions, but not variation in cooperation levels among individuals with similar partner choice options. Here, we explored how partner choice options and other key partner traits predict the quality of the service provided by cleaner fish, Labroides dimidiatus (the ‘chosen’ partner) to their client reef fishes (the ‘choosy’ partner). In this marine cleaning mutualism, cleaner fish sometimes cheat and eat mucus from their clients rather than cooperate by removing ectoparasites. We examined 13 client fish species and assessed the influence of nine variables on cooperation levels by cleaner fish (i.e. the frequency of cheating events) using information theory. Variables characterizing clients and the interactions themselves were used as proxies for a client's nutritional value, partner choice options, ability to punish a cheating partner and a cleaner's temptation to cheat its client (some of these traits were characterized by more than one variable). We found that six of the nine variables were equally important in explaining variation in client jolts: the duration of an interaction, and traits of client species including their body size (a proxy for the quality of a food patch), gnathiid ectoparasite load (nutritional value), partner choice options, agility (ability to control an interaction via avoidance or termination) and mucus characteristics (a measure of the cleaners' temptation to cheat). These findings suggest that numerous factors influence variation in cooperation levels, showing a parallel with sexual selection theory, where female mate choice does not depend on a single factor such as the operational sex ratio, but also on other key variables including a male's quality relative to other males, and a female's ability to avoid male coercion. © 2021 The Authors. Published by Elsevier Ltd on behalf of The Association for the Study of Animal Behaviour. This is an open access article under the CC BY license (http://creativecommons.org/licenses/ by/4.0/).

Cooperation and mutualism (interspecific cooperation) are interactions in which all participants gain higher payoffs than if they had acted alone (Bshary & Bergmüller, 2008;Lehmann & Keller, 2006). Biological market theory (BMT) is a form of game theory that aims to understand how partners share the gains of cooperation. According to BMT, many instances of cooperation and mutualism in nonhuman animals can be described as an exchange of goods or services between traders in the context of a 'biological market' (No€ e, Van Schaik, & van Hooff, 1991;No€ e & Hammerstein, 1995). Traders with goods and services in short supply (or high demand) can exert partner choice and choose from a pool of potential partners with which to interact (No€ e et al., 1991). Partner choice options lead to higher exchange rates in the market because of competition among the 'chosen' trader class through outbidding (No€ e & Hammerstein, 1994, 1995. A key example is delayed plumage maturation of males in some bird species such as the purple martin, Progne subis, and the lazuli bunting, Passerina amoena. In these species, young males (the 'chosen' partners) are junior partners to older males (the 'choosy' partners) and are allowed to co-breed in these older males' territories (Greene et al., 2000;Morton, Forman, & Braun, 1990). Females preferentially mate with males that show bright plumage. Therefore, old males prefer junior partners with a dull plumage to increase their chance of fertilizing the eggs of a junior partner's female and to reduce the risk of cuckoldry by their own female partner. BMT predicts two possible evolutionary outcomes in this situation: if young males are abundant and compete as co-breeders for access to the territories of older males (i.e. there is a high demand for territories), selection on outbidding competitors is strong and young males should uniformly display delayed plumage coloration to increase their chance of gaining access to a territory. In contrast, when young males are not abundant (i.e. there is a low demand for territories), they should all undergo rapid plumage maturation and increase their attractiveness to females since the older dominant males will prefer a brightly coloured co-breeder over breeding alone (No€ e & Hammerstein, 1994).
Partner choice is a powerful control mechanism that can stabilize cooperative behaviour by allowing partner switching in response to the chosen partner offering less than the market value during an interaction (i.e. cheating; Ferriere, Bronstein, Rinaldi, Law, & Gauduchon, 2002;Johnstone & Bshary, 2008;Izquierdo, Izquierdo, & Vega-Redondo, 2010;Wubs, Bshary, & Lehmann, 2016). Numerous different systems that have been studied in a biological market context, from bacteria to primates, have supported the predictions of BMT (sensu Noe & Voelkl, 2013). However, more nuanced outcomes are often predicted by sexual selection theory, which inspired the development of BMT (No€ e, 2017;No€ e et al., 1991;No€ e & Hammerstein, 1995). For instance, while changes in the operational sex ratio have long been recognized as a primary determinant of female mate choice (Emlen & Oring, 1977), additional factors such as a female's own quality (e.g. Kvarnemo & Simmons, 1999;Reinhold, Kurtz, & Engqvist, 2002) and ability to evade coercive mating by males (Clutton-Brock & Parker, 1995) will also affect her ability to choose with which partner(s) to mate. Thus, quantifying only partner choice options by the choosing partner (i.e. the female) would yield an incomplete picture of a mating market. Similarly, factors other than supply and demand may affect how partners share the benefits derived from cooperation or mutualism in a biological market context.
We studied a marine cleaning mutualism involving the cleaner wrasse, Labroides dimidiatus (hereafter 'cleaners') and its 'client' reef fishes (hereafter 'clients'). In this system, clients can be considered as the 'choosy' partner class and cleaners as the 'chosen' partner class. Given previous insights from sexual selection theory (see No€ e, 2017), we aimed to evaluate how a diversity of client traits, including their partner choice options, quality as a food patch and ability to control the course of an interaction, affect the outcome of cooperative exchanges in this mutualism (i.e. the payoff distribution). Clients visit cleaners at their territories (cleaning stations) to have their ectoparasites removed (reviewed in Cote, 2000): cleaners gain access to food (Grutter, 1996(Grutter, , 1999 and clients gain health benefits from parasite removal (Binning et al., 2018;Clague et al., 2011;Demair e et al., 2020;Ros, Nusbaumer, Triki, Grutter, & Bshary, 2020;Triki, Grutter, Bshary, & Ros, 2016). Cleanereclient interactions have frequently been used to test the predictions of BMT due to this exchange of goods and services (Bshary & Grutter, 2002b;. In addition to clients having different partner choice options, this mutualism is characterized by conflicts of interest because cleaners tend to prefer eating their clients' protective mucus over their ectoparasites (Grutter & Bshary, 2003). Mucus feeding constitutes cheating and clients must force cleaners to feed against their preference to receive a good service. To achieve this, clients respond to cheating by cleaners via various partner control mechanisms that reduce the cheater's payoff. Predatory clients (i.e. piscivores) could respond to cheating cleaner by trying to eat the cleaner, an option that seemingly leads to almost unconditional cooperative behaviour by cleaners (Bshary, 2001). Nonpredatory 'visitor' clients with large home ranges and access to several cleaning stations use their partner choice options to switch partners when cheated; in contrast, nonpredatory 'resident' clients with small home ranges and access to a single cleaning station tend to punish cheating cleaners by chasing them aggressively (Bshary & Grutter, 2002a). Such responses by clients cause cleaners to behave more cooperatively by feeding against their preference (Bshary & Grutter, 2005). Previous work has shown that cleaners prioritize visitor clients over resident clients to avoid missed foraging opportunities (Bshary, 2001). However, large unexplained variation exists in the quality of the cleaning services received by both resident and visitor clients (Figure 7.4, p 159 in Bshary, 2001). Variation among visitor clients is particularly intriguing as it is not linked to variation in client body size (Bshary, 2001), a proxy for the abundance of ectoparasites on fish hosts (Grutter, 1994(Grutter, , 1995. Game-theoretical models (Johnstone & Bshary, 2002 and laboratory experiments (Gingins, Werminghausen, Johnstone, Grutter, & Bshary, 2013) point to two factors that can affect the quality of services exchanged by cleaners beyond partner choice. First, the ability of a client to control the course of an interaction, either by avoiding or by terminating an exchange, should affect service quality irrespective of the client's partner choice options. Avoidance is partly a function of client home range size, as a large home range enables clients to stay out of reach of cleaners they want to avoid (by remaining outside their territory). The extent of a client's control of an ongoing interaction is linked to its swimming ability (speed and agility) relative to that of cleaners. If cleaners outperform clients in swimming speed and/or agility, they can prolong and/or initiate interactions against a client's will. The quality of cleaning services is expected to be low if cleaners can coerce their clients to engage in an interaction (i.e. when clients do not have the 'power' to avoid unwanted interactions; Bowles & Hammerstein, 2003). Coercion violates the fundamental assumption of classic BMT that exchanges are voluntary (No€ e & Hammerstein, 1995). Second, a cleaner's temptation to cheat its client will differ between species due to interspecific differences in ectoparasite load and mucus quantity/quality, with size positively affecting ectoparasite load (Arnal, Côt e, & Morand, 2001;Eckes, Dove, Siebeck, & Grutter, 2015;Gorlick, 1980;Grutter, 1994Grutter, , 1995Grutter & Bshary, 2004). As such, large clients tend to represent larger food patches and be more attractive to cleaners than smaller clients. Importantly, however, a high parasite load should encourage cooperation by cleaners but a high quantity or quality of mucus should increase a cleaner's temptation to cheat (Gingins et al., 2013;Grutter & Bshary, 2004;Johnstone & Bshary, 2002).
Here, we studied 13 client fish species (nonpredator residents and visitors) and measured interspecific differences in (1) home range size as a correlate of partner choice options, (2) fast start escape performance and agility (see Supplementary Material Video S1) as correlates of a client's ability to control the course of an interaction and (3) body surface area, parasite load and mucus characteristics as proxies of a client's quality as a food patch. We then examined how these traits relate to cheating by cleaners, which can be observed as client jolts in the wild. Jolts are manifested as an abrupt twitch of the client's body following contact with the cleaner's mouth and are indicative of cleaners eating mucus rather than ectoparasites (Bshary & Grutter, 2002a). Based on the tenets of BMT, we predicted that clients with greater partner choice options would benefit from better service quality (i.e. fewer jolts). In addition, we predicted that fast and agile species would receive better cleaning services than slower, less agile species with restricted abilities to avoid or terminate interactions. Finally, we predicted that cleaners would cheat more on client species with low parasite loads and a high quantity and/or quality of mucus.

METHODS
The study was carried out between July and October 2014 and July and September 2015 at the Lizard Island Research Station (LIRS; 14 40 0 S, 145 28 0 E) on the Great Barrier Reef, Queensland, Australia. Behavioural observations and animal collections were conducted on reefs surrounding Lizard Island (see Appendix Fig. A1). Thirteen client species (Table 1) were chosen based on body size, frequency of interactions with cleaners and estimates of variation in service quality from existing data from the Red Sea (Bshary, 2001(Bshary, , 2002, French Polynesia (Oates, Manica, & Bshary, 2010) and Lizard Island .

Quality of Cleaning Service
Field observations were carried out to determine the quality of the cleaning service provided by L. dimidiatus to the client species listed in Table 1. Observations were carried out on snorkel or shallow SCUBA diving (< 9 m depth) between 0700 and 1630 hours. All species were observed throughout the day to avoid biasing observations. Each cleaning station was observed only once, for a maximum of 20 min. The observers (M.J., D.G.R.) remained at a minimum of 2.5 m from the cleaner and recorded all cleaning events involving the species of interest using a pencil and underwater slate. If a cleaner repeatedly interacted with the same client, we recorded only the first interaction. Given the difficulty of uniquely identifying individuals of abundant species, some clients might have been double counted if they exited and re-entered the observer's field of view and interactions with the cleaner were separated by several minutes. For each cleanereclient interaction, we recorded five variables: (1) the client species, (2) an estimate of the client's total length (TL), (3) the duration of the interaction (s), (4) the number of jolts performed by the client and (5) how the interaction ended (the client swims away from the cleaner; the client jolts and flees; the client jolts and chases the cleaner; the cleaner switches partner and starts cleaning another client). Jolts indicate (i.e. are correlated with) instances when cleaners cheat a client and are manifested as an abrupt twitch of the client's body following contact with the cleaner's mouth (Bshary & Grutter, 2002a). A cleaning interaction started with the cleaner visually inspecting the client's body at close range (<10 cm) and ended when the cleaner switched partner or the client and cleaner became separated by a minimum of 30 cm for more than 2 s (Adam, 2010;Soares, Bshary, Cardoso, & Côt e, 2008). We monitored the duration of interactions using a waterproof stopwatch, wrote observations with a pencil on underwater paper and practised estimating client TL underwater using landmarks on the reef and cut PVC pipes of different lengths (range 3e60 cm).

Client Swimming Performance
Fish were captured with a barrier net (4.5 x 1.5 m; 1 cm stretched mesh) and hand nets on SCUBA and transported to the LIRS aquarium facilities in 20-litre buckets, within 1 h of capture (see Table A1 for sample sizes). Captured fishes were placed in large holding aquaria (1.1 x 0.6 m and 0.5 m deep) with sea water pumped directly from the reef. Shelters were provided in the form of cut PVC pipes. The water temperature was 24.1 ± 0.8 C (mean ± SD) and fishes were exposed to a natural photoperiod of approximately 12 h for 1e3 days prior to the experiments. Fishes were fed daily with prawn or fish flakes but fasted for 24 h prior to the experiments. Each fish's centre of mass (CoM) was marked dorsally with a piece of reflective tape placed on either side of the dorsal fin (Lefrançois, Shingles, & Domenici, 2005). The position of the CoM was determined for one euthanized fish per species, using cooled specimens; a sharp metal pin was inserted along the body midline until the point of balance was identified (Domenici, Standen, & Levine, 2004;Gotanda, Reardon, Murphy, & Chapman, 2012).
We used two large, opaque circular tanks (diameter ¼ 1.10 m, height ¼ 0.40 m) with flow-through sea water (6 litres/min) to maintain a stable temperature (24.1 ± 0.8 C, mean ± SD) and oxygen levels (>90% saturation) throughout the experiments (see Appendix Fig. A2). The water depth in the tanks was set to 20 cm to minimize movement by the fish in the vertical plane but allow full extension of their anal and dorsal fins. Floodlighting was supplied by two 500 W work lights placed at opposite ends of the tank. Escape responses were elicited by mechanoacoustic stimulation. Two weighted plastic vials with a tapered tip (height 12.0 cm, diameter 3.0 cm, weight 165 g) were suspended 70 cm above the water surface on either side of the tank, 13 cm from the edge of the tank (see Appendix Fig. A2; Turesson, Satta, & Domenici, 2009), and were released by remotely switching off an electromagnet using a switch connected by a cable. Using two stimuli per experimental tank facilitated standardizing the orientation of the fish relative to the stimulus. To avoid visual stimulation before it breached the water surface, the stimulus fell inside an opaque PVC cylinder (outer diameter 10.3 cm) suspended 1 cm above the water surface (see Appendix Fig. A2; Lefrançois et al., 2005). Escape responses were filmed at 240 Hz with a high-speed camera (GoPro Hero 3þ black; GoPro, San Mateo, CA, U.S.A.) placed above the centre of the experimental arena (see Supplementary Material Video S2). The camera was custom fitted with a commercial lens to avoid image distortion (4.14 mm f/3.0 86 HFOV 5MP GP41430; Peau Productions Inc, San Diego, CA, U.S.A.) and could be controlled remotely via WiFi using the GoPro App on an Ipad. The experimenter was shielded by a black tarp, behind which the camera and electromagnet were controlled.
Supplementary video related to this article can be found at https://doi.org/10.1016/j.anbehav.2021.03.005 A single fish was placed in the experimental arena and left undisturbed for a minimum of 30 min prior to testing. Fish were tested using both a lateral and a frontal stimulation: we used a lateral stimulation (87.8 ± 29.8 , mean ± SD standardized angle of the body relative to the stimulus) to obtain measures of escape distance (D esc ), maximum speed (U max ) and maximum acceleration (A max ); we used a frontal stimulation (19.8 ± 12.8, mean ± SD) to measure a fish's turning rate, an indicator of its agility. Fast starts were elicited at a fixed distance from the stimulus of approximately The table shows the number of cleaning interactions observed (N), the total length (TL) of clients (mean ± SD) and the duration of interactions (mean ± SD).
25 cm (25.44 ± 5.6 cm, mean ± SD). Some individuals and/or species frequently swam to the centre of the arena, whereas others tended to remain near the edges. In the latter case, we gently moved a ring of PVC pipe attached to strings or temporarily activated a bubble curtain along the walls of the arena, motivating fish to move away from the edges. If a fish moved and considerably changed its position immediately prior to the stimulation, another trial was conducted until a satisfactory response was obtained. We waited a minimum of 15 min between trials, a period deemed appropriate to minimize fatigue and habituation to the stimulus (Jornod & Roche, 2015). Videos of escape responses were analysed frame by frame in ImageJ (NIH, Bethesda, MD, U.S.A.), using the MTrackJ plugin to manually track the CoM (Meijering, Dzyubachyk, & Smal, 2012;Roche, 2021). Stage 1 of the escape response was defined as the time between the first movement of the head and the reversal of the angular motion of the head; stage 2 began with the end of stage 1 and ended with the straightening of the fish's body (Domenici, 2011;Domenici & Batty, 1997;Domenici & Blake, 1997). We measured the following variables: (1) stage 1 turning rate (calculated as the angle between the segment joining the CoM and the tip of the snout at the beginning and end of stage 1 divided by the duration of stage 1); and distanceetime variables including (2) (Lanczos, 1956) was then applied to obtain smoothed values of U max and A max (Lefrançois et al., 2005;Marras, Killen, Claireaux, Domenici, & McKenzie, 2011).

Client mucus Quantity and Nutritional Quality
We followed Gorlick (1980) to assess mucus quantity and quality for three to six fish per client species (see also . These fish were previously used to assess swimming performance (see Results for exact sample size for each species). Fish were euthanized by cervical dislocation and pithing of the brain, submerged in a freshwater bath for 5 s to remove excess sea water and then for 60 s in a 50 C water bath, inducing the mucus to denature and coagulate (Gorlick, 1980). Each fish was then suspended by the jaw on a wooden stand over a beaker and the mucus was then gently scraped off the epidermis with a clean microscope slide. The body of the fish was squirted with warm water to prevent drying. Mucus samples were then freeze dried for 30 h in a Virtis BenchTop Freeze Dryer (SP Industries Inc, Warminster, PA, U.S.A.) and weighed to the nearest 1 mg. To measure body surface area, we covered the fish's body with aluminium foil (excluding the fins), conforming to the body curvature, and cut out the outline. Pectoral fins were dissected at the base of the fin, spread out and pinned onto a foam sheet (Binning & Roche, 2015). Pictures of the pectoral fins, the flattened aluminium foil, and extended dorsal, anal and pelvic fins were used to measure total body surface area in ImageJ (see Supplementary Material Video S2). All pictures contained a scale.
We analysed mucus protein content with a Bradford protein assay (Bradford, 1976). We extracted proteins by adding 100 ml of cold Rensink buffer (50 mM Tris-HCl, pH 7.5, 100 mM NaCl, 0.5% (v/ v) Triton X-100, 10 mM bmercapto ethanol, 1 mM phenylmethylsulphonyl fluoride) to 1 mg of dried, pulverized mucus. We then added 200 ml of Bradford reagent and 796 ml of ddH2O to 4 ml of the supernatant from the centrifuged sample. After 5 min of incubation at room temperature, absorbance was measured at 595 nm with a spectrophotometer (Amersham Biosciences, Little Chalfont, U.K., Ultrospec 3100 Pro). Bovine serum albumin was used to construct a standard curve. Duplicates were obtained for each sample. Caloric content was estimated by CHN analysis, which provides percentages of carbon, nitrogen and ash Gorlick, 1980). We approximated mucus caloric content as: calories/g dry weight ¼ -227 þ 152 (% carbon) following Platt, Brawn, and Irwin (1969). Direct calorimetric measurement was not possible due to the low weights of dried mucus samples obtained from fishes.

Client Ectoparasite Load
We examined the ectoparasite load of five individuals per species following Sikkel, Cheney, and Côt e (2004). Five fish per species is deemed suitable to obtain a reliable estimate of ectoparasite abundance at a given location and time Grutter, 1994). Fish were collected on SCUBA using a barrier net and hand nets, and immediately transferred to individual, sealed 10-litre plastic bags filled with sea water. Fish were brought back to the boat within 15 min of collection and 1 ml of a 10% clove oil solution was added to the plastic bag to induce sedation. We determined that sedation was achieved 5 s after a fish exhibited loss of equilibrium (usually within 2 min); fish continued to ventilate their gills throughout. Each fish was then transferred to a freshwater bath for 5 min and the gills and body surface were then thoroughly rinsed for 1 min with freshwater using a wash bottle. A battery powered pump and air stone were used to aerate the water at each of these steps. Fish were photographed to measure body surface area and released at their site of capture after a 30 min recovery period in aerated sea water. The water from the plastic bag and freshwater bath was filtered with a plankton net (80 mm mesh size) to collect parasites. The samples were then brought back to the laboratory, where parasites were identified and counted with a binocular microscope. Total body surface area was calculated using a speciesspecific linear regression between body surface area (excluding the fins) and total body surface (see Client mucus quantity and nutritional quality). Since parasitic gnathiid isopods make up approximately 70% of the total number of items in the diet of L. dimidiatus (Grutter, 2000), we classified ectoparasites into two groups: gnathiid isopods and other ectoparasites. Other ectoparasites included copepod and monogenean parasites which are known to be consumed by L. dimidiatus (Grutter, 1997;Grutter & Bshary, 2003). We calculated a host's ectoparasite load as the number of (1) gnathiid isopods and (2) other ectoparasites per square centimetre of body surface area.

Client Partner Choice Options
We used home range size estimates from the literature (Table A2) as a correlate of partner choice options for our client species of interest (Table 1). Cleaner densities around Lizard Island typically vary between 0.33 and 2.33 individuals per 100 m 2 (Triki, Levorato, McNeely, Marshall, & Bshary, 2019). Taking into consideration that cleaners occasionally occur in pairs sharing a cleaning station, we sorted client species into three categories. Client species with a home range smaller than 100 m 2 (actual range 3e87 m 2 ) were considered as having few partner choice options; clients with a home range between 100 and 300 m 2 (actual range 100e245 m 2 ), as having intermediate partner choice options; and clients with a home range larger than 1000 m 2 (actual range 2300e15 300 m 2 ), as having many partner choice options (Table A2). This classification of partner choice options based on distinct home range sizes reflects categorizations of clients into 'residents' (small home range size providing access to a single cleaning station) and 'visitors' (medium or large home range size providing access to more than one cleaning station) used in several previous studies (Bshary & Grutter, 2002b;Bshary & Schaffer, 2002;.

Ethical Note
All field collections and experiments for this study were conducted under permits from the Great Barrier Reef Marine Park Authority (G14/37047.1) with approval from the Queensland Government (DAFF) Animal Ethics Committee (CA 2014/06/780).

Statistical Analysis
Given the exploratory nature of our analysis, we used the information-theoretic approach (model selection) to identify which variables best explained variation in jolts by client fishes (Burnham & Anderson, 2002;Grueber, Nakagawa, Laws, & Jamieson, 2011;Johnson & Omland, 2004). Collinearity among continuous predictors was assessed with the R function 'corrplot' (Fig. A3). We used two principal component analyses (R function 'PCA') to condense collinear predictors: one PCA was carried out for variables related to fast start performance (D esc , U max , A max ) and one for two of the three variables related to mucus properties (mucus quantity and caloric content). In both cases, the first principal component (PC1) was positively correlated with the variables included in the PCA and accounted for over 85% of the variance (Table A3).
We specified a full model as a generalized linear model (GLM; R function glm) with a Poisson distribution of error terms and nine predictor variables: duration of the interaction, client size (TL), fast start escape performance (PC fast start), turning rate, mucus quantity and caloric content (PC mucus), mucus protein content, number of gnathiid isopods/cm 2 , number of other ectoparasites/ cm 2 and home range size as a proxy for partner choice options (low, medium, high). The ratio of residual deviance over degrees of freedom (0.87) was close to 1, suggesting no over/underdispersion. Since continuous predictors were measured on different scales, predictors were z-standardized (i.e. mean ¼ 0, SD ¼ 1) prior to the analysis. Model assumptions were checked with plots of residuals versus linear predictors, a qqplot of residuals and plots of influence measures. Two observations with extreme interaction duration times greater than 150 s were excluded from the data set to satisfy model assumptions. We examined nonlinear relationships using pairs plots (R function pairs); quadratic terms were included for two predictors: interaction duration and client total length. We used the Akaike information criterion (AIC) to select the best candidate models with the R function 'glmulti' (Calcagno & de Mazancourt, 2010). Model averaging was performed with the R function 'model.avg' (Barton, 2015) on all models within the 95% confidence set (i.e. summed normalized Akaike weight values (w im ) of the best models ¼ 0.95). We calculated the percentage deviance explained to evaluate the model's goodness-of-fit with the R function Dsquared (Barbosa, Brown, Jimenez-Valverde, & Real, 2015). All analyses were done in R 3.1.2 (R Development Core Team, 2014).

Limitations
We specified the response variable as the number of client jolts per cleanereclient interaction (our unit of replication) rather than the number of client jolts per 100 s as has often been done in previous studies (Bshary, 2001;Bshary & Grutter, 2002a), allowing us to examine the effect of interaction duration on cleaning service quality. Nevertheless, our chosen response variable strongly correlates with the previously used measure of service quality given a strong relationship between the proportion of interactions without jolts and the number of jolts per 100 s of interaction for each species (see Appendix Fig. A4).
Cleaner fish identity was not included as a random effect in the model because this variable was unfortunately not recorded for most of the behavioural observations in the wild. Interaction terms were not included in the model for two reasons: first, we had no specific a priori predictions and the model would have had to include 36 two-way interactions given nine main effects; second, evidence suggests that estimating interaction terms requires 16 times the sample size needed to estimate main effects with adequate power (Gelman, 2018). We could not account for phylogenetic relatedness among client species due to inaccuracies when fewer than 20 taxa are included in phylogenetic analyses (Blomberg, Garland, Ives, & Crespi, 2003;Boettiger, Coop, & Ralph, 2012). Instead, we controlled for phylogeny by examining a single client species per fish family (Table 1).

RESULTS
We observed 1056 cleaning interactions over ca. 90 h of underwater observation. A minimum of 49 cleaning events were recorded per species except for A. nigropunctatus and L. obsoletus for which we observed 12 and 17 cleaning events, respectively (Table 1).
There was considerable variation in the quality of the cleaning services provided by cleaners to the 13 client reef fish species (Table 1). The percentage of interactions with at least one jolt, a correlate of jolts per time unit (see Appendix Fig. A4), ranged from 6.1% to 100% across species (Fig. 1). Similarly, the percentage of interactions comprising at least one jolt and terminated by a client chase ranged from 0 to 66% depending on the species (see Appendix Fig. A5). We also found notable variation in the explanatory variables we examined. For example, the mean duration of cleaning interactions per client species ranged from 1.3 to 31.1 s (Table 1); the client species that performed best in the fast start escape response experiments covered a distance (D esc ) four times greater than the worst performing species (Table A1); the turning rate (agility) of the best performing client species was over four times that of the worst (Table A4); mucus quantity varied from 0.12 to 0.89 g/cm 2 among client species (Table A5); highly parasitized clients harboured over 30 gnathiid isopods on average, whereas some had less than one (Table A6); and the estimated home range of clients ranged from 3 to 15 300 m 2 (Table A2).
Six models had support in explaining the number of jolts per cleaning interaction based on AIC scores and were included in the 95% confidence set of models (Table 2). Of the nine predictors included in the full model, six appeared in all models included in the 95% confidence set. The model with the highest support (w im ¼ 0.30) explained 27.7% of the deviance in jolt number (Table 2).
Six predictors had 95% confidence intervals that did not overlap zero and were equally important in determining the quality of services received by clients ( Table 2). The duration of a cleaning interaction was positively related to the number of jolts up to a plateau at 50 s ( Fig. 2a and b). Client size was also nonlinearly related to the number of jolts performed by a client: jolts increased with size for smaller clients but decreased with size for larger clients ( Fig. 2c and d). Jolts were positively related to mucus quantity and caloric content ( Fig. 2e and f), but negatively related to client turning rate and the number of gnathiid isopods infecting a client (Fig. 2gej). Finally, clients with no partner choice options ('resident' clients with a small home range and access to a single cleaner) received a cleaning service of lower quality (i.e. more jolts) than clients with access to more than one cleaner ('visitor' clients with a medium or large home range and access to more than one cleaner; 0 25 50 75 100

C . b a r o n e s s a H . m e l a p t e r u s L . o b s o l e t u s A . c u r a c a o S . d o l i a t u s S . c h r y s o p t e r u s L . q u i n q u e l i n e a t u s P . m u l t i f a s c i a t u s P . s e x s t r i a t u s C . s t r i a t u s S . b i l i n e a t a S . f l a v i p e c t o r a l i s A . n i g r o p u n c t a t u s
Percentage of cleaning interactions  Figure 1. The percentage of interactions with zero, one or more than one jolts for 13 client reef fish species from field observations (y axis on left). The data are ordered from highest to lowest percentage of jolts. The red horizontal bars represent the overall number of jolts per 100 s for each species (y axis on right). The number of observations per species is given under each bar. Data for the six explanatory variables deemed important in influencing service quality are shown above the graph. Interaction duration (s) and client size (TL; cm) are means ± SE. Values for mucus properties (PC mucus ¼ principal component positively related to mucus quantity and caloric content), turning rate ( /s) and gnathiid load are divided into five categories (very low, low, medium, high, very high) based on percentiles (0.2, 0.4, 0.6, 0.8). Colour gradient: dark red boxes indicate a negative association with the quality of cleaning services (i.e. more jolts) as predicted by model-averaged parameter estimates (Table 2). The model-averaged parameter estimate (b), unconditional standard error (SE), 95% confidence interval (95% CI) and normalized Akaike weight (w ip ) for each predictor are shown. Also indicated are models in which predictors were included (), the number of parameters included in each model (K), the percentage deviance explained, as well as each model's AIC, DAIC and normalized Akaike weight (w im ). Predictors are in order of importance (w ip ); those for which the 95% CI does not overlap zero are indicated in bold.
All models include a constant (intercept). TL ¼ total length. PC mucus ¼ principal component positively related to mucus quantity and caloric content. PC fast start ¼ principal component positively related to escape distance (D esc ), maximum velocity (U max ) and maximum acceleration (A max ).   Fig. 2k and l). Conversely, the 95% confidence interval of estimates for mucus protein content, ectoparasites other than gnathiids, and client fast start escape performance overlapped zero (Table 2). These predictors had 1.16, 1.69 and 3.57 times less support than previously listed variables in explaining variation in the number of jolts performed by clients (Table 2).

DISCUSSION
Our findings suggest that multiple factors, including traders' partner choice options, temptation to cheat and ability to control the course of an interaction, all influenced the quality of the cleaning services exchanged in our study system. Notably, the increase in cleaning service quality we observed with increasing home range size estimates, and hence correlated partner choice options, supports the predictions of classic BMT (No€ e & Hammerstein, 1994, 1995. Interestingly, our data suggest that, in contrast to having no partner choice options, having the choice to interact with even a relatively limited number of partners was sufficient to yield increased service quality; having access to many more partners did not appear to produce further benefits ( Fig. 2k and l). In addition to benefiting from partner choice options, clients that were agile (i.e. had a high turning rate) and had the ability to exert control over the course of an interaction were better at preventing cleaners from cheating ( Fig. 2g and h). Indeed, cleaners could more easily consume mucus against a client's will if the client lacked the agility needed to avoid or terminate an interaction by outswimming a cleaner. Such control was not considered in classic BMT, which assumes that all interactions are voluntary (No€ e & Hammerstein, 1995). However, the ability of partners to avoid or terminate an interaction is an integral part of human market theory since variation in the degree of control over when to interact with whom (termed 'power') is an intrinsic characteristic of human interactions (Bowles & Hammerstein, 2003). The possibility that nonhuman animals can also exert variable amounts of control over cooperative interactions should be considered as an important mechanism in future studies of BMT if we aim to develop more realistic game-theoretical models (McNamara, 2013). For instance, the results of a theoretical model (Johnstone & Bshary, 2008) and laboratory experiments (Gingins et al., 2013) already suggested the importance of variation in clients' power in determining the quality of cleaning services provided by L. dimidiatus to its clients.
Our data indicate that a cleaner's temptation to cheat might be the most important factor affecting cooperation levels since it consists of three variables deemed equally important in determining the quality of cleaning services (Table 2). First, client body size has repeatedly been correlated with indicators of cleaning service quality, yielding positive, negative or no relationships (Adam, 2010;Bansemer, Grutter, & Poulin, 2002;Bshary, 2001;Grutter & Poulin, 1998a;Oates et al., 2010;Poulin & Rohde, 1997;Soares, Bshary, & Cote, 2008). We found a nonlinear relationship between client body size and the occurrence of jolts, suggesting that this relationship is complex and requires further consideration. Second, we found a positive relationship between the energetic value of a client species' mucus and the number of jolts experienced by individuals of that species. This result was expected based on theoretical predictions (Johnstone & Bshary, 2002), a previous laboratory experiment (Gingins et al., 2013) and field observations of a Mediterranean cleaner wrasse Symphodus melanocercus . Third, we found a negative relationship between the mean number of gnathiid ectoparasites a species harboured and the number of jolts it experienced. Gnathiids are the most common ectoparasite of fishes in the Indo-Pacific Ocean (Grutter, 1994;Grutter & Poulin, 1998b), and cleaners are reported to selectively feed on them (Grutter, 1997). This negative correlation fits previous studies documenting an increase in the duration of cleaning interactions (Grutter, 1995;Sikkel et al., 2004;Soares, Bshary, & Cote, 2008) and a decrease in cleaners consuming client mucus (Grutter, 1997) with an increasing gnathiid ectoparasite load on clients.

Parallels to Sexual Selection Theory
Our results support earlier claims that a biological market approach to the study of cooperation will benefit from drawing parallels with sexual selection theory (Barclay, 2013;No€ e, 2017;No€ e & Hammerstein, 1995). Until now, theoretical models of cooperation in a biological market context have focused almost exclusively on the effects of partner choice options (Ferriere et al., 2002;Hoeksema & Schwartz, 2001;Johnstone & Bshary, 2008;No€ e & Hammerstein, 1994). Focusing solely on partner choice is equivalent to considering the operational sex ratio (Emlen & Oring, 1977) as the only factor influencing mate choice in sexual selection theory. Importantly, however, students of sexual selection have long been aware that female mate choice is affected by other factors, such as a female's own quality (Cotton, Small, & Pomiankowski, 2006) and her ability to evade coercive mating by males (Brennan & Prum, 2012). A female's quality is a composite trait that might depend on a female's fecundity, her experience or how often she has recently mated (Clutton-Brock, 2009): high fecundity and high experience improve a female's quality, whereas frequent prior matings diminish it. Furthermore, females may experience reduced partner choice in a so-called 'restricted market' when males are (partly) able to perform sneaky copulations or coerce females into mating (Bowles & Hammerstein, 2003). The combination of these variables will affect male willingness to invest in reproductive effort (from parading to potential paternal investment). In the cleaning mutualism, visitor clients are the choosy traders (like females), and cleaners are the chosen traders (like males). A client's quality as a food patch is a complex trait: having a large body size, a high parasite load and mucus with a high caloric content makes a client attractive to cleaners. However, unlike female attractiveness, a client being attractive does not necessarily translate into greater cooperation by cleaners. This is because parasite load and mucus caloric content have opposing effects on service quality: having more parasites promotes cooperation whereas having mucus with a high energy content encourages cheating. The concept of a restricted market also adds explanatory power to the understanding of cleaner service quality: clients with few partner choice options are often coerced into interactions with a low service quality unless they are agile, which allows them to exert control over the course of an interaction (and therefore achieve a service quality beyond their market value).
In conclusion, our study illustrates the need for BMT to extend beyond partner choice as the sole mechanism regulating exchanges of goods or services to capture the full richness of cooperative interactions. A challenge for future BMT models is to integrate evolutionary and ecological timescales. In nature, cleaner to client ratios can change over a cleaner's lifetime (Triki, Wismer, Levorato, & Bshary, 2018), affecting decisions by clients and cleaners in this biological market . These changes in decisionmaking are best captured by ontogenetic learning models (Quiñones, Leimar, Lotem, & Bshary, 2020). Combining learning parameters with individual differences in cognition and resulting social competence (Taborsky & Oliveira, 2012;Triki, Emery, Teles, Oliveira, & Bshary, 2020) can set the stage for evolutionary learning models. Our current study raises the question of how the learning rules of cleaners can be fine-tuned through natural selection such that cleaners can adjust the quality of their cleaning services to specific client species as a function of client parasite load, mucus load, home range size and agility. Evolutionary learning models have been developed for other cooperation games such as the iterated prisoner's dilemma and the snowdrift game (Dridi & Akçay, 2018;Dridi & Lehmann, 2014), and may therefore serve as a starting point for the next generation of BMT models.

Author Contributions
D.G.R., M.J., A.S.G. and R.B. designed the study. D.G.R. and M.J. collected the data and performed the analyses. V.D. helped with the mucus protein and caloric content analysis. D.G.R., M.J. and R.B. wrote the paper with input from all co-authors.

Data Availability
The data and analysis script for this study are archived in the repository figshare (Roche, Jornod, Douet, Grutter, & Bshary, 2021) and were made available to the editors and referees upon initial submission.   High speed camera Stimulus in pipe Figure A2. The experimental arena used for testing the fast start escape performance of 13 client reef fish species. Two stimuli were suspended 70 cm above the water surface, 13 cm from the edge of the tank, on either side of the arena. A stimulus was released by remotely deactivating an electromagnet that fell inside an opaque PVC pipe. Two mirrors placed under each PVC pipe at a 45 angle allowed the recording of the exact time of contact with the water surface via a high-speed camera placed above the centre of the arena.   Figure A3. Correlation coefficients (r) between variables considered for inclusion in the full model to investigate which variables influence cooperation levels by cleaners: interaction duration (Duration) and client total length (TL), cumulative escape distance (D esc ), maximum speed (U max ), maximum acceleration (A max ), turning rate, mucus weight, mucus protein content (mucus protein), mucus caloric content (mucus calories), ectoparasites/cm 2 (ectoparasites) and gnathiid isopods/cm 2 (gnathiids). Turning rate Partner choice Figure A5. The percentage of cleaning interactions resulting in one or more jolts across 13 client reef fish species. Interactions ended with one of the following: one or both partners swimming away; the cleaner switching client; the client escaping the cleaner; the client chasing the cleaner. Species are ordered from lowest to highest percentage of interactions that ended with the client chasing the cleaner following a jolt. The number of observations per species is indicated at the bottom of each bar. Client manoeuvrability (turning rate) and partner choice options for each species are shown above the bars. Turning rate ( /s) was divided into five categories (very low, low, medium, high, very high) based on percentiles (0.2, 0.4, 0.6, 0.8). Partner choice options were categorized as low, medium and high (Table A2).