Environmental and ecological factors mediate taxonomic composition and body size of polyplacophoran assemblages along the Peruvian Province

Intertidal communities’ composition and diversity usually exhibit strong changes in relation to environmental gradients at different biogeographical scales. This study represents the first comprehensive diversity and composition description of polyplacophoran assemblages along the Peruvian Province (SE Pacific, 12°S–39°S), as a model system for ecological latitudinal gradients. A total of 4,775 chitons from 21 species were collected on twelve localities along the Peruvian Province. This sampling allowed us to quantitatively estimate the relative abundance of the species in this assemblage, and to test whether chitons conform to elementary predictions of major biogeographic patterns such as a latitudinal diversity gradient. We found that the species composition supported the division of the province into three ecoregional faunal groups (i.e. Humboldtian, Central Chile, and Araucanian). Though chiton diversity did not follow a clear latitudinal gradient, changes in species composition were dominated by smaller scale variability in salinity and temperature. Body size significantly differed by ecoregions and species, indicating latitudinal size-structure assamblages. In some localities body size ratios differed from a random assemblage, evidencing competition at local scale. Changes in composition between ecoregions influence body size structure, and their overlapping produce vertical size segregation, suggesting that competition coupled with environmental conditions structure these assemblages.

The rarefaction curves obtained for each locality showed a saturation around 50 RA, reached between 6 to 10 species (Fig. S1); therefore, the estimated richness based on the RA matrix was strong enough to represent the diversity of chitons (Fig. S1).
Neither the ordinary least square regression analysis (OLS) nor the spatial autoregressive model (SAR) evidenced significant relationships between abundance, richness and eveness with the environmental variables evaluated (Table 2), due to a high spatial autocorrelation (ρ = 0.969). Only species composition (Jaccard) was significant related with temperature and salinity (Table 2).
Overall, most chitons were small (18 species; <50 mm TL), whith only three species (Acanthopleura echinata, Chiton magnificus and Enoplochiton niger; see Online Appendix Table S2, Fig. 7) corresponding to large chitons (>50 mm TL). Significant differences in chiton body lengths were found to be structured by species and ecoregions, with significant interactions between most factors ( Table 3). The year was the only factor without significant differences (Table 3). However, species was the only factor that explained a high proportion of deviance (56%; see Table 3). Mean minimum segment lengths in each locality were heterogeneous (59.1-88.8 mm), being higher at northern locations and decreasing through the south (Table S3). Maximum chitons body length at each locality was correlated with mean minimum segment length (r = 0.59, P < 0.05), suggesting that in places with larger chitons they are highly segregated by size. Body size ratios were highly variable (0.0-4.6) and showed significant differences from a random assemblage that was not structured by competition only in some localities, where we suggest some evidence of competition by size overlap at local scale (Table S3).

