Body size variability across habitats in the Brachionus plicatilis cryptic species complex

The body size response to temperature is one of the most recognizable but still poorly understood ecological phenomena. Other covarying environmental factors are frequently invoked as either affecting the strength of that response or even driving this pattern. We tested the body size response in five species representing the Brachionus plicatilis cryptic species complex, inhabiting 10 brackish ponds with different environmental characteristics. Principal Component Analysis selected salinity and oxygen concentration as the most important factors, while temperature and pH were less influential in explaining variation of limnological parameters. Path analysis showed a positive interclonal effect of pH on body size. At the interspecific level, the size response was species- and factor-dependent. Under the lack of a natural thermo-oxygenic relationship, the negative response of size to temperature, expected according to ‘size-to-temperature response’ rules, disappeared, but a positive response of size to oxygen, expected according to predictions selecting oxygen as a factor actually driving these rules, remained. Our results confirm the crucial role of oxygen in determining the size-to-temperature patterns observed in the field.

www.nature.com/scientificreports/ sister species in studies combining the selective and plastic responses to unravel the evolutionary processes. In this study, we examined inter-and intraspecific variability in body size in the Brachionus plicatilis (Rotifera) cryptic species complex in relation to environmental conditions in natural habitats. Comparison of the patterns at these two, inter-and intraspecific levels allows to distinguish between the genetic and phenotypic response. Rotifera is an especially interesting group in this regard because at least 42 species complexes of rotifers have been discovered 41 . Among these groups, the best known is the Brachionus plicatilis cryptic species complex. Currently, 15 Brachionus species have been recognized using molecular techniques, six of them have been formally described 42 , and four of the named ones are known to inhabit ponds in eastern and central Spain. The sympatric coexistence of these species in the well-documented system of brackish ponds in Spain is mediated by seasonal ecological specialization 43 related to factors such as salinity, temperature, resource use and vulnerability to predation as reviewed in 41,44 . A phylogenetic analysis showed signatures of coexistence in this region extending back to the Pleistocene 45 . According to paleolimnological studies, this pattern has persisted for several decades in single localities, at least in two of the large species in the complex (B. plicatilis sensu stricto and B. manjavacas) 46 . According to current knowledge, between-species gene flow is absent in this complex in the wild 47 . It is also important to mention that the members of the B. plicatilis species complex are all herbivorous and feed on algae with no specialization to specific type. The first indication that B. plicatilis was a cryptic species complex was the observation of three apparent size classes, initially referred to as the "L" (large), "SM" (medium) and "SS" (small) morphotypes 42,44 . Therefore, body size divergence and speciation are linked in the B. plicatilis species complex, but the causal relationship remains unknown. Regarding the response of size to experimental temperature, a size decrease with increasing temperature has been observed at the intra-and interspecific (congeneric) levels [48][49][50] . Additionally, an association between small-sized species and high temperature occurs in the wild at interspecific (congeneric) levels 49,51 . Moreover, comparison across three species from the B. plicatilis complex showed that species size affects the thermal dependence of diapause egg hatching 52 . Finally, B. plicatilis sensu stricto (the species that gives name to the complex) genetically adapts to low or high temperature relatively quickly through body (and egg) size adjustment 53 , revealing the crucial importance of temperature in the species' life history, which is consistent between its phenotypic and genetic background 53 .
Sampling methods for the present study took advantage of sediment egg banks, which were recently proposed as a newly emerging field of resurrection ecology 54 and provide a powerful approach for disentangling the plastic and genetic effects in adaptive evolution 55 . Being a cyclical parthenogen, the species in the B. plicatilis complex produce resting (diapausing) eggs which are a dormant, resistant stage in the bouts of sexual reproduction following periods of asexual proliferation 44 . Clones established from resting eggs deposited in the sediments of brackish ponds situated in eastern and central Spain were shown to belong to five species in the B. plicatilis complex differing in size.
Adopting an approach that stresses its evolutionary significance, our study aimed to reveal the role of environmental conditions in determining the body size response. Based on the extensive database of environmental conditions for 25 ponds in Spain-localities with the potential of harbor B. plicatilis species complex-, we selected 10 ponds representing the transect across the most important environmental parameters according to PCA analysis. We measured individuals of clones belonging to all five species, established from resting egg each, and we related the rotifers' body size to environmental conditions of the pond of clone origin. This approach enabled us to compare the environment-dependent size differences among species and within species (Fig. 1).

