Algal Epibionts as Co-Engineers in Mussel Beds : Effects on Abiotic Conditions and Mobile Interstitial Invertebrates

Mussels and macroalgae have long been recognized as physical ecosystem engineers that modulate abiotic conditions and resources and affect the composition of rocky shore assemblages. Their spatial distributions in the intertidal zone frequently overlap, as many algal species thrive as epibionts on mussel beds. Nonetheless, their potential for combined engineering effects has not been addressed to date. Here we illustrate that Porphyra sp.—a desiccation-resistant macroalga that develops mostly epiphytically onto mussel beds—affects temperature, desiccation levels, and mobile interstitial invertebrates in mussel beds. Specifically, we observed that Porphyra cover (a) reduced temperature at the surface of the mussel bed but not at their base, (b) reduced desiccation both at the surface and base of the mussel bed and, (c) increased the densities of an abundant interstitial species—the amphipod Hyale grandicornis—in several study sites/dates. Additionally, we found that the positive responses of these grazing amphipods to Porphyra were driven by physical habitat modification (engineering) rather than food availability. This suggests that co-engineering by Porphyra and mussels generates abiotic states and focal species responses that would not be predictable from their individual effects. We expect that increased appreciation of co-engineering aids our understanding of complex ecological dynamics.


Introduction
Physical ecosystem engineers are organisms that structurally modify the environment via their presence or their activities [1,2].In so doing, they affect the resources and abiotic conditions on which other organisms depend [1][2][3][4].The organisms affected by engineers often include other engineers (e.g., trees create living space for nesting woodpeckers which, in turn, create holes in the tree trunk).In such cases, primary and secondary engineers can collectively affect other organisms.These compound engineering effects are frequently sequential, as in the case of habitat cascades (sensu [5]; hierarchical facilitation or facilitation cascades are often used as synonyms e.g., [6,7], but see [5]) where a primary engineer (e.g., a tree) creates living space or ameliorates physical conditions for secondary engineer (e.g., woodpecker or epiphytes) that in turn creates living space or ameliorates physical conditions for other organisms (e.g., arthropods) (Figure 1a; see [5][6][7][8][9][10][11]).Alternatively, primary and secondary engineers could also have concurrent, combined impacts on one or a few specific abiotic conditions and resources, with secondary engineers either ameliorating or worsening the primary engineering impacts on a focal species or assemblage (Figure 1b; e.g., [12][13][14]).
Diversity 2018, 10, x FOR PEER REVIEW 2 of 19 conditions for other organisms (e.g., arthropods) (Figure 1A; see [5][6][7][8][9][10][11]).Alternatively, primary and secondary engineers could also have concurrent, combined impacts on one or a few specific abiotic conditions and resources, with secondary engineers either ameliorating or worsening the primary engineering impacts on a focal species or assemblage (Figure 1B; e.g., [12][13][14]).Here we test if epibiotic algae acts as co-engineer in mussel beds by modifying their abiotic environment and affecting the mobile interstitial macroinvertebrates (i.e., invertebrates retained on 500 μm mesh; hereafter invertebrates).Mussels are dominant ecosystem engineers in the intertidal zone of temperate rocky shores across the five continents (see [15][16][17] for reviews).They form dense three-dimensional beds, thus generating interstitial space where water flow is attenuated, sediments become trapped, and the impacts of predators, extreme temperatures, and desiccation are all reduced (see reviews in [16][17][18]).By way of these habitat changes, mussel beds enhance the abundance of small mobile macroinvertebrates in rocky intertidal areas and also allow for the occurrence of interstitial species that would not otherwise be able to colonize these habitats, thus increasing their overall species richness [19][20][21][22].
At the same time, mussel shell surfaces are often suitable for the attachment of sessile invertebrates and algae (see [15][16][17]), to the extent that many sessile species are more successful in mussel-covered areas than bare rock surfaces [23,24].For instance, red algae typically ascribed to the genus Porphyra (hereafter Porphyra; but see [25]) commonly succeed as epibionts on mussel beds (e.g., [24,[26][27][28]).As mussels, algae are also typical rocky shore engineers and modify the physical environment in comparable ways (i.e., colonizable surface, interstitial space, flow attenuation, and refugia; e.g., [29][30][31][32]).In the particular case of Porphyra, many species in this group thrive at mid to high intertidal elevations because they resist heat and desiccation levels that are usually stressful to other algae [33].This suggests that Porphyra cover could afford additional protection from extreme heat and desiccation to the mobile interstitial fauna of mussel beds.
Porphyra thalli occur in the mid intertidal zone of several rocky intertidal shores of Buenos Aires Province, Argentina, from late winter to late summer [34,35].The majority of Porphyra thalli (> 75 %) Figure 1.Compound ecosystem engineering impacts on other organisms: (a) A habitat cascade where a primary engineer (EE1) alters abiotic conditions or resources (AR), thus facilitating a secondary engineer (EE2) that in turn creates abiotic conditions or resource levels (AR) favorable to other organisms (Other); (b) Co-engineering where a primary engineer (EE1) facilitates a secondary engineer (EE2) via changes in abiotic conditions and resources (AR), and they both have concurrent, joint impacts on abiotic conditions and resources (AR) that affect other organisms (Other).
Here we test if epibiotic algae acts as co-engineer in mussel beds by modifying their abiotic environment and affecting the mobile interstitial macroinvertebrates (i.e., invertebrates retained on 500 µm mesh; hereafter invertebrates).Mussels are dominant ecosystem engineers in the intertidal zone of temperate rocky shores across the five continents (see [15][16][17] for reviews).They form dense three-dimensional beds, thus generating interstitial space where water flow is attenuated, sediments become trapped, and the impacts of predators, extreme temperatures, and desiccation are all reduced (see reviews in [16][17][18]).By way of these habitat changes, mussel beds enhance the abundance of small mobile macroinvertebrates in rocky intertidal areas and also allow for the occurrence of interstitial species that would not otherwise be able to colonize these habitats, thus increasing their overall species richness [19][20][21][22].
At the same time, mussel shell surfaces are often suitable for the attachment of sessile invertebrates and algae (see [15][16][17]), to the extent that many sessile species are more successful in mussel-covered areas than bare rock surfaces [23,24].For instance, red algae typically ascribed to the genus Porphyra (hereafter Porphyra; but see [25]) commonly succeed as epibionts on mussel beds (e.g., [24,[26][27][28]).As mussels, algae are also typical rocky shore engineers and modify the physical environment in comparable ways (i.e., colonizable surface, interstitial space, flow attenuation, and refugia; e.g., [29][30][31][32]).In the particular case of Porphyra, many species in this group thrive at mid to high intertidal elevations because they resist heat and desiccation levels that are usually stressful to other algae [33].This suggests that Porphyra cover could afford additional protection from extreme heat and desiccation to the mobile interstitial fauna of mussel beds.
Porphyra thalli occur in the mid intertidal zone of several rocky intertidal shores of Buenos Aires Province, Argentina, from late winter to late summer [34,35].The majority of Porphyra thalli (> 75%) develops epiphytically on mussels, Brachidontes rodriguezii [35], which occupy most of the primary rock substrate at the mid intertidal elevations of these shores [34,36,37].Beds of B. rodriguezii host high densities of mobile interstitial invertebrates [34,36,38].Porphyra thalli can cover as much as 85% of the area in these beds, with peaks in coverage taking place between late October and early January [35].Since Porphyra coverage peaks concurrently with high summer temperatures, maximum sun irradiance, and the longest photoperiods (see SoDa-Solar Radiation Data; www.soda-pro.com),we predicted that this alga might reduce heat and desiccation impacts on the interstitial macrofauna of mussel beds, thus acting as a co-engineer in habitats where mussels are the primary, basal engineers.
To test our prediction, here we evaluated: a. the effects of Porphyra cover on desiccation and temperature above and underneath the mussel bed-we compared experimentally created Porphyra patches with exposed mussel bed areas late in the Porphyra growing season (February-March), when temperatures are high but algal cover declines (<25%); b.
the effects Porphyra cover on interstitial mussel bed invertebrates-we conducted short-term, Porphyra-removal experiments during December-January, when Porphyra cover peaks and desiccation and heat are likely maximum (see above); and c.
the occurrence of Porphyra-related patchiness in interstitial invertebrate assemblages-we sampled Porphyra-covered and exposed mussel bed areas that naturally occur late in the Porphyra growing season, when algal cover is low and patchy (see (a)).
Additionally, we deployed structural Porphyra mimics in the field to test whether the positive responses of a dominant grazer to Porphyra cover are driven by engineering (i.e., physical habitat provision) or trophic effects (i.e., food availability).

Study Sites and Organisms
This study was conducted at five rocky intertidal sites located along a ca.40 km coastal range between the cities of Mar del Plata (38 • 00 S, 57  1 for details).Tides along this range are semidiurnal and microtidal (0.80 m mean amplitude; Servicio de Hidrografía Naval, Argentina, www.hidro.gov.ar).The substrates at these sites are either orthoquartzite or Pampean loess cemented by calcium carbonate [39] (see Table 1).Brachidontes rodriguezii forms dense beds (up to 2000 in.dm −2 ) at the mid intertidal elevation of rocky shores along this coastal range [34,37,40].This relatively small mytilid (up to 55 mm length, most individuals < 30 mm length) occurs on the Atlantic coast of South America from the state of Rio Grande do Sul (Brazil, 32 • S) to Punta Ninfas (Argentina, 42 • S) [41,42].Their beds are primarily single-layered with multilayered areas restricted to protected vertical rock surfaces and small-sized (<50 cm 2 ) sparse hummocks (<1 per m −2 ) [37].
Porphyra thalli proliferates at our study sites from late winter (August) to early summer (January) and then declines to near zero cover by early fall (March).They usually are less than 15 cm long (pers.obs.) and mostly occur as mussel epibionts (>75% of cases), covering up to 85% of the mussel bed when laying flat during low tides [35].During its growing season, Porphyra is the only epibiont that establishes significant cover in mussel beds (other abundant mussel epibionts are either small-sized or occur at other times of the year; e.g., barnacles, Balanus glandula, and ulvoid algae, respectively).In previous studies, Porphyra specimens sampled within our study range and nearby locations were either reported as Porphyra umbilicalis [34], P. leucosticta [43], or P. pujalsiae [44].Moreover, many species in the genus Porphyra were recently assigned to the genus Pyropia based on molecular evidence [25].However, specimens from our area were not yet analyzed in this regard.Given such taxonomic uncertainties as well as frequent phenotypic variation within species in this group (see [45] and citations therein), we conservatively term the specimens found at our sites as Porphyra sp.(see also [46]).

Effects of Porphyra Cover on Desiccation and Temperature
The potential of Porphyra to reduce temperature and desiccation levels within the mussel bed habitat during air exposure was evaluated in separate, although analogous experiments conducted at Las Brusquitas on 22 February 16 and 16 March 2018, respectively.Daytime air exposure periods (low tides at 1:59 p.m. and 1:37 p.m., respectively) coincided with hot (maximum air temperatures at the site during the experimental period were 30 and 29 • C, respectively) and unclouded weather in both experimental dates.Each experiment consisted of (a) artificial Porphyra-covered patches (ca.75% algal cover), which were created by manually translocating thalli into 10 experimental plots (25 cm side squares) and, (b) an equal number of exposed mussel bed areas, which largely prevailed over natural Porphyra-covered patches in the experimental dates (note that both experiments were conducted late in the summer and Porphyra cover was low; see 2.1.).We refrained from using natural Porphyra-covered patches in this comparison because this would not allow distinguishing the effect of Porphyra cover per se from the effect of other microhabitat factors potentially leading to the patchy persistence of Porphyra by the end of its growing season.
The desiccation and temperature of thin standard materials (plastic cards wrapped in black electric tape, and water-saturated absorbent cleaning cloth, respectively; see below) were measured above and underneath the mussel layer both in Porphyra-covered plots and exposed mussel bed areas (see Figure 2).To place the standard materials under the mussel layer, a circular piece of the mussel bed was removed with the aid of a cylindrical core (10 cm diameter) and then replaced in its original position.The standard materials placed above the mussel layer in Porphyra-covered plots were deliberately covered by algal thalli, as the experiments aimed to evaluate Porphyra effects on the underlying mussel bed habitat (see Figure 2).Separate experimental plots were used for above and underneath measurements to ensure independence (see Figure 2).This led to a crossed, factorial design with two levels each of Porphyra cover (Yes/No) and Position (Above/Underneath mussels).
Plastic cards (5 cm side squares) wrapped with black electrical tape (3M ® ) were placed in each experimental plot (n = 6) by noon (ca.11:30 a.m.).The temperature of these cards was then measured at 12:30, 2:00, and 3:30 p.m. by means of an infrared thermometer (Fluke 59 MAX, Fluke Corporation, Everett, WA, USA).This thermometer allows adjusting temperature readings to the emissivity of the target material.Black electric tape was used here as a standard material (see also [47]) because of its known emissivity (0.95) and its suitability for placement under the mussel layer.At the time of measurement, the algal or mussel layer, or both, was removed as necessary, the cards were kept in its original position (i.e., lying on the rock or mussel substrate), and the infrared sensor of the thermometer was placed 15 cm above and perpendicular to the card so as to ensure a circular field of measurement of that do not exceed the card area (i.e., 2 cm diameter, 8:1 distance to spot ratio; see www.fluke.com).Algae and mussels were returned to their original place immediately after each measurement.
Pre-weighed (precision: 0.1 g), water-saturated pieces of absorbent cleaning cloth (5 cm side squares) were placed in each experimental plot (n = 6) before noon (11:00 a.m.) and retrieved and weighed again 3 h later.Desiccation was then measured as the percent water loss from cloth pieces during the experimental period.Although sponges have been traditionally used for desiccation measurements in rocky shore studies (e.g., [29]), absorbent cleaning cloth was preferred here because is thinner and can be placed underneath the mussel layer causing negligible elevation.
Diversity 2018, 10, x FOR PEER REVIEW 5 of 19 measurements in rocky shore studies (e.g., [29]), absorbent cleaning cloth was preferred here because is thinner and can be placed underneath the mussel layer causing negligible elevation.Variations in plastic card temperature were evaluated with a repeated-measures ANOVA model [48] with Porphyra cover (Yes/No) and Position of the plastic card (Above/Underneath the mussel layer) as fixed factors, and the Time of measurement (12:30, 3:00, and 4:30 PM) as repeated measure.Variations in percent water loss from absorbent cloth pieces were evaluated with two-way ANCOVA model [48] with Porphyra cover (Yes/No) and Position (Above/Underneath the mussel layer) as fixed factors and Initial Water Content as covariate.Differences between Treatments at each level of the other factor were evaluated with the Bonferroni multiple comparison procedure [48].All these analyses were run with Stata, Release 14 [49].

Effects of Porphyra Cover on Interstitial Mussel Bed Invertebrates
To evaluate the effect of Porphyra cover on the abundance of interstitial mussel bed invertebrates, short-term (ca.3-week long) Porphyra removal experiments were conducted at Punta Cantera, Las Brusquitas, and Copacabana by the time when cover of this species peaks (6 December 2011 to 27 December 2011, 22 December 2015 to 14 January 2016, and 30 November 2017 to 21 December 2017, respectively).This time of the year is characterized by high temperatures (maximum records during each period were 30, 37, and 34 °C, respectively) and other conditions that favor desiccation (maximum sun irradiance, longest photoperiods).Sixteen 25-cm-side, square plots were demarcated at each site.Porphyra thalli were pruned at their base from half of these plots by means of multipurpose scissors.The remaining plots were left undisturbed as controls (i.e., n = 8).Porphyraremoval plots were checked periodically (every 3-5 days) to remove any newly recruited or regrowing thalli.After ca. 3 weeks, a cylindrical core sample of the mussel bed (10 cm diameter) was taken from each plot (samples from control plots included both mussels and epiphytic Porphyra).These samples were taken to the laboratory and then sieved and rinsed over a mesh (500 μm) to collect the interstitial mobile macroinvertebrates.
Additionally, the Porphyra thalli removed from each Porphyra-removal plot at the beginning of the Punta Cantera experiment were taken to the laboratory and carefully rinsed over a mesh (500 μm) to collect and quantify the mobile invertebrates that were jointly removed with the algae.This allowed evaluating whether the comparison between treatment and control plots can be affected by incidental invertebrate removal when applying the algal removal treatment.The few, isolated thalli that recruited or re-grew in the Porphyra-removal plots during the course of the experiment were not Variations in plastic card temperature were evaluated with a repeated-measures ANOVA model [48] with Porphyra cover (Yes/No) and Position of the plastic card (Above/Underneath the mussel layer) as fixed factors, and the Time of measurement (12:30, 3:00, and 4:30 p.m.) as repeated measure.Variations in percent water loss from absorbent cloth pieces were evaluated with two-way ANCOVA model [48] with Porphyra cover (Yes/No) and Position (Above/Underneath the mussel layer) as fixed factors and Initial Water Content as covariate.Differences between Treatments at each level of the other factor were evaluated with the Bonferroni multiple comparison procedure [48].All these analyses were run with Stata, Release 14 [49].

Effects of Porphyra Cover on Interstitial Mussel Bed Invertebrates
To evaluate the effect of Porphyra cover on the abundance of interstitial mussel bed invertebrates, short-term (ca.3-week long) Porphyra removal experiments were conducted at Punta Cantera, Las Brusquitas, and Copacabana by the time when cover of this species peaks (6 December 2011 to 27 December 2011, 22 December 2015 to 14 January 2016, and 30 November 2017 to 21 December 2017, respectively).This time of the year is characterized by high temperatures (maximum records during each period were 30, 37, and 34 • C, respectively) and other conditions that favor desiccation (maximum sun irradiance, longest photoperiods).Sixteen 25-cm-side, square plots were demarcated at each site.Porphyra thalli were pruned at their base from half of these plots by means of multipurpose scissors.The remaining plots were left undisturbed as controls (i.e., n = 8).Porphyra-removal plots were checked periodically (every 3-5 days) to remove any newly recruited or re-growing thalli.After ca. 3 weeks, a cylindrical core sample of the mussel bed (10 cm diameter) was taken from each plot (samples from control plots included both mussels and epiphytic Porphyra).These samples were taken to the laboratory and then sieved and rinsed over a mesh (500 µm) to collect the interstitial mobile macroinvertebrates.
Additionally, the Porphyra thalli removed from each Porphyra-removal plot at the beginning of the Punta Cantera experiment were taken to the laboratory and carefully rinsed over a mesh (500 µm) to collect and quantify the mobile invertebrates that were jointly removed with the algae.This allowed evaluating whether the comparison between treatment and control plots can be affected by incidental invertebrate removal when applying the algal removal treatment.The few, isolated thalli that recruited or re-grew in the Porphyra-removal plots during the course of the experiment were not considered in this analysis as they were too small when removed (<3 cm length) and, then, unlikely to host a significant number of mobile invertebrates.
Differences in the assemblages associated to Porphyra-removal and control plots (Treatment) at each experimental Site/Date were evaluated with non-metric multidimensional scaling (nMDS) followed by two-way PERMANOVA (9999 permutations) and using the Bray-Curtis index as measure of dissimilarity [50,51].One-way PERMANOVA was used for individual, pair-wise comparisons of Porphyra-removal and control plots at each Site/Date (p-values adjusted by the sequential Bonferroni method; see [51]).Both analyses were run with PAST, Version 3.19 [52].
In addition, generalized linear models (GLM) [53] with a negative binomial distribution and a log link function were used to evaluate whether Porphyra removal (Treatment) and experimental Site/Date explain variations in the overall density of invertebrates and the density of common species (i.e., defined here as those that were present at all sites and showing mean overall densities in excess of 1 individual per sample).Model parameters were estimated with the maximum likelihood method using Stata, Release 14 [49].We used a negative binomial distribution as our data were counts showing either overdispersion or a relatively large proportion of zeroes (see [54]).The Deviance statistic was tested against the chi-square distribution to evaluate the effect of each explanatory variable and their interaction on the dependent variable [55].The explained deviance (D 2 , i.e., an equivalent to R 2 for GLM) was used as a measure of the amount of deviance accounted for by each variation of the model [56].Differences between Treatments at each Site/Date were evaluated with Bonferroni-adjusted contrasts of marginal linear predictors [55].

Porphyra-Related Patchiness in Interstitial Invertebrate Distribution
Invertebrate densities were compared between Porphyra-covered and exposed mussel bed areas.Samples (n = 6) were taken at Punta Cantera (05 February 2016), Los Acantilados (09 March 2018), Copacabana (08 March 2018), and Punta Hermengo (11 February 2016, Table 1) by means of a cylindrical core (10 cm diameter).The samples were then transported to the laboratory and sieved and rinsed over a mesh (500 µm) to obtain and quantify the macroinvertebrates.The samples of exposed and Porphyra-covered mussel bed areas respectively consisted of (a) mussels and their associated sediments and, (b) mussels, sediments, and Porphyra thalli.As mentioned above, we sampled our sites in the summer but late in the Porphyra growing season, as this is when Porphyra cover becomes patchy (<25%) and, thus, more likely to influence mobile invertebrate distributions.
Variations related to Patch Type and experimental Site/Date were evaluated here with the multivariate and univariate statistical methods used in 2.3.Again, a negative binomial distribution (log link function) was used in GLM, either to cope with overdispersion or a relatively large proportion of zero counts in the data [54].

Responses of a Dominant Grazer to Porphyra Thalli and Structural Porphyra Mimics
The amphipod, Hyale grandicornis, was positively associated to Porphyra both in the Porphyra removal experiment (2.3) and when comparing Porphyra-covered and exposed mussel bed areas (2.4) (see Results).To test whether the positive response of this grazing species to Porphyra is driven by engineering or trophic mechanisms, a short-term experiment (11 days) was conducted at Las Brusquitas from 12 December 2018 to 24 December 2018 (i.e., the time when Porphyra cover peaks and temperatures and desiccation are high; see 2.3).Twenty-one 25 cm side, square plots were demarcated and assigned to the following treatments (n = 7): (a) Porphyra removal and introduction of Porphyra mimics, (b) Porphyra removal, and (c) Porphyra-covered controls.Porphyra thalli in treatments (a) and (b) were pruned as in 2.3.Additionally, Porphyra mimics made of greenish brown, "rip-stop" nylon cloth were introduced to treatment (a).Foliose algae mimics made of nylon cloth were successfully used in previous studies of amphipod habitat selection (see [57]).Five to six mimics (cloth strips 8-12 cm long and 4-6 cm wide) were tied to a base made of plastic mesh (30 cm side squares, 2 cm mesh size), and the base with the mimics was nailed to the substrate from each corner.Mimics were mounted to a mesh base instead of individually driven in the plots to avoid disturbing mussel cover in the sampling area at the middle of each plot.The mimics and the mesh added to a coverage similar to that of naturally occurring Porphyra thalli at the time of the experiment (60-75% total cover, ca.10% mesh cover).The mesh base was not added to Porphyra-covered controls because this would lead to disturbance of naturally occurring Porphyra thalli.Mesh percent cover in the mimic treatment (ca.10%) was also minor relative to mimic (50-65%) and mussel cover (100%), and mesh threads were largely embedded within the mussel layer, which suggests negligible mesh effects via shading or structural addition.
As in 2.3, a cylindrical core sample of the mussel bed (10 cm diameter) was taken from each plot, and then sieved and rinsed over a mesh (500 µm) to collect the interstitial mobile macroinvertebrates.Variations in H. grandicornis densities across Treatments were again evaluated with a Generalized Linear Model (negative binomial distribution and log link function) followed by Bonferroni contrasts (see 2.3, 2.4).

Effects of Porphyra Cover on Desiccation and Temperature
Porphyra cover caused significant (5-7 • C) reductions in the temperature of plastic cards placed above the mussel bed but had no significant effect on the temperature of those placed underneath the mussel layer (Table 2, Figure 3).Plastic cards placed underneath the mussel layer showed lower temperatures than those placed atop, irrespective of Porphyra cover (Figure 3).These results were consistently observed across the three repeated measurements (Figure 3).
Porphyra cover significantly reduced water loss from absorbent cloth pieces both above and under the mussel layer (Table 3, Figure 4).The cloth pieces in Porphyra-covered plots showed larger water losses when placed underneath the mussel layer (Figure 4).On the other hand, water loss in exposed mussel bed areas was not significantly affected by the placement of cloth pieces relative to the mussel layer (i.e., above and underneath; Figure 4).Porphyra cover significantly reduced water loss from absorbent cloth pieces both above and under the mussel layer (Table 3, Figure 4).The cloth pieces in Porphyra-covered plots showed larger water losses when placed underneath the mussel layer (Figure 4).On the other hand, water loss in exposed mussel bed areas was not significantly affected by the placement of cloth pieces relative to the mussel layer (i.e., above and underneath; Figure 4).

Effects of Porphyra Cover on Interstitial Mussel Bed Invertebrates
The composition of interstitial invertebrate assemblages in Porphyra-removal and control plots (Treatments) at each experimental Site/Date showed varying degrees of separation after sample ordination (nMDS; Figure S1).Accordingly, there were significant interactive effects of Treatment and Site/Date on assemblage composition (Table 4).Pairwise analyses indicated that significant compositional differences among Treatments occurred at just one Site/Date (Las Brusquitas/2015; p = 0.003).

Effects of Porphyra Cover on Interstitial Mussel Bed Invertebrates
The composition of interstitial invertebrate assemblages in Porphyra-removal and control plots (Treatments) at each experimental Site/Date showed varying degrees of separation after sample ordination (nMDS; Figure S1).Accordingly, there were significant interactive effects of Treatment and Site/Date on assemblage composition (Table 4).Pairwise analyses indicated that significant compositional differences among Treatments occurred at just one Site/Date (Las Brusquitas/2015; p = 0.003).Total invertebrate densities varied mostly as a function of the Treatment X Site/Date interaction (64% of explained deviance; see Table 5).A posteriori pairwise contrasts indicated a negative effect of Porphyra removal on invertebrate densities at Las Brusquitas and lack of significant effects at the two other sites (Figure 5).

Porphyra-Related Patchiness in Interstitial Invertebrate Distribution
As in 3.2, the composition of interstitial invertebrate assemblages of Porphyra-covered and exposed mussel bed areas showed varying degrees of separation across experimental Sites/Dates after sample ordination (Figure S2).In concordance, there were significant interactive effects of Patch Type and experimental Site/Date on assemblage composition (Table 6).Pairwise analyses indicated significant compositional differences among Patch Types at Punta Cantera/2016 and Copacabana/2018 (p = 0.005 and 0.015, respectively), but not at the two other sites.7).Patch Type explained a minor amount of variation in invertebrate densities (6 % explained deviance) and the Patch Type X Site/Date interaction explained similar variation than Site/Date (see Table 7, Figure 6).The pulmonate limpet, Siphonaria lessonii, the amphipod, Hyale grandicornis, and the polychaete, Perinereis anderssoni, were the only species found at all experimental Sites/Dates and at overall densities in excess of one individual per sample (seven other species were found either at one or two Sites/Dates, low overall densities, or both; see Table S1).H. grandicornis densities varied mainly due to Porphyra removal (46% of deviance explained by Treatment, no difference in explained deviance with the Treatment X Site/Date interaction; see Table 5).Pairwise analyses indicated a negative effect of Porphyra removal on H. grandicornis densities at all Sites/Dates (Figure 5).On the other hand, variation in the densities of S. lessonii and P. anderssoni were best explained as a function of Site/Date (64 and 59% of explained deviance, respectively), with Treatment negligibly contributing to variation in the density of both species (<2% explained deviance in both cases, explained deviance similar for Treatment X Site/Date interaction and Site/Date as single explanatory term; see Table 5, Figure 5).
Porphyra removal at the onset of the Punta Cantera experiment just led to the incidental removal of 1.25 (SD = 1.28) and 3.37 (SD = 1.92) individuals of S. lessonii and H. grandicornis per plot (625 cm 2 ), respectively, or an equivalent to 0.16 (SD = 0.16) and 0.42 (SD = 0.24) individuals per core sample (78.5 cm 2 ).This is an order of magnitude lower than the final densities of both species in the samples from Porphyra-removal plots (see Figure 5) and, thus, a negligible influence on the experimental outcome.

Porphyra-Related Patchiness in Interstitial Invertebrate Distribution
As in 3.2, the composition of interstitial invertebrate assemblages of Porphyra-covered and exposed mussel bed areas showed varying degrees of separation across experimental Sites/Dates after sample ordination (Figure S2).In concordance, there were significant interactive effects of Patch Type and experimental Site/Date on assemblage composition (Table 6).Pairwise analyses indicated significant compositional differences among Patch Types at Punta Cantera/2016 and Copacabana/2018 (p = 0.005 and 0.015, respectively), but not at the two other sites.Total invertebrate densities varied mainly as a function of experimental Site/Date (55% of explained deviance; Table 7).Patch Type explained a minor amount of variation in invertebrate densities (6% explained deviance) and the Patch Type X Site/Date interaction explained similar variation than Site/Date (see Table 7, Figure 6).Table 7. Analyses of deviance for generalized linear models relating invertebrate counts (total invertebrates and three common species) with Patch Type (Porphyra-covered vs. Exposed) and experimental Sites/Dates (Punta Cantera/2016, Los Acantilados/2018, Copacabana/2018, and Punta Hermengo/2016).Models were fit using a negative binomial distribution and a log link function.Significant deviance values (Chi-squared test, p < 0.05, see asterisks) indicate lack of effect of explanatory variables or their interaction on the dependent variable.D 2 is the proportion between the deviance explained by the model in question and the deviance of the null, intercept model.

Responses of a Dominant Grazer to Porphyra Thalli and Structural Porphyra Mimics
Variations in the densities of the grazing amphipod, Hyale grandicornis, across Treatments significantly adjusted to a general linear model with a binomial error distribution and log link function (Deviance = 3.12, df = 18, p = 0.99, D 2 = 0.63).The density of this species did not differ between the Porphyra mimic treatment and the controls with natural Porphyra cover, whereas it was higher in these two treatments relative to Porphyra-removal plots (Figure 7).Invertebrate species other than H. grandicornis were found at very low densities (< 1 ind.per sample across treatments; see Table S3).Siphonaria lessonii, Hyale grandicornis, and Perinereis anderssoni were again the only species found at all experimental Sites/Dates and at overall densities in excess of one individual per sample (Eleven other species were found either at one or two Sites/Dates, low overall densities, or both; see Table S2).Variations in the density of H. grandicornis were best explained by the Patch Type X Site/Date interaction (56% of explained deviance; see Table 7).Densities of this species were higher in Porphyra-covered than exposed mussel bed areas both at Punta Cantera and Copacabana, whereas differences in density among Patch Types were not significant at the two other sites (Figure 6).Variations in the density of S. lessonii and P. anderssoni were best explained as a function of Site/Date (62 and 59% of explained deviance, respectively), with Treatment negligibly contributing to variation in the density of both species (1% explained deviance in both cases, explained deviance similar for Patch Type X Site/Date interaction and Site/Date as single explanatory term; see Table 7, Figure 6).

Responses of a Dominant Grazer to Porphyra Thalli and Structural Porphyra Mimics
Variations in the densities of the grazing amphipod, Hyale grandicornis, across Treatments significantly adjusted to a general linear model with a binomial error distribution and log link function (Deviance = 3.12, df = 18, p = 0.99, D 2 = 0.63).The density of this species did not differ between the Porphyra mimic treatment and the controls with natural Porphyra cover, whereas it was higher in these two treatments relative to Porphyra-removal plots (Figure 7).Invertebrate species other than H. grandicornis were found at very low densities (<1 ind.per sample across treatments; see Table S3).
Variations in the densities of the grazing amphipod, Hyale grandicornis, across Treatments significantly adjusted to a general linear model with a binomial error distribution and log link function (Deviance = 3.12, df = 18, p = 0.99, D 2 = 0.63).The density of this species did not differ between the Porphyra mimic treatment and the controls with natural Porphyra cover, whereas it was higher in these two treatments relative to Porphyra-removal plots (Figure 7).Invertebrate species other than H. grandicornis were found at very low densities (< 1 ind.per sample across treatments; see Table S3).

Discussion
The impacts of mussels on sessile epibiotic organisms and mobile interstitial invertebrates have traditionally been treated independently (see [15][16][17][18] and citations therein).However, they can be construed as interrelated in this example.Here, mussels and their algal epibionts act as primary and secondary engineers having combined effects and, thus, co-engineering the interstitial mussel bed habitat with consequences for some of their mobile invertebrate inhabitants (Figure 1b).The effects of epiphytic Porphyra cover on the mussel bed habitat-i.e., reduced temperature and desiccation-and their mobile invertebrate inhabitants-i.e., increased densities of an abundant amphipod-are discussed in detail below.

Effects of Porphyra Cover on Desiccation and Temperature
Temperature and desiccation were differently affected by Porphyra and mussel cover.Porphyra cover reduced temperatures at the mussel bed surface, which suggests that it effectively insulates mussels from direct solar radiation.Yet, temperatures were even lower under the mussel bed and irrespective of Porphyra cover, which suggests that mussel cover alone suffices to buffer heat levels near the rock surface.Gradients of decreasing temperature from the surface to the base of mussel beds were recently described for multilayered beds (12 cm depth) of a larger mussel species (Mytilus californianus, see [58]).From our findings, it becomes apparent that these vertical gradients also occur in thinner (<3 cm depth), single-layered beds of smaller mussel species (Brachidontes rodriguezii in this study) and that algal cover (here Porphyra sp.) can reduce temperatures at the upper end of the gradient and, thus, its overall amplitude.
On the other hand, Porphyra cover reduced desiccation both above and under the mussel bed (cf., temperature reductions only at the surface of the mussel bed).Interestingly, water losses from absorbent pads in Porphyra-covered areas were lower above than underneath the mussel bed.This can be counterintuitive if we consider that temperatures at the surface of the mussel bed are higher than at its base (see above).This could simply be attributed to water dripping from the alga onto the pads, though Porphyra thalli is known to rapidly dehydrate when air exposed (40% water loss after 30 min; see [59]) and they were externally dry when the experiment was set up.As a possible alternative explanation, we posit that Porphyra thalli laying flat at the surface of the mussel bed during low tide could be acting as "condensation traps", analogous to those made by humans for water collection and purification in arid regions (e.g., [60,61]).Specifically, we predict that moisture will evaporate as the air within the mussel bed heats up leading to vapor accumulation and condensation at the undersides of Porphyra thalli.In this way, Porphyra thalli could have helped maintain moisture in the absorbent pads placed immediately underneath them (i.e., beneath Porphyra thalli and the mussel bed, see Figure 2).A similar condensation-trapping mechanism has been proposed to contribute to self-irrigation in a prostrate, broad-leaved desert plant (the desert rhubarb, Rheum palaestinum; see [62]).
Here, condensation at the undersides of Porphyra thalli is expected to enhance the local retention of moisture, which may ultimately buffer the alga and underlying invertebrates from desiccation.

Mobile Invertebrate Responses to Porphyra Cover
The amphipod Hyale grandicornis was the only abundant mobile invertebrate species affected by Porphyra cover.The densities of this species decreased after Porphyra removal at all experimental sites/dates (Figure 5) and were higher in Porphyra patches than exposed mussel bed areas at two out of four sampling sites/dates (Figure 6).Differences in the invertebrate assemblages associated to Porphyra-covered and exposed mussel bed areas at some sites/dates (i.e, Porphyra removal experiment at Las Brusquitas and comparisons of natural Porphyra-covered patches and exposed mussel bed areas at Punta Cantera/2016 and Copacabana/2018; Tables 4 and 6, Figures S1 and S2) can largely be attributed to H. grandicornis, as they disappear when this species is excluded from the analysis (see Tables S4 and S5).Similar considerations apply to the decreased total invertebrate densities observed at Las Brusquitas/2015 after Porphyra removal.In this case, excluding H. grandicornis from the analysis leads to mean invertebrate densities of 2.50 (SD = 1.60) and 3.63 (SD = 2.33) individuals per sample in Porphyra-removal and control plots, respectively (cf.5.25 and 17.38 ind.per sample, respectively, when H. grandicornis is included; see Figure 5).
H. grandicornis densities were enhanced to a similar degree by natural Porphyra cover and by structural Porphyra mimics.This suggests that the positive response of H. grandicornis to Porphyra-covered areas is primarily driven favorable environmental conditions under algal thalli rather than food availability.This is not surprising if we consider that our samples were collected during daytime low tides and that hyalid amphipods mostly graze during high tides or at night (see [63,64]).Although H. grandicornis has been reported to graze on Porphyra spp. in South Africa [65], our findings align to those of previous studies suggesting that hyalid amphipods congregate under macroalgal canopies at low tide to avoid detrimental temperature and desiccation levels rather than feeding (e.g., [63,66]).
The positive responses of H. grandicornis to structural Porphyra mimics indicate that engineering mechanisms are implied in the facilitative effect of Porphyra on this species (see also [67,68] for examples of structural effects of macroalgae on intertidal amphipod densities).Considering the intertidal elevation where our experiment was carried out (i.e., the mid to high intertidal elevations where Porphyra occurs), we can expect this effect to be mainly a response to reduced heat and desiccation levels under Porphyra thalli (see also [63,66]).This positive engineering effect of Porphyra on amphipods should become less prevalent as the summer advances due to gradual reductions on Porphyra cover (see [35]).This suggests that extreme heat and desiccation would be more detrimental to amphipods in the late summer.Yet, the facilitative effects of Porphyra on H. grandicornis observed here to occur in late spring and early summer might still have important indirect effects on other ecosystem components (e.g., reductions in the biomass of algal food species, increases in the prey base available to littoral fishes and the overall transfer of energy to higher trophic levels (see [64,69]).
It is uncertain if engineering mechanisms other than the modulation of heat and desiccation have contributed to increased amphipod densities under Porphyra thalli.Known effects of epibiotic algal cover on the underlying mussel bed habitat include reduced impact of epibenthic predators, decreased water flow, and increased pH [70][71][72], all of which occur or prevail during tidal submergence, when hyalids are expected to leave the mussel matrix to forage (see [64]).However, protection from low tide predators (e.g., passerine birds such as the great kiskadee, Pitangus sulphuratus; pers.obs.) might possibly have contributed to this pattern.

Concluding Remarks: Co-engineering, Habitat Cascades, and More
A major research focus since the inception of the ecosystem engineering concept were the effects of individual species on their environment and other organisms (see [73] and citations therein).Yet, virtually all organisms engineer their environment to a higher or lesser degree [3], which implies that engineering species and their impacts rarely occur in isolation.In recent years, there has been growing interest in studying the compound impacts of multiple engineers on their environment and other organisms.A major emphasis in this regard has been the study of habitat cascades, that is, the positive effects on organisms that result from successive habitat creation or amelioration by hierarchy of primary to n th order engineers (sensu [5], e.g., [5][6][7][8][9][10][11]).In these cases, the responses of focal organisms are largely driven by the structural effects of the intermediate engineer.For example, phytotelmata communities in epiphytic tank bromeliads occur because of their ability to trap water (secondary engineering) rather than any structural attribute of the basal trees (i.e., the primary engineers; see Figure 1a).
Our example, however, departs from typical examples of habitat cascades as our secondary engineer-i.e., Porphyra-was relatively unimportant per se as physical habitat for interstitial mobile invertebrates (see 3.1).In this case, secondary engineers ameliorated physical conditions and facilitated other organisms within the primary engineered habitat-i.e., the mussel bed.Thus, abiotic conditions (temperature, desiccation) and interstitial invertebrate densities (specifically H. grandicornis) can be construed here as controlled by combined impacts of primary and secondary engineers (Figure 1b) that would not be predictable from their individual effects (cf., habitat cascades).We expect these non-additive, co-engineering impacts (termed also "cooperative engineering", see [13]) to be commonly driven by epibiotic engineers and their basibionts (see [12,74] this study) though not strictly so (e.g., combined effects of salt marsh plants and ribbed mussels on sediment properties and fiddler crab densities, [14]).We also expect that co-engineering cause both positive and negative effects on focal organisms (i.e., secondary engineers ameliorating or worsening the abiotic states created by the primary ones; cf., inherent focus of habitat cascades on positive impacts).
We certainly recognize that the above distinctions and generalizations are tentative and that they would benefit from further scrutiny, especially as regards to comprehensiveness and utility.Clearly, as the study of compound engineering impacts is still in its infancy, there is room for interactions and mechanisms that may not necessarily conform to the above categories (e.g., secondary engineers that modify or rearrange the structures made by primary ones, thus affecting other species; e.g., [75]).Our overall intent here is to stimulate further thought regarding co-engineering, habitat cascades, and other possible kinds of compound engineering impacts.In the long run, we hope that increased appreciation and conceptualization on this topic help scientists understand, integrate, and predict complex ecological dynamics.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1424-2818/11/2/17/s1, Figure S1.nMDS ordination of invertebrate assemblages in Porphyra-removal and control plots, Figure S2.nMDS ordination of invertebrate assemblages in Porphyra-covered patches and exposed mussel bed areas, Table S1: Invertebrate species found at low densities in the Porphyra-removal experiment, Table S2: Invertebrate species found at low densities in Porphyra-covered and exposed mussel bed areas, Table S3: Invertebrate species found at low densities in the experiment with structural Porphyra mimics, Table S4: Two-way PERMANOVA comparing the composition of interstitial invertebrate assemblages in the Porphyra-removal experiment when excluding the amphipod Hyale grandicornis from the analysis, Table S5: Two-way PERMANOVA comparing the composition of interstitial invertebrate assemblages in Porphyra-covered and exposed mussel bed areas when excluding the amphipod Hyale grandicornis from the analysis, Table S6.Names and taxonomic authorities of the species analyzed in this study.
Author Contributions: J.L.G. and M.G.P. conceived and designed the study.J.L.G. and M.B. carried out the field research and sample analysis.J.L.G. was primarily responsible for data analysis and the writing of manuscript drafts.M.B. and M.G.P. contributed to interpretation of findings and commented on manuscript drafts.All authors read and approved the final manuscript.

Figure 1 .
Figure 1.Compound ecosystem engineering impacts on other organisms: (a) A habitat cascade where a primary engineer (EE1) alters abiotic conditions or resources (AR), thus facilitating a secondary engineer (EE2) that in turn creates abiotic conditions or resource levels (AR) favorable to other organisms (Other); (b) Co-engineering where a primary engineer (EE1) facilitates a secondary engineer (EE2) via changes in abiotic conditions and resources (AR), and they both have concurrent, joint impacts on abiotic conditions and resources (AR) that affect other organisms (Other).

Figure 2 .
Figure 2. Treatments in the factorial experiments conducted to evaluate the effect of Porphyra cover (Yes/No) and place of measurement (Above/Under the mussel layer) on temperature and desiccation.The bold black line indicates the position of the standard material used to measure temperature and desiccation (Plastic cards wrapped with black electrical tape and water-saturated pieces of absorbent cleaning cloth, respectively).

Figure 2 .
Figure 2. Treatments in the factorial experiments conducted to evaluate the effect of Porphyra cover (Yes/No) and place of measurement (Above/Under the mussel layer) on temperature and desiccation.The bold black line indicates the position of the standard material used to measure temperature and desiccation (Plastic cards wrapped with black electrical tape and water-saturated pieces of absorbent cleaning cloth, respectively).

Figure 3 .Figure 3 .
Figure 3. Mean (SE) temperature of plastic cards placed above and underneath the mussel layer in Porphyra-covered and exposed control areas.Three measurements were made during a daytime, air

Figure 4 .
Figure 4. Mean (SE) water loss from water-saturate absorbent cloth pieces placed above and underneath the mussel layer in Porphyra-covered and exposed control areas.Cloth pieces were left in the field for three hours (11:00 AM to 2:00 PM) during a daytime, air exposure period (Low tide at 1:37 PM).Different letters above bars indicate significant differences.

Figure 4 .
Figure 4. Mean (SE) water loss from water-saturate absorbent cloth pieces placed above and underneath the mussel layer in Porphyra-covered and exposed control areas.Cloth pieces were left in the field for three hours (11:00 a.m. to 2:00 p.m.) during a daytime, air exposure period (Low tide at 1:37 p.m.).Different letters above bars indicate significant differences.

Figure 7 .
Figure 7. Mean (SE) densities of the grazing amphipod Hyale grandicornis in Porphyra-covered areas (control), Porphyra-removal plots, and removal plots where structural Porphyra mimics ("rip stop" nylon cloth strips) were added to a similar cover than controls (60-75%).The experiment was conducted at one site (Las Brusquitas) during December 2018, coinciding with peak Porphyra cover.Different letters above bars indicate significant differences.

Figure 6 .
Figure 6.Mean (SE) total invertebrate density and densities of common species (i.e., species present at all sites at mean overall densities in excess of 1 individual per sample) in Porphyra-covered patches and exposed mussel bed areas.Samplings were conducted during February-March when Porphyra cover becomes relatively low (<30% cover) and patchy, and at three sites/dates: Punta Cantera/2016 (PC), Los Acantilados/2018 (LA), Copacabana/2018 (CO), and Punta Hermengo/2016.Asterisks indicate significant different between treatments at a site/date.

Figure 7 .
Figure 7. Mean (SE) densities of the grazing amphipod Hyale grandicornis in Porphyra-covered areas (control), Porphyra-removal plots, and removal plots where structural Porphyra mimics ("rip stop" nylon cloth strips) were added to a similar cover than controls (60-75%).The experiment was conducted at one site (Las Brusquitas) during December 2018, coinciding with peak Porphyra cover.Different letters above bars indicate significant differences.

Figure 7 .
Figure 7. Mean (SE) densities of the grazing amphipod Hyale grandicornis in Porphyra-covered areas (control), Porphyra-removal plots, and removal plots where structural Porphyra mimics ("rip stop" nylon cloth strips) were added to a similar cover than controls (60-75%).The experiment was conducted at one site (Las Brusquitas) during December 2018, coinciding with peak Porphyra cover.Different letters above bars indicate significant differences.

Table 1 .
Study site coordinates, rock type, and field activities.

Table 2 .
Repeated-measures ANOVA model relating the temperature of a standard material (plastic cards wrapped in black electric tape) with Porphyra cover (Yes/No), Position of the plastic card (Above/Underneath the mussel layer), and the Time of measurement (12:30, 2:00, and 3:30 p.m.).Asterisks indicate significant effects.

Table 3 .
Two-way ANCOVA model relating water loss from pre-wetted absorbent cloth pieces (i.e., a measure of relative desiccation) with Porphyra cover (Yes/No), Position of the cloth pieces (Above/Underneath the mussel layer), and their Initial Water Content (continuous covariate, range 30.00-35.50 mL).Asterisks indicate significant effects.

Table 3 .
Two-way ANCOVA model relating water loss from pre-wetted absorbent cloth pieces (i.e., a measure of relative desiccation) with Porphyra cover (Yes/No), Position of the cloth pieces (Above/Underneath the mussel layer), and their Initial Water Content (continuous covariate, range 30.00-35.50 ml).Asterisks indicate significant effects.

Table 5 .
Analysis of deviance for generalized linear models relating invertebrate counts (total invertebrates and three common species) with Treatment (Porphyra removal vs. Controls) and experimental Sites/Dates (Punta Cantera/2011, Las Brusquitas/2015, and Copacabana/2017).Models.were fit using a negative binomial distribution and log link function.Significant deviance values (Chi-squared test, p < 0.05, see asterisks) indicate lack of effect of explanatory variables or their interaction on the dependent variable.D 2 is the proportion between the deviance explained by the model in question and the deviance of null, intercept model.
Figure 5. Mean (SE) total invertebrate density and densities of common species (i.e., species present at all sites at mean overall densities in excess of 1 individual per sample) in the Porphyra-removal treatment and control plots with natural Porphyra cover.Experiments were conducted during December-January when Porphyra cover peaks (60-80% cover) and at three sites/dates: Punta Cantera/2011 (PC), Las Brusquitas/2015 (LB), and Copacabana/2017 (CO).Asterisks indicate significant different between treatments at a site/date.

Table 7 .
Analyses of deviance for generalized linear models relating invertebrate counts (total invertebrates and three common species) with Patch Type (Porphyra-covered vs. Exposed) and experimental Sites/Dates (Punta Cantera/2016, Los Acantilados/2018, Copacabana/2018, and Punta Hermengo/2016).Models were fit using a negative binomial distribution and a log link function.Significant deviance values (Chi-squared test, p < 0.05, see asterisks) indicate lack of effect of * * * * Figure 5. Mean (SE) total invertebrate density and densities of common species (i.e., species present at all sites at mean overall densities in excess of 1 individual per sample) in the Porphyra-removal treatment and control plots with natural Porphyra cover.Experiments were conducted during December-January when Porphyra cover peaks (60-80% cover) and at three sites/dates: Punta Cantera/2011 (PC), Las Brusquitas/2015 (LB), and Copacabana/2017 (CO).Asterisks indicate significant different between treatments at a site/date.