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.

Introduction

Understanding the relationship between features of living individuals and their abiotic environment is a fundamental step in understanding ecology. However, this scientific issue remains elusive and has been quite recently even cited as one of the grand challenges in biology[1]. Special interest is expressed toward the combination of selection and phenotypic plasticity in understanding the organismal response to changing environment[2-4].

Temperature is regarded as the most important abiotic factor in ecological niches and, consequently, species geographical ranges. It affects organisms at all levels of life organization, ranging from physical thermal effect on biological processes[5,6], to ecological trophic and other interspecific relationships[7,8]. Organisms respond to temperature in many different ways. Among these responses, the most easily recognizable and crucial for life strategies is the body size response. Size decrease with increasing temperature is observed at the phenotypic level[9,10] and reflected in the genetic background[11-13]. However, this pattern is ecologically puzzling (e.g.,[14]), and the underlying proximate and ultimate mechanisms are not conceptually easy to test[15]. The difficulty arises from the occurrence of environmental factors correlated with temperature[16], which may or may not be temperature induced. The factors that are noted in the literature as interfering with the general organismal response to temperature, which differs across habitats, are precipitation[17], resource availability[18,19], seasonality[20,21], generations number[22], the ability to disperse[23] and oxygen availability[24-27]. The last is of special interest. A decline in oxygen availability at high temperature, especially in aquatic systems[28], has been suggested to be a selective factor driving decreases in body size through cell shrinkage to meet aerobic metabolic demands under a diminishing oxygen supply[29,30]. There are a number of studies confirming the size response to temperature-dependent oxygen conditions, either at the long-term genetic level or the short-term plastic phenotypic level, both indirectly[24,25,31-34] and directly[27,35,36]. However, the empirical evidence confirming the adaptive significance of this pattern is very limited[37]. The temperature-size rule (TSR) is a widespread phenomenon describing the phenotypic plastic response of size to temperature[9]. In this perspective, TSR could be explained in terms of anticipatory vs. responsive plasticity toward oxygen and temperature, respectively, following the reasoning of Whitman and Agrawal[38]. Provided that the role of oxygen in the response of size to temperature is correctly predicted, the question arises: what is the actual cue for the size response.

Cryptic species complexes provide a promising study system for ecological hypothesis testing because the high similarity of the species in such complexes simplifies the inference of the patterns and processes involved in evolutionary ecology. Cichlids[39] and sticklebacks[40] are the most pronounced examples of the potential of the 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-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).

Figure 1

Map of Spain with the locations of the 10 ponds from which the populations of the Brachionus plicatilis species complex were sampled (see Table 1 for the acronyms).

Table 1

Mean values of the limnological parameters (collected at least three times in distinct periods of the year) of 10 ponds from which the rotifer populations were studied, selected after PCA. Sampling details are shown. Numbers after each species name are the number of clones studied.

Pond	Acronym	N	Temp (°C)	Salinity (g/L)	O2 conc. (mg/L)	pH	Sediment collection date	Species isolated (number of clones)
El Bassaret de L’Altet	ALT	7	19.4	6.7	7.5	7.8	05/10/2017	B. ibericus (13)*, B. plicatilis (6), SM-X (7)*
Camino de Villafranca	CVF	5	20.7	48.8	10.4	9.2	22/05/2013	B. ibericus (10), B. manjavacas (10)
Clot de Galvany	CdG	5	18.1	14.1	9.3	8.3	05/10/2017	B. ibericus (9), B. plicatilis (8)*
Hondo Norte	HON	21	20.3	10.3	4.7	7.8	05/10/2017	B. ibericus (6), B. rotundiformis (2), SM-X (9)
Hondo Sur	HOS	24	20.6	10.8	4.4	7.8	05/10/2017	B. ibericus (7), B. manjavacas (1), B. plicatilis (3), B. rotundiformis (2), SM-X (2)
Ontalafía	ONT	3	14.2	5.4	9.3	8.6	22/05/2013	B. ibericus (7), B. manjavacas (2), SM-X (5)
Pétrola	PET	15	15.7	52.7	7.7	8.3	14/09/2017	B. ibericus (10), B. manjavacas (8), B. plicatilis (1)
Torreblanca Poza Norte	TON	38	17.8	10.0	5.0	7.6	06/06/2017	B. plicatilis (10), B. rotundiformis (6)
Torreblanca Poza Sur	TOS	106	18.6	25.3	10.9	8.0	06/06/2017	B. plicatilis (4), B. rotundiformis (12)
Salobrejo	SAL	12	14.3	12.9	5.3	8.9	24/09/2013	B. manjavacas (10), B. plicatilis (3)

*Sampled from the water column.

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.

Figure 2