Results
Limnological parameters and population selection. When applied to the limnological parameters recorded in ponds of Central and Eastern Spain, including those where our specimens were collected, PCA showed that 35.7% of the variance was explained by the first principal component (PC1), 27.6% by PC2, 20% by PC3 and 17% by PC4. The limnological parameters associated with PC1, PC2 and PC3 were oxygen concentration, salinity, and temperature, respectively ( Fig. 2A,B). pH was of secondary importance in PC1 and PC3. The mean scores for the 25 ponds are presented in Fig. 2C. From among them, we selected 10 ponds representing a transect across the parameters which drove the first two PCs. These were the ponds from which the sediment was collected to establish the clonal populations of Brachionus rotifers.

Populations and clones.
In total, we measured 4150 females carrying a single egg and belonging to 178 clones, with 22 ± 6 individuals being measured per clone on average. In the cases where clones were not numerous enough, to increase the sample size the measurements were supplemented with few females carrying two eggs. The highest proportion of such females per clone was 4 out of 25. The distribution of the clones across the species is provided in Table 1. The most frequent species in the studied system were B. ibericus and B. plicatilis, each of which was present in seven ponds, and the least common was B. rotundiformis, which was found in four ponds (Table 1). In the majority of ponds, two or three species were collected, with Hondo Sur showing the highest recorded species richness (Table 1).

Path analysis.
Causal relationships between limnological parameters and body size were inferred using path analysis. From an initial (a priori) model (Fig. 3A), in the first step the following paths were removed: Temperature to Salinity (p = 0.99 for the linear relationship estimation), Temperature to Body size (p = 0.81), and Temperature to Oxygen concentration (p = 0.20). The output was referred to as the Version 1 model. In the next step, we removed the following paths: Oxygen to Body size (p = 0.17) and Salinity to Body size (p = 0.10), thus generating the Version 2 model. The goodness-of-fit indices were slightly worse for Version 2 than for the Version 1 model (Table 2). Therefore, we provided the final path coefficients and drew inferences for the Version 1 www.nature.com/scientificreports/ model (Fig. 3B, Table 2). The qualitative results did not differ between the three models. All indices obtained for the Version 1 model indicated an overall good model fit 56 . According to the results, (i) the body size of Brachionus sp. clones increased with increasing pH and was not affected by any other parameter; (ii) the oxygen concentration was positively affected by pH and salinity, with no effect of temperature; and (iii) pH was negatively affected by temperature and positively affected by salinity (Fig. 3B, Table 2B). The R 2 value was 0.10 for body size, 0.32 for the oxygen concentration and 0.40 for pH.  Fig. 4A, together with the pond environment described by three limnological parameters: temperature, oxygen concentration and salinity (Fig. 4B). Body size was species and parameter dependent. B. ibericus, B. manjavacas and B. plicatilis were larger in the presence of higher oxygen concentrations, and signatures of the same trend were observed for SM-X (Table 3). B. manjavacas, B. plicatilis and B. rotundiformis were larger in the presence of higher temperatures, whereas the opposite relationship was identified for SM-X. B. ibericus, B. manjavacas and SM-X which were larger in the presence of a higher pH (Table 3). No species showed a correlation between its body size and salinity (data not shown in Table 3). The regression plots for each species-parameter combination are presented in the supplementary materials (Fig. S2).

