Ircinia_Community_manuscript.Rmd

---
title: "Sponge presence increases the diversity and abundance of fish and invertebrates in a subtropical seagrass bed"
output:   
  # bookdown::pdf_document2:
  bookdown::word_document2:
    number_sections: false
    reference_docx: "template.docx"
editor_options: 
  chunk_output_type: console
bibliography: Paper2.bib
csl: estuaries-and-coasts.csl
---
```{r setup, include=FALSE}
source("scripts_comm/02_community_data_org.R")
source("scripts_comm/08_import_results.R")

if(!require(kableExtra))install.packages("kableExtra");library(kableExtra)
```

Finella M. Campanino^1^, Philina A. English^2^, Craig A. Layman^3^, Stephanie K. Archer^1^*
  
1. Louisiana Universities Marine Consortium, Chauvin, Louisiana, USA
2. Pacific Biological Station, Fisheries and Oceans Canada, Nanaimo, British Columbia, Canada
3. Center for Energy, Environment, and Sustainability, Department of Biology, Wake Forest University, Winston-Salem, North Carolina, USA
  
*corresponding author: sarcher@lumcon.edu  

# Abstract  

Ecosystem engineers can have profound effects on biodiversity and community structure.
Seagrasses act as ecosystem engineers and support a wide variety of organisms.
Sponges possess many traits (e.g., create complex structure) that suggest they may be another example, but we know little about their relationships with animal communities in seagrass systems.
This study explored the effects of the marine sponge, *Ircinia felix*, on fish and invertebrate taxa richness, abundance, turnover, and community composition in a subtropical seagrass bed through a 1-year field experiment performed in The Bahamas.
We recorded the fish and invertebrate communities present in 5 × 5 m plots with the addition of either a live sponge, a polypropylene sponge replica (structure), or no additional structure (control).
Sponge presence alone explained the most variance in fish and invertebrate communities.
Fish and invertebrate taxa richness and abundance significantly increased over time in the presence of *I. felix*, but decreased or remained stable in the other plot types.
Sponge presence increased the number of new taxa gained, and reduced the number lost, for both fish and invertebrates.
Finally, we found that sponge presence shifted the entire fish community over time and across taxa space. 
Our study shows that sponges can act as ecosystem engineers in seagrass systems; however, additional research is needed to determine the full extent and implications of their effects.\par 

##### break

# Introduction 

Ecosystem engineers [*sensu* @Jones1994] are organisms that impact ecosystem structure and function by altering the physical and/or biological aspects of a system. 
They can create or modify entire habitats, ultimately influencing the distribution, diversity, and abundance of other organisms [@Byers2006]. 
@Jones1994 coined the term ecosystem engineer almost 30 years ago, however, we still lack an understanding of which taxa fill this role in many systems [@Byers2006]. 
More so, we lack an understanding of the multitude of ecosystem engineers within one system and their potential to affect the community and ecosystem characteristics interactively.
\par 

Seagrasses are ecosystem engineers that form meadows (or beds) supporting a wide range of organisms [@Hori2009]. 
They provide shelter to numerous species at multiple life stages [@James1994; @Nakamura2012], are a food source for animals [@Vicente1980], and play a key role in nutrient cycling [@Hemminga1991].
*Thalassia testudinum*, turtle grass, positively affects crustacean and fish abundance, species diversity, and species evenness when compared to sand habitat [@Arrivillaga1999; @Duffy2006; @Heck1980; @Orth1984].
<!-- Are we capitalizing common names? Journal may have a policy on this? -->
<!-- So after digging into this plants and invertebrates are lower case and the official common name (as indicated on fishbase.org) should be capitalized. The journal doesn't have a policy on this that I can find, but these are the policies set by the relevant societies (e.g., American Fisheries Society) -->
In addition, fish and invertebrate diversity and abundance are influenced by the morphological components and structural complexity of seagrass beds [@Duffy2006; @Hori2009].
Sponges are among the many organisms seagrass beds support and possess many traits (e.g., create complex structure) that suggest they may also act as ecosystem engineers in seagrass systems [@Bell2008a].\par

Sponges can add to the structural complexity of a habitat, alter localized water flow, stabilize substrate, and act as a settlement substrate for sessile organisms [@Bell2008a; @Wulff1984].
Some sponges have extensive microbial communities and macrofauna in and on their tissue [@Wulff2006; @Taylor2007; @Bell2008a] and can provide camouflage and protection to other organisms [@McClay1983; @Martin1992]. 
Sponges and their associated microbial communities can also act as a source of bioavailable forms of phosphorus and nitrogen in seagrass beds, especially important under oligotrophic conditions [@Archer2017; @Southwell2008; @Diaz1997].
Likely as a result of the this nutrient cycling, at least one sponge, *Ircinia felix* facilitates primary producers (seagrasses and macroalgae) in sub-tropical seagrass beds [@Archer2021a].
Combined, these traits suggest sponges have considerable potential to enhance the biodiversity of seagrass communities.\par