Discussion
Our results showed an absence of a latitudinal gradient in abundance and diversity of coastal polyplacophorans in the Peruvian Province, supporting prior results on this coastline [21][22][23][24] . However, the size structure abruptly changed along the coast as a consequence of changes in the latitudinal taxonomic composition. These findings agree with previous descriptions [37][38][39] and suggest that species assemblages do conform to the biogeographic province they inhabit, similar to what has been described in many other SE Pacific marine invertebrate species such as crustaceans, molluscs and echinoderms 12,15,16 .   www.nature.com/scientificreports www.nature.com/scientificreports/ would be associated with the low diversity of seaweeds and the scarce rocky ecosystems 42 . It is well known that chitons prefer hard substrata on rocky shores, where some species may be found in high abundances 29,37 .
As shown here, nearby localities tended to have similar species, with composition varying significantly as the distance between localities increased. This same pattern has been observed in subtidal and intertidal fishes 42 as well as in benthic communities [21][22][23][43][44][45][46] . In this context, Beta diversity (in addition to have the component of species replacement) can also be evaluated from the perspective of similarity between habitats or species' nestedness 47,48 . In the case of the SEP chitons, replacement did not show a latitudinal pattern since these species rarely have small geographical ranges, and dissimilarity between close localities remained low through the latitudinal gradient. The highest species replacement occurred towards the northern part of Central Chile ecoregion  www.nature.com/scientificreports www.nature.com/scientificreports/ (Huasco and Coquimbo) and the southern end of the Peruvian Province (Valdivia), being evidenced as a direct relationship between geographic distance and Beta diversity.
The taxonomic composition variation of chiton assemblages along the sampled localities was reflected in the faunistic groups supported by the nMDS analysis. Variation in the assemblages' composition was associated to temperature and salinity gradient along the Peruvian Province. Also, the concordance of three faunal groups with some biogeographic breaks was evident, similar to what has been documented in other coastal invertebrates 27 . Some chiton species showed strong endemism, being unique in the Humboldtian ecoregion (e.g. Callistochiton pulchellus, Chaetopleura hennahi, Ischnochiton punctulatissimus, Tonicia swainsoni), while others were exclusive to the Araucanian ecoregion (e.g. Chaetopleura benaventei, Gallardoia valdiviensis). Those species strongly contributed to the differences in species composition. Contrary to such patterns, some chiton species also inhabited different ecoregions (e.g. Acanthoplerura echinata, Chiton granosus, C. cumingsii, T. calbucensis and Chaetopleura peruviana) maintaining the low levels of replacement observed among localities.
For most intertidal invertebrates along the SEP, sea surface temperature and salinity have been shown to play important roles 15,20,21,49 , contrasting to our findings where no significant correlations with overall chiton diversity where observed. The Humboldt ecosystem is characterized by coastal upwelling in Peru and Chile, and is tightly related to the spatial and temporal dynamics of coastal environments 15,50 . In Chile, two main regions of coastal upwelling have been detected, one in the Mejillones Peninsula (northern Antofagasta; ~23°S) and the other between Punta Lengua de Vaca and Punta Pájaros (southern Coquimbo; ~30°S) 15,50 . Accordingly, these sites evidenced the higher richness and diversity of chitons found between Antofagasta and Valparaiso, localities exhibiting maximum levels of upwelling and productivity. However, chitons diversity was not correlated with chlorophyll-a (an indicator of upwelling events in the HCS) but would be more related to the local availability of favorable habitats (i.e. rocky shores) and food 42 , rather than to environmental variables. Along the Humboldt Current System, the sandy beaches are also frequent 15 , representing unsuitable habitat for chitons. These heterogeneous stretches of beaches, combined with rocky shores, could explain the heterogeneous pattern and the absence of a cline of diversity of chitons found in this study. In fact, it has been suggested that local and small-scale processes could be determinants of the diversity and composition of intertidal communities in Chile, at least between Arica and Valparaíso 23 .
Chiton body lengths showed a high spatial variability and a latitudinal cline in assamblage structure. In northern localities of the Humboldtian and Central Chile ecoregions, assamblages were composed by large and small species, but in southern localities of the Araucanian ecoregion assamblages were composed only by medium and small species. These changes in body size structure were related to changes in species composition. The largest species (>100 mm; E. niger and A. echinata, Table S2) inhabit from Coquimbo and Talcahuano through northern localities, respectively 39 . Intertidal chitons are generally thought to be larger than subtidal and deep-sea species 51 , and Southeastren Pacific chitons larger than closely-related species from other ecosystems elsewere, excepting for North Pacific 52 . The evolutionary colonization of the intertidal by chitons increased the competition among them and among other movile invertebrates such as limpets, needing bigger sizes to competitively exclude interspecifics 52,53 . In central Chile, the largest chitons inhabit zones of high energy on the lower intertidal, while small and juveniles coexist in protected zones in high intertidal zones 37 . This pattern suggests certain kind of competition or niche partitioning among chitons and among other invertebrates, since species of similar body size might not be able to coexist because they overlap too much in the use of shared resources 11 . However, it should be noted that different species of chitons, are known to have different diets and niches. For example, there is evidence to www.nature.com/scientificreports www.nature.com/scientificreports/ suggest that the radiation of large bodied Mopalia species in the NE Pacific was driven by niche partitioning to avoid direct competition 54 . Our results of body size ratio and minimum segment length overlap agreed with this hypothesis, since were different from body size simulated data in absence of competition. Another evidence of competition among chitons comes from diet, as several Chilean chitons overlap their diet resources among species and between other intertidal grazers, suggesting a strong competition (e.g. Fissurella spp. and Scurria spp. 33 ).  www.nature.com/scientificreports www.nature.com/scientificreports/ Overall, this work represents the first faunistic and ecological description of polyplacophoran assemblages inhabiting along the Humboldt Current System, providing critical information for an adequate characterization of biodiversity, conservation and potential management strategies in this area, as 3 of the 21 species recorded correspond to artisanally-exploited resources in Chile and Peru.