Discussion
Based on long-term data describing the environmental characteristics of 25 shallow, brackish ponds in central and eastern Spain, we found that oxygen and salinity were the most important environmental factors, followed by pH and temperature. Five Brachionus rotifer species differed in body size across the ponds, with considerable variability among populations of the same species and clones of the same population. This suggests that the egg banks that we examined covered a wide range of variability in pond-specific environments. These results are in accordance with previous molecular marker and life history analyses revealing high population differentiation 45,57 and high within-population genetic diversity in species of the B. plicatilis complex 58 . Notably, the size variation found in our study was not due to the laboratory conditions during development, as individuals were measured  Table 1 for the acronyms). www.nature.com/scientificreports/ under standard laboratory conditions, but it was associated to differences between the ponds of clones origin. Size at maturity at the interclonal level (regardless species) correlated with pH, with larger clones being observed at a higher pH. Path analysis did not reveal a size response to the thermo-oxygenic conditions or salinity. A different pattern was unmasked when body size was analyzed separately for the studied species. In this analysis, three species exhibited larger individuals in better-oxygenated ponds, and one additional species exhibited signatures  Table 1). The arrows show the direction of the increase in the oxygen concentration (horizontal) and salinity (vertical). www.nature.com/scientificreports/  Table 2). Arrow thickness indicates the importance of a given path; path and covariance coefficients are provided when significant. *** p < 0.0001, ns -no significant difference. www.nature.com/scientificreports/ of this pattern (p = 0.052). In three species body size was larger, while one species had smaller clones in warmer ponds. Finally, three species had larger individuals in the ponds with a higher pH, while none of the species showed a change in size associated with salinity. The positive response of size to the oxygen level supports the general prediction of small size being an adaptation to low oxygen availability, while the positive response to temperature contradicts the theory 9 and the considerable amount of empirical evidence of a decrease in body size with increasing temperature. The patterns that we found are likely driven by many abiotic and biotic factors that interact in the examined shallow, brackish water bodies, which we will discuss below.
The effect of salinity on pH and oxygen. Salinity variation is very important in our study system, as shown by the PCA of the 25 ponds and its high relative variability across the 10 ponds in which rotifer populations were sampled (CV = 88%). Salinity affected the pH and the oxygen concentration. However, our results show no correlation of the body size of Brachionus species and salinity, confirming previous findings 50 , although this variable has been reported a key factor in the ecological specialization of Brachionus populations 57,59 . The dependence of salinity on temperature would mean that the concentrations of all the ionic constituents dissolved in pond water would increase with desiccation caused by high temperature, as observed under summer conditions. The lack of a temperature effect on salinity led us to conclude that idiosyncratic differences in salinity between the ponds resulted from geological and geographic characteristics (e.g., proximity to the sea; 60 ) and  www.nature.com/scientificreports/ played a dominant role over seasonal changes. On the other hand, it is generally assumed that in aquatic systems, pH and oxygen concentrations are correlated with photosynthesis 60,61 . This might explain why both the PCA (the effect explained by PC1) and path analysis showed a positive relationship between pH and oxygen concentration. Additionally, high photosynthetic activity indicates a high microalgal density, meaning that more food is available to rotifers. This might explain the positive effect of pH on the body size of Brachionus rotifers. Speciation and adaptive radiation caused by trophic conditions is known for the system of African cichlids 39 .
Thermo-oxygenic environment. The oxygen concentration was the second most variable parameter, and temperature was the third across the investigated ponds. The overall mean oxygen concentration in the studied ponds was relatively high. Path analysis showed no effect of temperature on the oxygen concentration. This unexpected result could also be observed in a bivariate regression analysis, as temperature was unrelated to oxygen (slope = -0.06, R 2 = 0.004, p = 0.86), contrary to the general pattern observed in aquatic systems 60,62 . Such a negative relationship is a necessary condition for size adjustment to be responsive with regard to temperature and anticipatory with regard to oxygen 37 . The lack of a relationship between temperature and oxygen seems to be sufficient to explain the lack of an influence of temperature or oxygen on rotifer body size that was found when individuals of the different species were merged. The importance of the link between temperature and the oxygen concentration in modulating species size was previously shown for Icelandic diatoms 26 .
Other possible influential abiotic factors. The genetic diversity of the egg banks of species of the B.
plicatilis complex was previously found to depend on the pond area 63 , while the species life history strategy is affected by environmental unpredictability 64,65 ; specifically, interannual fluctuations in the length of the planktonic growing season 66 . The size response to these conditions has not been systematically tested so far. We show some patterns of species-specific environmental preferences. For example, SM-X is absent at high salinity, and B. plicatilis and B. ibericus are euryoic (but see below), while the occurrence of the smallest species, B. rotundiformis, is limited to ponds with an average temperature higher than 17 °C. This uneven species distribution may be responsible for the pattern that we found in the path analysis, showing the most apparent dependence of body size on pH, which is a proxy for resource availability.

Sampling effects and species distribution.
When data analyses were performed separately for each species, the two largest species, B. plicatilis and B. manjavacas, were found to share a similar response of exhibiting a larger body size in warmer and betteroxygenated ponds. One difference between these species was that B. manjavacas also was larger in ponds with a higher pH than in ponds with a lower pH. A common feature of the two SM species was positive size dependence on pH (possibly due to the direct relationship of pH with the food of these herbivorous rotifers), but size of B. ibericus increased with O 2 conditions, and not changed with temperature; SM-X presented a similar tendency in regard to oxygen but was smaller in warmer ponds. The smallest species responded positively in size solely to temperature. Such a pattern suggests that speciation within the B. plicatilis species complex has been driven to some extent by diverging tolerance to crucial environmental factors. It was previously found that the response of the largest B. plicatilis species to temperature is affected by salinity 48,69 and that temperature-driven egg size adjustment occurs only at intermediate salinities 48 , while unidentified SM species showed an inverse response of size to temperature regardless of salinity conditions 48 . Therefore, speciation toward the "L" and "SM" morphotypes could have been associated with differential vulnerability to temperature and salinity. Our results are in line with this speculation.