Despite their potential to significantly affect the distribution, diversity, and abundance of other organisms within a system, the effects of sponges on communities are often underappreciated [@Bell2008a].
Seagrass beds, where sponges can be common, are influenced by anthropogenic factors that threaten their biodiversity and function world wide [@Tin2021; @Orth2006; @Waycott2009]. 
Therefore, it is important to understand what drives diversity in these systems so that we can facilitate management and restoration. 
To this end, we investigated the effect of sponge presence on biodiversity within a relatively undisturbed seagrass ecosystem.\par

We investigated the impact of the sponge, *Ircinia felix*, on seagrass bed-associated communities. 
Using a 1 yr field-based experiment, we examined how sponge presence influenced fish and invertebrate taxa richness, abundance, taxa turnover, and community composition.
Because seagrass bed structural complexity and productivity affect fish and invertebrate communities, and *I. felix* increases growth and abundance of seagrasses and macroalgae, we further explored direct and indirect affects of sponge presence on the community.
We tested if the fish and invertebrate communities were influenced by the presence of the sponge alone (i.e., a direct effect), the enhancement of the dominant primary producers (macroalgae and seagrasses; i.e, an indirect effect), or a combination of the two.
We predicted that seagrass plots with *I. felix* would have higher taxa richness and abundance than seagrass plots without *I. felix* and that this would be mostly driven by indirect effects on the colonization and retention of species resulting from increases in structural complexity and productivity of the primary producers caused by sponge presence. 
As a result we hypothesize that the presence of *I. felix* would result in a change in community structure over time.


# Methods  

## Study site and experimental design  

This study was conducted in a continuous subtidal seagrass bed (~1 m low tide depth) on the southern region of Great Abaco Island, The Bahamas (26.02610 N, 77.37408 W). 
Within this seagrass bed, we delineated fifteen 5 × 5 m plots by placing wooden stakes at the corners and center of each plot on June 9, 2013. 
All plots were separated >2 m. 
We sampled the fish and epibenthic macroinvertebrate communities, as well as all explanatory variables (see below), once before the beginning of the experiment and again at 1 and 12 months after the treatments were established. 
Sponge (n=5), structure (n=5), and control (n=5) treatments were randomly assigned to plots. 
At the center of each sponge treatment plot, we placed one living sponge (*I. felix*, average volume ± standard deviation, 2.5 ± 0.75 L).
The structure treatment received a polypropylene model of a sponge, also placed at the center of each plot.
Both live and model sponges were contained within a cage to keep the sponge stationary.
The control plots remained unaltered.
Within the first month of the experiment, we replaced 3 live sponges when they showed signs of dying.
Live sponges were monitored throughout the rest of the experiment but no additional replacements were necessary.\par

## Faunal community

We monitored the response of the fish and epibenthic macroinvertebrate communities during each summer sampling event.
On each sampling occasion, all fish were identified to the lowest possible taxonomic unit, counted, and recorded within each plot during a five minute period.
Small silvery pelagic fishes which form large schools (e.g. fishes in the families Atherinopsidae and Clupeidae) were excluded from analyses following @Peters2015.
The macroinvertebrate community was quantified within a 3 m × 1 m band extending from the center of the plot (Figure S1).
All macroinvertebrates were identified to lowest possible taxanomic level and recorded.
If identification was not possible *in situ*, a representative sample was photographed and subsequently collected.
Colonial tunicates were excluded from analysis, as counts are not an accurate reflection of their abundance. 
Fish and macroinvertebrates were considered separately in all analyses. \par


## Primary producer covariates

@Archer2021a recently showed that *I. felix* facilitated primary producers in these same plots.
Seagrass growth, and seagrass and macroalgal abundance from @Archer2021a were used to calculate seagrass productivity, and seagrass and macroalgal structure for our analysis.
Seagrass productivity was calculated by multiplying the average growth per shoot by the total shoot density for *T. testudinum* for each plot at each sampling point.
Seagrass structure was calculated as the mean number of seagrass shoots (i.e., shoot density) for all seagrass species (*Halodule wrightii, Syringodium filiforme, and T. testudinum*) in twelve 20 × 20  cm “sub-quadrats” per plot at each sampling event.
Macroalgal structure was calculated as the total macroalgal abundance across all taxa at each sampling event.
Drift algae (e.g., *Karenia brevis*) were not included. 
Seagrass productivity, seagrass structure, and macroalgal structure were centered and standardized by two standard deviations to allow effect size comparisons between continuous and factor covariates within models [@Schielzeth2010]. \par

## Statistical Analyses

All statistical analyses were conducted in R version 4.1.0 [@RCoreTeam2021]. 

*Total abundance and taxa richness*