PCA of the major limnological parameters of 25 brackish ponds in eastern-central Spain where the occurrence of the Brachionus plicatilis species complex has been reported. (A) – PC1-PC2 scores (the 10 ponds selected to study rotifer populations are individually marked); (B) – factor loadings for PC1-PC4, where the highest load for each PC is indicated in bold; (C) – pond-specific loadings (mean ± SD) for the 10 ponds where rotifer populations were sampled (acronyms in Table 1). The arrows show the direction of the increase in the oxygen concentration (horizontal) and salinity (vertical).

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 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 R2 value was 0.10 for body size, 0.32 for the oxygen concentration and 0.40 for pH.

Figure 3

Model for the path analysis of limnological parameters affecting rotifer body size. (A) – the initial (a priori) model, E – endogenous variables; (B) – path correlations with the selected model (Version 1; 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.

Table 2

Path analysis of rotifer clone (N = 178) body size as related to the limnological parameters of 10 ponds. A – goodness-of-fit indices chosen according to recommendations of O'Rourke and Hatcher[56]; B – standardized path coefficients for the Version 1 model.

Model	χ2	df	P > χ2	CFI	SRMR	RMSEA	(RMSEA CL90)
A
Initial	1.14	1	0.28	1.00	0.014	0.029	(0.000–0.204)
Version 1	2.85	4	0.58	1.00	0.026	0.000	(0.000–0.097)
Version 2	7.85	6	0.25	0.99	0.047	0.042	(0.000–0.112)

Note: Initial (a priori) model as illustrated in Fig. 3A; Version 1 – the model without the paths of Temperature to Salinity, Temperature to Body size and Temperature to Oxygen concentration; Version 2 – reduction of Version 1 by removing the paths of Oxygen concentration to Body size and Salinity to Body size.

χ2 – chi-square; df – degrees of freedom; CFI – Comparative Fit Index; SRMR – Standardized Root Mean Square Residual; RMSEA – Root Mean Square Error of Approximation; RMSEA CL90 – RMSEA 90% Confidence Limits.

Species-specific body size variation

The mean body sizes of the species at maturity (from the smallest to the largest species) were 15 700 ± 100 for B. rotundiformis, 19 700 ± 100 for B. ibericus, 26 200 ± 200 for SM-X, 44 800 ± 400 for B. manjavacas and 50 800 ± 400 for B. plicatilis (µm2 ± SE). The GLM analysis showed that body size differed across species (F(4,146) = 394.44; p < 0.0001) and ponds (F(9,146) = 15.84; p < 0.0001), their interaction (F(13,146) = 2.96; p = 0.0007), and clones (i.e., within-population genetic variation; Z-value = 7.84; p < 0.0001). The variation in body size among all five species across the studied ponds is shown in 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).

Figure 4

Illustration of the environmental background for body size performance across five species from the Brachionus plicatilis cryptic species complex from Spain. (A) – mean computed by averaging the clonal mean values (least square means ± SE) of the body size of five species in 10 studied ponds. (B) – limnological parameters of the ponds (ordered as in A) described by temperature (x-axis), the oxygen concentration (y-axis) and salinity (bubble area).

Table 3

Bivariate linear regression analyses of body size on limnological parameters.

	B. ibericus			B. manjavacas			B. plicatilis			B. rotundiformis			SM-X
	O2	Temp	pH	O2	Temp	pH	O2	Temp	pH	O2	Temp	pH	O2	Temp	pH
Slope	0.02	0.01	0.10	0.06	0.04	0.33	0.04	0.07	− 0.10	− 0.01	0.09	0.12	0.02	− 0.03	0.19
Adj. R2	0.11	0.01	0.13	0.48	0.31	0.48	0.18	0.30	0.02	− 0.03	0.36	− 0.02	0.13	0.24	0.24
p-value	0.004	0.22	0.002	< 0.0001	0.0007	< 0.0001	0.006	0.0004	0.20	0.57	0.002	0.42	0.052	0.011	0.01

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 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 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 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, R2 = 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.

Sampling effects and species distribution

Our study might be affected by incomplete sampling of the species of the B. plicatilis complex in the 10 ponds. For example, previous studies showed the presence of B. ibericus at Torreblanca Poza Sur TOS;[52,67], while B. rotundiformis was previously found in El Bassaret de L’Altet[68], and the “L” morphotype (B. plicatilis or B. manjavacas) was previously known from Hondo Norte[68]. On the other hand, there was no previous record of B. rotundiformis at Hondo Sur. The most likely reason for such differences is the habitat heterogeneity of the sediments in some ponds. Therefore, the resting eggs of some species may be found only at specific sites within a pond sediment. Nevertheless, our results show relatively good coverage of the general distribution of B. plicatilis cryptic species as reviewed by Lapesa[68].

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.

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 better-oxygenated 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 O2 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. 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.

Methods

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-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 × 105 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]. Therefore, we designated this species as SM-X. Five clones (two B. plicatilis or B. manjavacas from Hondo Sur and three B. ibericus or B. rotundiformis from Poza Norte) remained unidentified (not included in 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 µm2); 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.

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.

The difference in body size between ponds was tested at the intra- and 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.

Supplementary Information 1.