Body size and species interactions.
The species in the B. plicatilis complex exhibit somewhat different but overlapping niches, resulting in an overlap of their seasonal distributions if they co-occur in a locality. Niche differentiation implies temperature and salinity differences 47,51 , which allow seasonal succession and the partitioning of resources 43 . The most interesting case is, however, the coexistence of B. plicatilis and B. manjavacas. These two largest Brachionus species from the complex are not distinguishable morphologically 70 and differ in size by only 6% 71 . Their ecological niches are intriguingly close. Empirical evidence suggests that under environmental fluctuations involving salinity, differential adaptation to salinity and divergence in life history traits associated with different levels of opportunistic strategies are relevant to their coexistence 69 . The results of our study indicate a new candidate for a crucial parameter affecting ecological divergence: oxygen availability. Both species adjust their size to oxygen conditions, but this relationship is steeper and stronger in B. manjavacas. Interestingly, this species had the smallest body size in the least-oxygenated pond (HOS) and its largest body size in the best-oxygenated ponds (ONT and CVF). These observations may indicate high sensitivity of B. manjavacas to oxygen availability.
The temperature-oxygen effect on body size. Among the five Brachionus species collected from 10 ponds, we found only one species that exhibited a smaller body size in warmer ponds, consistent with the expected size-to-temperature response. Moreover, this species, SM-X was not found in ponds with high salinity. www.nature.com/scientificreports/ The considerably variable salinity conditions could affect the response of size to temperature in the other more euryhaline species. However, the most likely factor responsible for the reversal of the response of size to temperature in our study system is the lack of the common, expected negative temperature-oxygen relationship. The increase in body size with increased nutrition (availability of algae; pH as a proxy in our case), which was the most apparent result of our path analysis, could mask the possible pattern of decreasing size with increasing temperature, rather than explaining the generally reversed pattern that we found. As noted elsewhere, "the stronger effect of nutrition than of oxygen may be observed when temperature-dependence of the former is steeper than that of the latter" 26 . Therefore, in our study system, salinity (indirectly) and pH (directly) affected the response of size to temperature, causing the absence of the expected decrease with increasing temperature. Intriguingly, under the presence of interfering factors, the rotifer species exhibited smaller sizes at lower oxygen levels, as predicted by the theory, confirming the crucial role of the oxygen concentration in driving body size patterns 26 . The only exception to the "when there is less oxygen, grow smaller" response was observed for the smallest species, B. rotundiformis. However, as mentioned above, this species occurred only in warm waters, and possibly because of its small size, it is equipped with physiological mechanisms for dealing with hypoxia, making the body size plasticity found in other species unnecessary. The relationship of the strength of the response of size to temperature with the level of species thermal specialization was noted in another study on three Brachionus species 49 , in which only B. plicatilis, the largest species, which is euryoic, showed a phenotypic size decrease with increasing temperature, while two other less-temperature-tolerant species, B. ibericus and B. rotundiformis, showed no such pattern. Although we found B. ibericus to be the most frequently occurring species across ponds, this species shows thermal specialization by occurring in the water column within the narrowest time window during the year 51 . Previous studies confirm that our results do not violate the general size-to-temperature rules for Brachionus rotifers. Clear size-dependent temperature preferences were shown in three Brachionus species originating from the same pond system: (i) smaller species generally presented a higher optimal temperature in relation to the population growth rate 49 ; (ii) smaller species preferred higher temperatures for hatching from resting eggs 52 ; and (iii) Brachionus plicatilis s. s. decreased in size with increasing temperature, which was reflected at the levels of both short-term phenotypic, nongenetic plasticity 49,50 and genetics 53 .
In the previous studies, the joint contribution of selection and plasticity as affecting the morphological traits was examined for systems with high gene flow 3 , while in this study we unravel the respective patterns for the system with limited gene flow, referring to body size, the organismal trait which informs about morphology, physiology and life strategy.

Limnological parameters. The system of brackish ponds in eastern and central Spain inhabited by the
Brachionus plicatilis species complex has been sampled since the 1990s for research projects performed at the Cavanilles Institute of Biodiversity and Evolutionary Biology (ICBiBE; University of Valencia). A total of 44 ponds were seasonally inspected, and major limnological parameters (temperature, conductivity, salinity, oxygen concentration, oxygen saturation and pH) were recorded 43,45,46,51,65,[71][72][73][74][75] . Among the sampled ponds, 25 ponds were inspected at least three times in distinct periods of the year. Based on the resulting database for these ponds, Principal Component Analysis (PCA) was conducted for temperature, oxygen concentration, pH and salinity using Canoco 5.0 76 . Conductivity and oxygen saturation were excluded from the analysis because of their very high correlations with salinity and oxygen concentrations, respectively. The mean parameter values for the 25 ponds are presented in the supplementary materials (Table S1). The PCA results were used to select the ponds from which the rotifers were isolated and studied.