We analyzed total abundance and taxa richness of two faunal communities: fish and invertebrates.
First, we used analysis of variance (ANOVA) to confirm that none of the faunal communities differed in total abundance or taxa richness between the treatment plot types at the start of the experiment (sampling month = 0).
Next, we tested for differences in overall abundance and taxa richness of each faunal community using univariate Generalized Linear Mixed Models (GLMMs), with sponge presence as the reference level when treatment was included in the model.
We chose multiple model configurations to represent the biological hypotheses that could best explain community changes within seagrass plots.
Specifically, sponge presence, seagrass productivity, seagrass structure, and macroalgal structure were each modeled separately and additively to represent sponge presence alone, three characteristics of primary producers, and the combined effects of changes in primary producers and the presence of the sponge (Table S1-3).
All models included a random effect for plot, sampling month as a factor, and, for models including treatment, an interaction between treatment and time.
All models were fitted with a Conway–Maxwell–Poisson distribution, parametrized via the mean and a log link and implemented using the `glmmTMB` package [@Brooks2017a; @Huang2017].
We reviewed the residuals of all model configurations to confirm that GLMM assumptions were met.
Specifically, residuals were tested for compliance with model assumptions using the DHARMa package and a custom-built function to test for over-dispersion.
We used the package `AICcmodavg` [@Mazerolle2020] to create model selection tables based on corrected Akaike Information Criterion scores (AICc).
All models within two $\delta$ ***is this correct?*** AICc of the best-supported model are considered to be well-supported.\par

*Taxa Turnover*

Increases in taxa richness through time can result from either increases in colonization rates, or a reduction in the loss of taxa. 
To assess the relative rates at which taxa were gained or lost between sampling events, we used the `turnover` function in `codyn` package [@Hallett2020]. 
Using analysis of variance, we assessed both components of taxa turnover separately for each faunal community in each plot type, and compared patterns of turnover among treatments at 1 and 12 months into the experiment.
When significant differences were detected we compared treatments using a Tukey Honest Significant Difference test.\par

*Compositional Vectors*

To assess community changes in a way that simultaneously accounts for changes in taxa richness, taxa identity, and abundance, we calculated compositional vectors for each faunal community in each plot [@McCune2002].
We used Hellinger's transformation to help minimize the effect of large differences in abundance between taxa [@Legendre2001].
We then conducted a principal components analysis (unconstrained by environmental predictors) using the `rda` function in the `vegan` package [@Oksanen2020].
The scores from this analysis were used to calculate a vector length and angle from 0 to 12 months into the experiment for each plot as per @McCune2002. 
We performed an ANOVA on each faunal community vector length and angle to test for significant differences among treatments.
When there were significant differences, we used a post-hoc analysis via a Tukey Honest Significant Difference test.
For visualization purposes, all vectors were standardized to the same initial position and then plotted to show how much a community changed over time (i.e., vector length) and in which direction across taxa space the community shifted (i.e., vector angle).\par

We present effect sizes ($\beta$) and 95% confidence intervals (CI). 
All test statistics and p-values can be found in the Supplemental Material.\par

# Results  
```{r results function,include=FALSE}
decimalplaces <- function(x) {
  if ((x - round(x)) != 0) {
    strs <- strsplit(as.character(format(x, scientific = F)), "\\.")
    n <- nchar(strs[[1]][2])
  } else {
    n <- 0
  }
  return(n) 
}

# summary tables
aov.res<-function(rtable,row){
  tv<-round(rtable[row,4],2)
  pv<-ifelse(rtable[row,5]< 0.001, "p < 0.001", paste("p =", signif(rtable[row,5],1)))
  return(paste("=",tv,",",pv))
}
aov.df<-function(rtable,row){
  return(paste0(round(rtable[row,1],0),",",round(rtable[nrow(rtable),1],0)))
}

efsize<-function(rtable,row){
  betas<-signif(rtable$coefficients[row,1], 2)
  betas<-ifelse(abs(betas)<10,
    formatC(signif(betas,2), digits=2, format="fg", flag="#"), 
    round(betas))
  betas2 <- as.numeric(betas)
  decicount <- ifelse(betas2>1 & betas2 <10, 1, decimalplaces(betas2))
  se<-rtable$coefficients[row,2]
  cil<-rtable$coefficients[row,1]-1.96*se
  cil<-format(round(cil, digits=decicount), scientific=F)
  cih<-rtable$coefficients[row,1]+1.96*se
  cih<-format(round(cih, digits=decicount), scientific=F)
  return(paste0(betas,", ",cil," to ",cih))
}

efsize.tmb<-function(rtable,row){
  betas<-signif(rtable$coefficients$cond[row,1], 2)
  betas<-ifelse(abs(betas)<10,
    formatC(signif(betas,2), digits=2, format="fg", flag="#"), 
    round(betas))
  betas2 <- as.numeric(betas)
  decicount <- ifelse(betas2>1 & betas2 <10, 1, decimalplaces(betas2))
  se<-rtable$coefficients$cond[row,2]
  cil<-rtable$coefficients$cond[row,1]-1.96*se
  cil<-format(round(cil, digits=decicount), scientific=F)
  cih<-rtable$coefficients$cond[row,1]+1.96*se
  cih<-format(round(cih, digits=decicount), scientific=F)
  return(paste0(betas,", CI = ",cil," to ",cih))
}

# overall anovas
aov.res2<-function(rtable,row){
  tv<-round(rtable[row,5],2)
  pv<-ifelse(rtable[row,6]< 0.001, "p < 0.001", paste("p =", round(rtable[row,6],2)))
  return(paste("=",tv,",",pv))
}
aov.df2<-function(rtable,row){
  return(paste0(round(rtable[row,1],0),",",round(rtable[row+1,1],0)))
}

# tmb models
# overall anovas
tmb.res<-function(rtable,row){
  tv<-round(rtable$coefficients$cond[row,3],2)
  pv<-ifelse(rtable$coefficients$cond[row,4]< 0.001, "p < 0.001", paste("p =", round(rtable$coefficients$cond[row,4],2)))
  return(paste("=",tv,",",pv))
}

#Tukey results
tukey.res<-function(rtable,row){
  return(paste(signif(rtable[row,1],2),",",
               signif(rtable[row,2],2),"to", 
               signif(rtable[row,3],2),", p =",
               ifelse(rtable[row,4]<0.001,0.001,signif(rtable[row,4],2))))
}
```
  

