Influence of ecological traits on spatio-temporal dynamics of an elasmobranch community in a heavily exploited basin

Elasmobranchs, which include sharks and batoids, play critical roles in maintaining the integrity and stability of marine food webs. However, these cartilaginous fish are among the most threatened vertebrate lineages due to their widespread depletion. Consequently, understanding dynamics and predicting changes of elasmobranch communities are major research topics in conservation ecology. Here, we leverage long-term catch data from a standardized bottom trawl survey conducted from 1996 to 2019, to evaluate the spatio-temporal dynamics of the elasmobranch community in the heavily exploited Adriatic Sea, where these fish have historically been depleted. We use joint species distribution modeling to quantify the responses of the species to environmental variation while also including important traits such as species age at first maturity, reproductive mode, trophic level, and phylogenetic information. We present spatio-temporal changes in the species community and associated modification of the trait composition, highlighting strong spatial and depth-mediated patterning. We observed an overall increase in the abundance of the dominant elasmobranch species, except for spurdog, which has shown a continued decline. However, our results showed that the present community displays lower age at first maturity and a smaller fraction of viviparous species compared to the earlier observed community due to changes in species’ relative abundance. The selected traits contributed considerably to explaining community patterns, suggesting that the integration of trait-based approaches in elasmobranch community analyses can aid efforts to conserve this important lineage of fish.

. Variance partitioning among the explanatory variables included in the models. The panels show the results (a) for the presence-absence model and (b) for the conditional abundance model. The explanatory power is measured by Tjur R 2 for the presence-absence (PA) model and R 2 for the conditional abundance (ABU) model. The mean variance proportions averaged over the species are reported in brackets alongside the legend. 'Random: haul' indicates haul-level spatial random effects, while 'Random: year' represents temporal random effects. Figure S4. Predicted species mean abundance across the study area. Predictions refer to the year 2019. In the predictions, the swept area value is fixed to the mean of the trawl hauls used in the study (i.e., 0.047 km 2 ). Figure S5. Heatmap of the estimated parameters (regression slopes) which link species traits to species niches. The left panel represents the presence-absence (PA) model while the right panel the conditional abundance (ABU) model. Blue color indicates parameters that are estimated to be negative with at least 95% posterior probability, while responses that did not gain strong statistical support are shown in white. No parameter is estimated to be positive with at least 95% posterior probability.

Models' equations
We modeled the occurrence and abundance of each species (denoted as ) in each sampling unit (denoted as ) using a generalised linear model where is the statistical distribution (i.e., probit for the presence-absence model and normal for the conditional abundance), is the linear predictor, and 2 is the variance term (which is excluded for the probit model).
The linear predictor is modeled as the sum of fixed and random effects The fixed effects were modeled as a regression

= ∑ (3)
Here, represents the covariate measured at site , and represents the response of species to covariate , where 1 = 1 denotes the intercept.
The species' response to covariates is assumed to follow a multivariate normal distribution Here, .. denotes the vector of regression coefficients for all species' response to the covariates, which can be interpreted as species environmental niches. The symbol ⊗ denotes the Kronecker product, and 0 ≤ ≤ 1 measures the strength of the phylogenetic signal.
The expected niche models the influence of species-specific traits on species' responses, with where is the value of trait for species (with 1 = 1 denoting the intercept), and measures the effect of trait on the response to covariate . In this equation, denotes the phylogenetic covariance matrix, and denotes the identity matrix. When = 0, the residual variance is independent among the species, implying that closely related species do not have more similar environmental niches than distantly related ones. When approaches = 1, species' environmental niches are fully structured by their phylogeny, with related species having more similar niches than expected by chance, implying niche conservatism.
The random effects are modeled as where , is the linear predictor related to the random effect . In this case, = 1 represents the random effect of haul, and = 2 represents the random effect of year.
We define "units" as the sampling units at the haul and year level and denote the unit behind sampling unit as ( ). The random effect number is then defined as , = ∑ ( )ℎ ℎ=1 ℎ (7) Here, the summation goes over , which are the number of factors included for random effect . The ( )ℎ are the "unit loadings," and the ℎ are the "species loadings". The unit loadings ( )ℎ have an exponentially decaying correlation structure, where .ℎ ∼ (0, ) (8) and ∑ = 1 2 (− 1 2 / ℎ ) (9) Here, 1 2 is the distance in space or time between the units 1 and 2 , and ℎ is the spatial or temporal scale associated with the factor number ℎ of the random effect .