Establishment of clones and body measurements.
From the 25 ponds referred to above, we chose 10 ponds (Fig. 1) to study the corresponding populations of the B. plicatilis species complex. According to the PCA, the selected ponds represented a gradient across the two first principal components (Fig. 2). Clones of the species were established from hatched resting eggs collected from the pond sediments. We used ponds' sediments that were either previously obtained by ICBiBE or freshly collected in the field for this study. For two ponds, we established clones from both individuals collected in the water column and resting eggs isolated from the sediment. The details are provided together with the statistics of the pond limnological parameters in Table 1. Resting eggs were obtained from 30 g sediment samples using a modified sucrose flotation technique 58 . Hatching was induced at a salinity of 6 ppt in Petri dishes exposed to light at 25 °C. Hatchlings were individually transferred from Petri dishes to the wells of a 24-well plate with 1 mL of microalgal suspension containing approximately 3 × 10 5 cells/mL of Tetraselmis suecica as food (under the same culture conditions as for resting egg hatching). Density of microalgae was estimated using an automated cell counter based on image analysis (Celeromics Technologies; Valencia, Spain). After clonal proliferation (25 °C, 12 ppt salinity, continuous light of approximately 75 µmol quanta m -2 s -1 ), individuals were fixed with 40 µL of Lugol solution for size measurements. B. plicatilis and B. manjavacas were identified via the PCR-RFLP technique 71 in a subsample that was not fixed with Lugol. B. ibericus and B. rotundiformis were identified according to spine morphology 77 . A fifth morphotype was morphologically a medium-sized Brachionus species (i.e., SM clade; 42 in which females carried resting eggs inside their body (supplementary materials, Fig. S1). However, it was easily distinguishable from B. ibericus because of its larger size and distinct spine morphology. The spines were of blunt triangle shape and relatively large, while in B. ibericus they are clearly sharper and smaller. Most likely, this morphotype was one of the not yet formally named species B. ' Almenara' or B. 'Tiscar' , both of which are known to occur in the examined study area 46,78  www.nature.com/scientificreports/ Table 1). These clones were included in the path analysis but not in the analyses in which the species factor was involved (see below). For size measurements, several female rotifers from each clone carrying a single egg were photographed under 20 × magnification using a Nikon Eclipse E800 microscope equipped with a Nikon DS-Ri1 camera, assisted by NIS-Elements BR software. The length and width of the lorica, the external cuticle, were measured individually in ImageJ 1.46r software (NIH, USA), and the product of these measurements served as a body size estimate (in µm 2 ); this method has been applied previously to estimate the size of Brachionus 49 and Lecane 37 rotifers. Females carrying a single egg were assumed to have matured very recently.

Path analysis.
To test for the dependence of rotifer body size on limnological parameters, a path analysis was conducted on the measurements of a total of 178 rotifer clones, without reference to species. We used PROC CALIS 79 with the maximum likelihood method of coefficient estimation based on the variance-covariance matrix 56 . First, we assumed an a priori model (referred to as the Initial model; Fig. 3A) based on common limnological knowledge 60 . This model accounts for the direct effects of pH, oxygen concentrations and salinity on body size and for the direct and indirect effects of temperature. Therefore, the first four factors were endogenous variables in our model, while temperature was an exogenous variable (Fig. 3A). To assess the model, we followed the rules recommended by O'Rourke and Hatcher 56 , removing the least statistically significant paths identified by the Wald test when the path-analysis model did not show proper goodness-of-fit for any of the most important indices, which are provided in Table 2A. This model selection process is a stepwise process and allowed us to achieve the reduction of the initial model and to compute the importance of each path.
Species-specific body size variation. The difference in body size between ponds was tested at the intraand interspecific levels using a generalized linear mixed model (Method = REML) in PROC MIXED (SAS, 2013). The model included 'pond' and 'species' as fixed factors and 'clone' , nested in the species and pond combination, as a random factor. Additionally, PROC REG 79 was used for bivariate linear regression analysis to assess the dependence of body size on temperature, the oxygen concentration, pH and salinity, addressing each species separately, with clonal mean measures as the input.

Data availability
The data are available from an open repository of the Jagiellonian University at a https:// doi. org/ 10. 26106/ h0a7-3c58.