## Total abundance and taxa richness 

*Fish*  

Throughout the experiment, a total of 15 fish taxa were observed.
At the start of the experiment between 0 and 2 individual fish, never of the same taxa, were found in each plot (mean = `r round(mean(fish.uni[fish.uni$sampling == 0, ]$f.abund,na.rm = FALSE),1)`, sd= `r round(sd(fish.uni[fish.uni$sampling == 0, ]$f.abund,na.rm = FALSE),1)`). 
*Halichoeres bivittatus* (Slippery dick) was the most common taxa. 
One month later, the maximum number of individuals per plot was 4, belonging to as many as 4 different taxa (mean abundance = `r round(mean(fish.uni[fish.uni$sampling == 1, ]$f.abund,na.rm = FALSE),1)`, sd = `r round(sd(fish.uni[fish.uni$sampling == 1,]$f.abund,na.rm = FALSE),1)` and mean richness = `r round(mean(fish.uni[fish.uni$sampling == 1, ]$spr,na.rm = FALSE),1)`, sd = `r round(sd(fish.uni[fish.uni$sampling == 1,]$spr,na.rm = FALSE),1)`) with *H. bivittatus* remaining the most abundant taxa across treatment types.
By the second summer of the experiment, the maximum number of individuals had increased to 7, belonging to as many as 6 different taxa (mean abundance = `r round(mean(fish.uni[fish.uni$sampling == 12, ]$f.abund,na.rm = FALSE),1)`, sd = `r round(sd(fish.uni[fish.uni$sampling == 12,]$f.abund,na.rm = FALSE),1)` and mean richness = `r round(mean(fish.uni[fish.uni$sampling == 12, ]$spr,na.rm = FALSE),1)`, sd = `r round(sd(fish.uni[fish.uni$sampling == 12,]$spr,na.rm = FALSE),1)`).
Juvenile grunts (*Haemulidae* spp.) were slightly more abundant than *H. bivittatus*.\par

Changes in fish abundance and fish taxa richness were best explained by sponge presence (Table 1, S1-3), such that in the presence of a sponge, both increased during the experiment (abundance: $\beta$ = `r efsize.tmb(fa.treat.sum,4)` and $\beta$ = `r efsize.tmb(fa.treat.sum,5)` at 1 and 12 months respectively, Fig. \@ref(fig:ab-plots); richness: $\beta$ = `r efsize.tmb(fspr.treat.sum,4)` and $\beta$ = `r efsize.tmb(fspr.treat.sum,5)` at 1 and 12 months respectively, Fig. \@ref(fig:turnover-plots)A).
By contrast, changes were non-significant for both control and structure plots (Figs. \@ref(fig:ab-plots) and \@ref(fig:turnover-plots)A).
Only control plots differed significantly from sponge plots after 12 months into the experiment (abundance: $\beta$ = `r efsize.tmb(fa.treat.sum,8)`; richness: $\beta$ = `r efsize.tmb(fspr.treat.sum,8)`).\par

*Invertebrates*