Methods
Polyplacophoran samplings were carried out in 12 localities during 2013-2016, between Callao (Peru; 12°S) and Valdivia (Chile; 39°S) ( Fig. 1a), encompassing more than 2,500 km of coastline including three of the four previously mentioned ecoregions of the Peruvian Province (i.e. Humboldtian, Central Chile and Araucanian). Chitons were collected during low tide for a period of one to three hours, sampling on average 4 days per locality. Intertidal species were sampled by 1 to 3 people, while subtidal species by 1 to 2 divers through snorkelling at depths up to 5 m. The collected species were relaxed with 5% ethanol with sea water, flattened and preserved in 96% ethanol for further genetic analyses. Once in the laboratory, the total length (mm, including the girdle) of each chiton was measured to perform comparissons among species, ecoregions and years. In order to standardize the sampling effort, a relative abundance (RA) index was estimated, corrected by sampling time and number of people (RA = n°chitons/hour/person). Species identification was performed visually under stereomicroscopy using descriptions and taxonomic keys [55][56][57][58] .
To determine whether the RA matrix of species richness is representative of local diversity, a rarefaction curve 59 was calculated for each locality. The expected number of species was estimated using the Chao-1 algorithm 59 , based on relative abundance.
In each locality, species richness (S) and evenness (J) were estimated using the RA data. To evaluate the replacement degree in species composition between localities and ecoregions, beta diversity was estimated using the Williams index 60 defined as changes species composition with distance 61 . Correlation between geographical distance (km) and beta diversity was determined using the Mantel test with 5,000 permutations 62 .
Similarity between chiton assemblages was evaluated using a non-metric Multidimensinal Scaling (nMDS) 63 performed with presence-absence data through Jaccard similarities. Differences in diversity variables and species composition between ecoregions were evaluated with one-way PERMANOVA using 10,000 permutations 64 . Subsequently, to determine which species differed significantly between-groups, a Percentage Similarity Analysis (SIMPER) was used 63 . Environmental variables such as chlorophyll-a (as a proxy of productivity [mg/m 3 ]), sea surface temperature [SST, °C] and salinity [PSU] were used to test their relationship with different diversity variables and composition using ordinary linear regressions (OLS) and spatial autoregressive models (SAR), both performed in the software SAM v.4.0 51 . These environmental data (monthly averages 2000-2014) were obtained from the Bio-Oracle v2.0 database (http://www.oracle.ugent.be/) through the DIVA-GIS v.7.5 software 65 . Rarefaction curves, Diversity, Mantel test, nMDS, SIMPER and PERMANOVA analyses were performed in PAST v.3.22 statistical package 66 .
Chiton body length structure was compared between species, ecoregions and years using Generalized Linear Model (GLM) implemented in the statistical software R v 3.6.1 67 In order to explore the possibility of competition, we estimated body size overlap between species in each locality, we generated 5,000 random matrices to create a random size assemblage in the software ECOSIM v7.72 68 . As size overlap metric, we used the minimum segment length (calculated as the difference in body sizes of two adjacent species) and absolute size differences for each locality. The minimum segment length corresponded to Hutchinson's hypothesis, which stated that a minimum spacing between species is necessary for coexistence. The standardized effect size (SES) of the body size ratio, was calculated as: observed index -mean (simmulated indices)/standard deviation (simulated indices). If SES was greater than 2 or less than -2, it was considered as statistically significant with a probability <0.05 68,69 .