A total of 32 invertebrate taxa were observed in the experiment.
At the start of the experiment, between 0 and 4 individual invertebrates from up to 3 taxa were found in each plot (mean abundance = `r round(mean(inv.uni[inv.uni$sampling == 0, ]$i.abund),1)`, sd = `r round(sd(inv.uni[inv.uni$sampling == 0,]$i.abund),1)` and mean richness = `r round(mean(inv.uni[inv.uni$sampling == 0, ]$spr),1)`,sd = `r round(sd(inv.uni[inv.uni$sampling == 0,]$spr),1)`).
*Phrontis alba*, variable dog whelks, were the most abundant taxa. 
One month later, the maximum number of individuals per plot was up to `r range(inv.uni[inv.uni$sampling == 1, ]$i.abund)[2]` individuals belonging to up to `r range(inv.uni[inv.uni$sampling == 1, ]$spr)[2]` different taxa (mean abundance = `r round(mean(inv.uni[inv.uni$sampling == 1, ]$i.abund),1)`, sd = `r round(sd(inv.uni[inv.uni$sampling == 1,]$i.abund),1)` and mean richness = `r round(mean(inv.uni[inv.uni$sampling == 1, ]$spr),1)`, sd = `r round(sd(inv.uni[inv.uni$sampling == 1,]$spr),1)`) with more *Viatrix globulifera* (turtle grass anemones) and *Cerithium spp.* (ceriths) recorded than *P. alba*. 
By the second summer of the experiment, the maximum number of individuals was `r range(inv.uni[inv.uni$sampling == 12, ]$i.abund)[2]`, belonging to as many as `r range(inv.uni[inv.uni$sampling == 12, ]$spr)[2]` different taxa (mean abundance = `r round(mean(inv.uni[inv.uni$sampling == 12, ]$i.abund),1)`, sd = `r round(sd(inv.uni[inv.uni$sampling == 12,]$i.abund),1)` and mean richness = `r round(mean(inv.uni[inv.uni$sampling == 12, ]$spr),1)`, sd = `r round(sd(inv.uni[inv.uni$sampling == 12,]$spr),1)`) and *P. alba* were once again the most abundant taxa per plot, followed by *Callinectes sapidus* (blue crabs). \par

Changes in invertebrate abundance and taxa richness were also best explained by sponge presence (Table 1), but other models were equally well-supported for both abundance and taxa richness  (Table 1; S1).
Abundance of invertebrates significantly increased after 1 and 12 months in the presence of a sponge ($\beta$ = `r efsize.tmb(ia.treat.sum,4)` and $\beta$ = `r efsize.tmb(ia.treat.sum,5)` respectively, Fig. \@ref(fig:ab-plots)).
Only after 12 months were abundances significantly lower in both control and structure plots than in sponge plots ($\beta$ = `r efsize.tmb(ia.treat.sum,8)` and $\beta$ = `r efsize.tmb(ia.treat.sum,9)` respectively).
Seagrass productivity appeared to explain some of this change ($\beta$ = `r efsize.tmb(ia.treat.prod.sum,6)`), but there was still a significant difference between structure plots and sponge plots after accounting for the differences in seagrass productivity ($\beta$ = `r efsize.tmb(ia.treat.prod.sum,10)`).
Likewise, the model including seagrass structure was well-supported (Table 1); however, it's contribution did not reach significance ($\beta$ = `r efsize.tmb(ia.treat.struct.sum,6)`) suggesting that new growth rather than shoot density alone had some impact on invertebrate abundance. 
Invertebrate taxa richness also increased over time in the presence of a sponge ($\beta$ = `r efsize.tmb(ispr.treat.sum,4)` at 1 month and $\beta$ = `r efsize.tmb(ispr.treat.sum,5)` at 12 months, Fig. \@ref(fig:turnover-plots)).
After 12 months, taxa richness was lower in both control and structure plots than in the sponge plots, although only significantly so for the latter (control: $\beta$ = `r efsize.tmb(ispr.treat.sum,8)`; structure: $\beta$ = `r efsize.tmb(ispr.treat.sum,9)`).
However, the other well-supported model did not include a treatment effect, suggesting that the effect of the sponge might be accounted for through differences in macroalgal structure facilitated by the presence of a sponge ($\beta$ = `r efsize.tmb(ispr.alg.sum,4)`; Fig. \@ref(fig:inverts-by-algae)).

## Taxa Turnover  

Between the start of the experiment and the first sampling event, most plots across treatment types gained more taxa than they lost (Fig. \@ref(fig:turnover-plots) B left panel), and this resulted in no difference between treatments for both fish (gained: F~`r aov.df(fishgain1[[1]], 1)`~ `r aov.res(fishgain1[[1]], 1)`; lost: F~`r aov.df(fishloss1[[1]], 1)`~ `r aov.res(fishloss1[[1]], 1)`) and invertebrates (gained: F~`r aov.df(invgain1[[1]], 1)`~ `r aov.res(invgain1[[1]], 1)`; lost: F~`r aov.df(invloss1[[1]], 1)`~ `r aov.res(invloss1[[1]], 1)`).
After 12 months, however, sponge treatment plots tended to gain the most new fish (F~`r aov.df(fishgain12[[1]], 1)`~ `r aov.res(fishgain12[[1]], 1)` with sponge-control = `r tukey.res(fishgain12.tuk[[1]], 2)` and sponge-structure = `r tukey.res(fishgain12.tuk[[1]], 3)` while structure-control = `r tukey.res(fishgain12.tuk[[1]], 1)`) and invertebrate (F~`r aov.df(invgain12[[1]], 1)`~ `r aov.res(invgain12[[1]], 1)` with sponge-control = `r tukey.res(invgain12.tuk[[1]], 2)` and sponge-structure = `r tukey.res(invgain12.tuk[[1]], 3)` while structure-control = `r tukey.res(invgain12.tuk[[1]], 1)`) taxa.
Sponge plots also retained more fish taxa (F~`r aov.df(fishloss12[[1]], 1)`~ `r aov.res(fishloss12[[1]], 1)` with sponge-control = `r tukey.res(fishloss12.tuk[[1]], 2)` and sponge-structure = `r tukey.res(fishloss12.tuk[[1]], 3)` while structure-control = `r tukey.res(fishloss12.tuk[[1]], 1)`); Fig. \@ref(fig:turnover-plots) B right panel).
This pattern is similar for both fish and invertebrates, although overall the proportion of invertebrates lost did not differ between treatments (F~`r aov.df(invloss12[[1]], 1)`~ `r aov.res(invloss12[[1]], 1)`). \par


## Compositional Vectors

By 12 months into the experiment, the fish communities differed significantly between treatments for both vector length and angle, meaning that treatments differed in the how much they changed over time and in the direction of change across taxa space (vector length:  F~`r aov.df2(fish.vl.aov.sum,1)`~ `r aov.res(fish.vl.aov.sum,1)`; vector angle: F~`r aov.df2(fish.angle.aov.sum,1)`~ `r aov.res(fish.angle.aov.sum,1)`). 
The post-hoc analysis of vector lengths showed a significant difference in how much the fish communities changed over time between the control and sponge presence treatment plots and between the structure and sponge presence treatment plots, with control and structure treatment plots both having less change in community when compared to the sponge presence treatment plots (control-sponge: `r tukey.res(fish.vl.aov.tuk,1)`; structure-sponge `r tukey.res(fish.vl.aov.tuk,2)`). 
However, there was no difference in the amount of change in the fish community between structure and control treatment plots (`r tukey.res(fish.vl.aov.tuk,3)`). 
The post-hoc analysis of the vector angles showed no significant difference in the direction in which the fish communities shifted across taxa space between the control and sponge presence treatment plots (control-sponge: `r tukey.res(fish.angle.aov.tuk,1)`; Fig. \@ref(fig:vector-plots)), although this may be due to the small change in the control plot communities over time.
By contrast, there were differences in the direction of change in taxa space between the structure and sponge presence treatment plots, as well as between the structure and control treatment plots, with structure treatment plots having less change in communities than both sponge presence or control (structure-sponge: `r tukey.res(fish.vl.aov.tuk,2)`; structure-control `r tukey.res(fish.vl.aov.tuk,3)`). \par

By 12 months into the experiment, the invertebrate communities showed no significant differences between treatments for vector length and angle (vector length:  F~`r aov.df2(inv.vl.aov.sum,1)`~ `r aov.res(inv.vl.aov.sum,1)`; vector angle: F~`r aov.df2(inv.angle.aov.sum,1)`~ `r aov.res(inv.angle.aov.sum,1)`).
However, while control and structure control plots shifted in variable directions and amounts, all sponge plots tended to shift a similar amount in a somewhat more consistent direction (Fig. \@ref(fig:vector-plots)).
\par


# Discussion  

Ecosystem engineers can affect habitat structure, resource availability, and community structure and dynamics [@Jones1994; @Anderson1992; @Byers2006].
Here we present, to the best of our knowledge, the first *in situ* experiment examining if sponges act as ecosystem engineers in seagrass beds.
We did this by testing the affect of *I. felix* presence in subtropical seagrass beds on taxa richness, abundance, turnover, and community composition during a yearlong experiment.
We found that sponge presence alone explained the most variance in fish and invertebrate communities.
Overtime, sponge presence was associated with the addition of new fish and invertebrate taxa and the retention of previously present fish taxa.
This resulted in a shift in overall taxa space for the fish community, but not for the invertebrate community.\par

Although the sponge presence model explained the most variation in invertebrate taxa richness and abundance, other models fell within 2 $\delta \Delta$ ***both symbols are here*** AICc values and thus cannot be ruled out: macroalgal structure, sponge presence + seagrass productivity, and sponge presence + seagrass structure (Table 1). 
This suggests that sponge presence may impact the invertebrate community indirectly via the enhancement of primary producers. This is consistent with the presence of *I. felix* increasing the abundance macroalgae and both the abundance and productivity of seagrasses [@Archer2021a].  
Many studies have reported increasing invertebrate taxa richness and diversity related to the productive and structural complex habitat dense seagrass provides [@Heck1977;@Orth1984;@Edgar1992;@Webster1998;@Barry2021].
Our results are consistent with this well-supported paradigm.
However, the inclusion of sponge presence in two of three best supported models suggests that sponges may contribute directly to supporting invertebrate biodiversity, in addition to indirectly doing so through increasing the structural complexity of the seagrass bed.
Sponges host diverse invertebrate communities such as shrimps, amphipods, copepods, and brittle stars within their aquiferous systems and tissues [@Pearse1950;@Duffy1996;@Ribeiro2003].
In fact, *I. felix* can host over 6 species of infaunal invertebrates at densities up to 33 individuals per cm^-3^ of sponge [@Greene2008].
Although we did not account for these infaunal invertebrate communities in this study, these infaunal residents play a role in the link between sponges and the abundance of larger fauna.
For example, the snapping shrimp living within large sponges contribute to the auditory cues that juvenile spiny lobsters (*Panulirus argus*) use to identify suitable settlement habitat [@Butler1995; @Butler2016a].
<!-- PE: Should we stay consistent with putting the common names in () here? -->
<!-- I left this as is because we refer to them by the common name in the next sentence. I don't think it causes confusion for the reader, but if others feel strongly I'm happy to change it. -->
While spiny lobsters were not a large contributor to the patterns observed in this study (only one was observed--in a sponge plot), the use of sound as a settlement cue is common within marine taxa [@Mann2007;@Lillis2014;@Radford2007;@Vermeij2010].
Sponge presence can result in increased soundscape complexity in nearshore environments [@Butler2016a], but it is unclear if this mechanism influenced our results. 
<!--Marine sponges can have a range of microbes and macrofauna (e.g., invertebrates) living on and in their tissues[@Wulff2006; @Taylor2007; @Bell2008a] which could explain the increase in invertebrate taxa richness and abundance via habitat availability.
SA: this sentence feels out of place - how is this related to how invertebrates respond to increased seargass structure? -->
<!-- Notably, we did not measure sediment infauna invertebrate communities nor invertebrates residing within the sponge, so the influence of the sponge on invertebrate communities could be much larger than what we present. -->
<!-- SA: this sentence feels out of place - how is this related to how invertebrates respond to increased seargass structure? Same comment as previous sentence. Make sure each paragraph has a topic sentence and a related conclusion and everything in that paragraph expands on or provides more detail about those points. This paragraph drifts from the impact of sponges on the invertebrate community may be indirect to sponge infauna to sediment infauna to we need more research.-->
Regardless of the mechanism responsible, the presence of a sponge contributed to increased richness and diversity of invertebrates in this study. \par

As for invertebrates, the presence of a sponge increased taxa richness and abundance of fishes.
@Ferrari2018 showed that fish abundance could be best explained by benthic biota and complexity metrics in subtropical reefs. 
@Gratwicke2005 also found similar results with fish species richness being positively correlated with rugosity in subtropical habitats, including seagrass beds.
However, our results suggest *I. felix* effects on fish taxa richness and abundance exceed what would be expected from structure alone as the presence of a sponge replica did not result in increased fish richness or abundance.
It is possible that the increase in fish richness and abundance is related to the increase in invertebrates.
@Brook1977 suggested that fish populations in *T. testudinum* dominant seagrass beds were limited by abundance of polychaetes and crustaceans. 
@Yeager2012 results parallel the hypothesis that higher food availability associated with increase seagrass cover positively correlates with abundance of white grunts (*Haemulon plumierii*).
Where sponges were present, the fish community included species typically associated with reefs, such as grunts, damselfish (*Pomacentridae* spp.), and Nassau grouper (*Epinephelus striatus*)[@Nagelkerken2008; @CocheretdelaMoriniere2003].
These sponge plot-associated taxa all rely, at least in part, on invertebrate prey.
Juvenile grunts use seagrass as nursery habitat [@Arrivillaga1999; @Yeager2012; @CocheretdelaMoriniere2003], and they and beaugregory damselfish (*Stegastes leucostictus*, observed in sponge plots) partially rely on invertebrates as a food source [@Cervigon1966]. 
Nassau groupers associate with structure and use seagrass beds as feeding grounds with juvenile and midsize Nassau grouper feeding crustaceans [@Eggleston1998].
Likely, the increased taxa richness and abundance of invertebrates with sponge presence contributed to the colonization of structure-associated fish species that could take advantage of increased invertebrate abundance.\par

It is important to understand the effects of ecosystem engineers on the organisms around them because it can result in changes in the community structure and dynamics, ultimately influencing ecosystem functions and services provided by the community.
We studied the effect of the presence of an ecosystem engineer, *I felix*, on low-impact seagrass beds and found that sponge presence can increase fauna richness and abundance, reduce turnover rates, and shift fish communities.
This study is another *in situ* example of how an ecosystem engineer can affect an entire community, even when other, seemingly dominant, ecosystem engineers exist within that system.
Overall, we found that sponge presence was correlated with taxa retention and taxa gain.
For fish and invertebrates, we saw a shift towards taxa gained after 12 months (Fig. \@ref(fig:turnover-plots)), while control and structure treatment plots showed loss in taxa (Fig. \@ref(fig:turnover-plots)). 
This indicates that sponges may assuage conditions that lead to taxa loss in sub-tropical seagrass beds and suggests that sponges should be considered when managing seagrass ecosystems. 
<!-- SA: This sentence needs to be followed up with a sentence or two on what conditions might lead to taxa loss in sub-tropical seagrass beds - there's lots of research out there about drivers of biodiversity in seagrass beds -  and then with a couple sentences addressing how sponges might alleviate those conditions -->
Where sponges are abundant they enhance diversity of the community around them [@Maldonado2017; @Beazley2013]. 
Here we show that even a single sponge may sufficiently alter local conditions such that shifts in the community structure and ecosystem processes occur.
Thus, their role in providing structural complexity, increasing habitat availability, and enhancement of primary producers may help alleviate taxa loss in sub-tropical seagrass beds [@Ferrari2018;@Gratwicke2005;@Wulff2006; @Taylor2007; @Bell2008a; @Archer2021a]. \par


# Acknowledgements  
We would like to thank Friends of the Environment (NGO, Abaco, The Bahamas), Diane Claridge and Charlotte Dunn for their logistical support, Erik Archer, Elizabeth Whitman, and Ryann Rossi for their assistance in the field, and Katie Lewia and Jillian Tucker for their assistance in the lab. This work was supported by donations from Win and Tana Archer, North Carolina State University, and NSF OCE 1405198.\par
<!-- and the reviewers for their help improving this manuscript -->
<!-- This gets added after the first round of revision -->
  
# References  
<div id="refs"></div>  

##### break  
Table 1: Best supported models (those within two *∆ AICc* of the top model) of changes in abundance and taxa richness.
<!-- Need to change above to knit to pdf in order to get the table -->
<!-- Table 1: Best supported models (those within two *$delta$ AICc* of the top model) of changes in abundance and taxa richness. -->

```{r table, echo=FALSE, message=F, warning=FALSE, results='asis'}
# best.models<-best.models%>%
#   filter(Community!="Clonal Invertebrates")%>%
#   mutate(Community=ifelse(Community=="Fish","Fish","Invertebrates"))
# 
# kable(best.models[,c(-1:-2,-7,-10)],
#                           format = "latex",
#     booktabs=T,
#     escape = FALSE,
#     col.names = c("","$k$","$AIC_c$","$\\Delta AIC_c$",
#                   "$Akaike$ $Weight$","$LL$")#,
#     # caption="Best models, i.e., those within two
#     # $\\Delta AIC_c$ of the top model, explaining changes in fish
#     # and invertebrate abundance and taxa richness as well
#     # as clonal invertebrate taxa richness."
#     )%>%
#   # kableExtra::kable_styling(latex_options = "scale_down",font_size=11)%>%
#   kableExtra::pack_rows("Abundance",1,4,bold=T)%>%
#   kableExtra::pack_rows("Fish",1,1,bold=F,italic=TRUE)%>%
#   kableExtra::pack_rows("Invertebrates",2,4,bold=F,italic=TRUE)%>%
#   kableExtra::pack_rows("Taxa richness",5,7,bold=T)%>%
#   kableExtra::pack_rows("Fish",5,5,bold=F,italic=TRUE)%>%
#   kableExtra::pack_rows("Invertebrate",6,7,bold=F,italic=TRUE)#%>%
#   #save_kable(file = "figures/table.1.pdf")
# #   as_image()

knitr::include_graphics("figures/Table1-new.png")
```

##### break  

<!-- Figures -->

```{r ab-plots, fig.align='center', out.width="85%", echo=FALSE, message=FALSE, fig.cap="Change in abundance for fish (top panel) and invertebrates (bottom panel) at 0, 1, and 12 months into experiment."}

knitr::include_graphics("figures/fish_ncinvert_abundance figure.png")
```

##### break  

```{r turnover-plots, fig.align='center', out.width="85%", echo=FALSE, message=FALSE, fig.cap="Change in taxa richness (A) and the components of taxa turnover (B) for fish (upper panel of A and fish shapes in B) and invertebrates (lower panel of A and crab shapes in B) from start of experiment to the first experimental sample (1 month into experiment and left panel of B), and between that first experimental sample and one year later (12 months into experiment and right panel of B)."}

knitr::include_graphics("figures/Species_Richness_Turnover_A_B.png")
```

##### break  

```{r inverts-by-algae, fig.align='center', out.width="75%", echo=FALSE, message=FALSE, fig.cap="Invertebrate taxa richness can be explained by algae abundance almost as well as by experimental treatment. Line and shading represent the predicted relationship at 12 months after initiation of experiment from a GLMM with only algae abundance and time as fixed effects."}

knitr::include_graphics("figures/invert_spr_by_algae.png")
```

##### break  

```{r vector-plots, fig.align='center', out.width="55%", echo=FALSE, message=FALSE, fig.cap="Compositional vectors showing the direction and magnitude of change in the fish (top panel) and invertebrate (bottom panel) communities. The taxa contributing the most to the ordination results are shown in the biplots to the left of the compositional vector plots."}
# knitr::include_graphics("figures/community_vector_plots.png")
knitr::include_graphics("figures/community_vector_plots_withspp.png")
```