Thermal and oxygen conditions during development cause common rough woodlice (Porcellio scaber) to alter the size of their gas-exchange organs

Abstract Terrestrial isopods have evolved pleopodal lungs that provide access to the rich aerial supply of oxygen. However, isopods occupy conditions with wide and unpredictable thermal and oxygen gradients, suggesting that they might have evolved adaptive developmental plasticity in their respiratory organs to help meet metabolic demand over a wide range of oxygen conditions. To explore this plasticity, we conducted an experiment in which we reared common rough woodlice (Porcellio scaber) from eggs to maturation at different temperatures (15 and 22 °C) combined with different oxygen levels (10% and 22% O2). We sampled animals during development (only females) and then examined mature adults (both sexes). We compared woodlice between treatments with respect to the area of their pleopod exopodites (our proxy of lung size) and the shape of Bertalanffy’s equations (our proxy of individual growth curves). Generally, males exhibited larger lungs than females relative to body size. Woodlice also grew relatively fast but achieved a decreased asymptotic body mass in response to warm conditions; the oxygen did not affect growth. Under hypoxia, growing females developed larger lungs compared to under normoxia, but only in the late stage of development. Among mature animals, this effect was present only in males. Woodlice reared under warm conditions had relatively small lungs, in both developing females (the effect was increased in relatively large females) and among mature males and females. Our results demonstrated that woodlice exhibit phenotypic plasticity in their lung size. We suggest that this plasticity helps woodlice equilibrate their gas exchange capacity to differences in the oxygen supply and metabolic demand along environmental temperature and oxygen gradients. The complex pattern of plasticity might indicate the effects of a balance between water conservation and oxygen uptake, which would be especially pronounced in mature females that need to generate an aqueous environment inside their brood pouch.

Terrestrial isopods have evolved pleopodal lungs that provide access to the rich aerial supply of oxygen. However, isopods occupy conditions with wide and unpredictable thermal and oxygen gradients, suggesting that they might have evolved adaptive developmental plasticity in their respiratory organs to help meet metabolic demand over a wide range of oxygen conditions.
To explore this plasticity, we conducted an experiment in which we reared common rough woodlice (Porcellio scaber) from eggs to maturation at different temperatures (15 and 22 � C) combined with different oxygen levels (10% and 22% O 2 ). We sampled animals during development (only females) and then examined mature adults (both sexes). We compared woodlice between treatments with respect to the area of their pleopod exopodites (our proxy of lung size) and the shape of Bertalanffy's equations (our proxy of individual growth curves).
Generally, males exhibited larger lungs than females relative to body size. Woodlice also grew relatively fast but achieved a decreased asymptotic body mass in response to warm conditions; the oxygen did not affect growth. Under hypoxia, growing females developed larger lungs compared to under normoxia, but only in the late stage of development. Among mature animals, this effect was present only in males. Woodlice reared under warm conditions had relatively small lungs, in both developing females (the effect was increased in relatively large females) and among mature males and females.
Our results demonstrated that woodlice exhibit phenotypic plasticity in their lung size. We suggest that this plasticity helps woodlice equilibrate their gas exchange capacity to differences in the oxygen supply and metabolic demand along environmental temperature and oxygen gradients. The complex pattern of plasticity might indicate the effects of a balance between water conservation and oxygen uptake, which would be especially pronounced in mature females that need to generate an aqueous environment inside their brood pouch.
The evolution of gas exchange structures, called pleopodal lungs, was apparently a key innovation among land-colonizing isopods (Csonka et al., 2013). Isopod lungs are formed by pleopod exopodites and are located on the ventral side of the body. Depending on the level of terrestrial specialization, different species of isopods have different lung structures (Hoese, 1982;Hsia et al., 2013;Schmidt and W€ agele, 2001). Species less specialized for dry terrestrial environments have gas exchange surfaces directly exposed to the external environment. In drought-tolerant species, the gas exchange surface is located within lungs that are formed via invagination of the pleopod exopodite surfaces (Babula, 1981;Bielawski and Babula, 1980;Schmidt and W€ agele, 2001) ( Fig. 1 A, C). The lungs limit the direct contact of the respiratory organs with ambient air, which helps to conserve water (Hornung, 2011;Schmidt and W€ agele, 2001). The lungs consist of a large number of branched tubules (pseudotracheae) formed by a thin respiratory epithelium (Bielawski and Babula, 1980;Wright and Ting, 2006). The apical side of the epithelial cells is covered with a cuticle (Bielawski and Babula, 1980). In relatively xeric isopod species, air enters the lungs through one or more openings called spiracles (Hoese, 1982;Hsia et al., 2013;Warburg, 1993). Interestingly, the cuticle that covers cells surrounding the spiracles can create microfolds (Fig. 1A, B) that trap water; these humidify the entering air and serve as additional protection from drought (Babula, 1981). In insects with a tracheal system, oxygen diffuses directly into tissues and cells without mediation by oxygen-binding proteins (Klok et al., 2004;Wright and Ting, 2006). In contrast, terrestrial isopods possess a hemocyanin protein that binds oxygen in the lungs and then transports it via the hemolymph to the tissues (Klok et al., 2004;Wright and Ting, 2006).
Isopods that utilize ambient oxygen might be expected to be less oxygen-limited than aquatic species. However, even terrestrial isopods are likely challenged by poor oxygen (hypoxic) conditions, especially if they inhabit periodically flooded burrows, areas with high microbial respiration, such as decaying wood, habitats covered by snow, intertidal zones or simply high elevations (Harrison et al., 2018a;Hoback and Stanley, 2001;Paim and Beckel, 1964;Wright and Ting, 2006). Like other ectotherms (Verberk et al., 2011), isopods also face a risk of an increasing incongruity between oxygen supply and demand with increasing temperatures. This incongruity was postulated to set thermal limits in ectotherms and might lead to functional hypoxia even in normoxic conditions (Harrison et al., 2018a;P€ ortner, 2001). To address the effects of this imbalance, we studied the common rough woodlouse (Porcellio scaber), a detritivorous terrestrial isopod, which inhabits leaf litter and compost, where it is likely to be challenged with hypoxic conditions (Wright and Ting, 2006). There is also evidence of P. scaber inhabiting locations with elevations up to 2100 m a.s.l. (Beron, 2008). Originally, P. scaber woodlice occurred throughout Europe, but after spreading to other continents, the species is now distributed globally (Schmalfuss, 2003). Indirectly, such a wide geographic distribution suggests that P. scaber has high colonization potential and may adapt to a range of environmental conditions. The lungs of P. scaber represent the evolutionarily advanced type, characteristic of species well adapted to dry habitats, and they occur in two sexually dimorphic pairs of pleopods ( Fig. 2 A-D).
Our study explored developmental changes in the lung size in P. scaber exposed to a gradient of thermal and oxygen conditions. We hypothesized that the increased demand for oxygen combined with an oxygen supply shortage will favor an increased gas exchange area. To address this hypothesis, we reared individuals of P. scaber from egg to an adult in four combinations of temperature and oxygen levels. If our hypothesis was correct, then poor access to oxygen, especially in combination with increased metabolic demands at high temperatures, would result in increased lung size. To assess links between changes in lung size and the life history strategy, we additionally estimated growth curves for the woodlice exposed to our treatments. If developmental plasticity in lung size helps to meet the metabolic demand for oxygen, then we A structure of an exopodite with pleopodal lungs in Porcellio scaber (imaging with scanning electron microscopy). A: A right female exopodite of the first pair of pleopodsa general ventral (internal) view. B: Cells of the perispiracular area. C: Surface of the pleopodal lungs. Abbreviations: cl, connection of exopodite lamella with the base of pleopod; dist, distal part; lu, lung branches; mf, microfolds of the external surface of the perispiracular cells; pa, cuticulized perispiracular area; prox, proximal part; s, funnel-shape area leading into a spiracle. expected no dramatic differences in the growth patterns in those exposed to different oxygen treatments.

Parental animals
This work was a part of a long-term study on P. scaber previously described elsewhere (e.g., Horv� athov� a et al., 2015a). In autumn 2013 and spring 2014, we collected specimens from a population of P. scaber in an old backyard in Krak� ow, Poland (50 � 04'15.9"N 19 � 56'21.9"E) and transferred them to the Institute of Environmental Sciences (Jagiellonian University, Krak� ow). The field-collected animals served as the parents for a new generation, which was produced under controlled conditions and used in laboratory experiments (see below). The field-collected animals underwent at least a two-week acclimation to common laboratory conditions before being subjected to any experimental treatments. The animals were fed ad libitum with a mixture of dried leaves of black alder (Alnus glutinosa) and European ash (Fraxinus excelsior) that were collected in a nearby forest. The females and males were maintained separately in plastic boxes that contained moist sand and pieces of broken clay pots as shelter. The animal boxes were maintained in temperature-controlled cabinet (Pol-Eko, Poland), with a 12 L:12 D photoperiod and a temperature of 15 � C during the day and to 8 � C at night, which simulated typical autumn/spring conditions in the Krak� ow area (short-day conditions). Acclimation was the initial step in our procedure for stimulating synchronous offspring production among the parental generation (see below).

Experimental setup
The experiment was performed on a new generation of woodlice, which were obtained under controlled laboratory conditions. The experiment involved four developmental treatments created by combining two temperatures (15 and 22 � C) with two oxygen concentrations (10% and 22%) in two thermostatic cabinets (Pol-Eko, Poland) set to either 15 � C or 22 � C and four 110-liter plexiglass chambers (40 � 50 � 55 cm; YETI, Agencja Reklamy, Kryspin� ow, Poland), with two chambers placed in each cabinet. One of the two chambers per cabinet was set to normoxic conditions (22%) and the other to hypoxic conditions (10%). Oxygen conditions inside the chambers were controlled by a four-channel ROXY-4 gas regulator (Sable Systems International, USA). The regulator was connected to gas tanks with nitrogen and oxygen (Air Products, Poland), supplying the normoxic chambers with oxygen (normoxia) and the hypoxic chambers with nitrogen. Inside the chambers, the air circulation was maintained with fans, and the relative . The photos were used to measure the area of flat-laying exopodites, our proxy of lung size. Abbreviations: cl, connection of exopodite lamella with the base of the pleopod; dist, distal part; la, lung area; pa, cuticulized perispiracular area; prox, proximal part; s, funnel-shape area leading into a spiracle. humidity was controlled by four dew point generators DG-4 (Sable Systems International, USA), which were set to 75% relative humidity and connected to chambers via PP2 pumps (Sable Systems International, USA). The temperature and humidity inside the chambers were monitored with Hygrochron iButtons (Maxim/Dallas Semiconductor, USA). The relative humidity inside the boxes with the experimental animals reached 98%, which closely resembled conditions in the microhabitat occupied by the source population (Horv� athov� a et al., 2015a).
To stimulate synchronous mating and egg production, the parental animals were allocated to 28 boxes for mating (40 males and 50 females in each box; in total, 1120 males and 1400 females), and the boxes were placed in the Plexiglas chambers that maintained the specific experimental conditions. This ensured that the new generation of woodlice experienced experimental conditions from the very beginning of their life cycle. Following McQueen and Steel (1980), the stimulation of offspring production involved exposing the short-day acclimated animals to a photoperiod with a prolonged light phase (16 L:8 D). Initially, all four experimental treatments involved descendants of the parental animals collected in autumn 2013, but due to the failure of one of our thermostatic cabinets, we lost all animals from the warm treatment. Given an expected developmental lag between cold-and warm-reared woodlice, we decided to continue the cold treatment and run the warm treatment again later, using descendants of the parental animals collected immediately after winter (spring 2014). There is no reason to believe that the parental animals give birth to phenotypically different offspring in autumn vs spring, which would bias our comparisons of cold vs warm animals. All parental animals used in the experiment originated from the same population and represented the same generation/cohort. Importantly, the experiment was performed on a new generation of woodlice, which was produced under controlled laboratory conditions. Moreover, prior to the stimulation of mating in the laboratory, all parental animals experienced prolonged acclimation to the common laboratory environment and underwent a shift from short-day to long-day conditions, which should homogenize the phenotypic effects (if any) of prior environmental experience in the field-collected individuals.
The boxes with mating animals were checked weekly for the presence of gravid females. In isopods, egg laying and early juvenile development take place in a brood pouch called a marsupium, which develops on the ventral side of the female body during the parturial molt (a molt preceding reproduction) and is shed during the molt following reproduction. The offspring undergo twenty marsupial developmental stages (Milatovi� c et al., 2010;Wolff, 2009). After leaving the marsupium, they go through two manca stages and twelve juvenile stages (Tomescu and Craciun, 1987;Zimmer, 2002b). In our experiment, gravid females were individually transferred to 100-ml plastic boxes with a layer of wet sand, a piece of a clay pot for shelter and dry leaves for food. The boxes with gravid females and subsequently the boxes with offspring were maintained under the same experimental conditions (in the same thermal cabinets with oxygen regulated chambers). The boxes were checked weekly, and if free-living mancae were observed, the female was removed from the box. The boxes with the offspring were sprayed weekly with water and supplemented with our standard mixture of tree leaves (see Animals section). During the first nine weeks of postmarsupial life, the offspring were fed with only leaves. Two-week-old juveniles were additionally fed one dead conspecific adult to provide them gut microbiota, which are important in this species for juvenile growth (Horv� athov� a et al., 2015b). To provide enough protein, after the 4th week of life, the animals were fed one dried mealworm (Tenebrio molitor) once a week.
Individuals of the same sex from different clutches were pooled at the age of 20 weeks to form new single-sex groups, with 20 animals per group. Each group was placed in a separate box. When the majority of females reached a body mass within the body mass range of mature woodlice in the wild (minimum 21.68 mg and average weight 76.83 mg following Antoł and Czarnoleski (2018)), we added several males to each box to provide conditions for physiological maturation (Kight, 2008).

Dissections and measurements
To represent the respiratory capacity of P. scaber, we measured the areas of four flat-lying exopodites, and each woodlice was represented by the total area of exopodites (hereafter exopodite area), a proxy of lung size. According to Hoese (1982;Figs 5, 9) and Schmidt and W€ agele (2001; Fig. 4), the pseudotracheal system occupies a significant portion of the exopodites in P. scaber and, compared to Armadillo and Armadillidium spp. (Isopoda), the exopodites of P. scaber have the most two-dimensional (flat) structure. To validate our proxy of lung size, we measured the areas of the whole exopodites and a dark area of these exopodites, which corresponds to a visible net of pseudotracheae (personal observations), and then we correlated the two areas. Note that the validation involved a subsample of 90 exopodites (spread evenly between sexes and the two pairs of exopodites) for which the dark areas with pseudotracheae were sufficiently visible for the measurements.
We started measuring the exopodites and body mass from the 5th (22 � C) and 8th (15 � C) months of postmarsupial life, when the animals reached a similar body size; the difference in time reflected the temperature-dependent difference in development. The measurements were conducted every two months over a 17-month period. The sex ratio was female biased, with approximately only 30% males; therefore, we performed the measurements only on females, which were randomly sampled from each experimental treatment. Note that the measurements of pleopods required sparing animals for dissections; therefore, each animal was measured only once in a lifetime and was characterized by age, pleopod size and body mass upon being sampled. Additionally, at the end of the experiment, we sampled exopodites from individuals of both sexes (assumed to be fully mature animals). Each animal was weighed to the nearest 0.001 mg on a microbalance before dissection (Mettler Toledo XP 26, Mettler Toledo, Switzerland). During dissection, each individual was decapitated with a scalpel on a Petri dish. The body was submerged in 1x PBS (Avantor Performance Materials, Poland), and the exopodites were removed with forceps, placed in an Eppendorf-type tube and fixed in 10% buffered neutral formalin (Krakchemia, Poland). Upon collection, samples were washed in 1x PBS, placed in a 1x PBS drop on a microscope slide and covered with a cover slip. Photos of the exopodites were taken in a bright field under a light microscope (Eclipse 80i, Nikon, Japan) equipped with a digital camera (Axio Cam MRc5, Zeiss, Germany) and image acquisition software (ZEN, Zeiss) under 4x objective magnification (Fig. 2). Photos were used to measure the areas of the exopodites as well as their dark areas containing pseudotracheae. Additionally, to visualize the structure of the exopodites, one pair of male and one pair of female formalin-fixed exopodites were prepared for scanning electron microscopy (SEM). After hydration, the exopodites were fixed in 1% osmium tetroxide solution (Sigma-Aldrich Co., USA) in 1.6% potassium hexacyanoferrate (Chempur, Poland) for 2 h at 4 � C. Then, they were washed in distilled water, dehydrated in a series of ethanol concentrations and acetone, dried in a CO 2 critical point drier (LADD, USA), sputter-coated with gold and viewed and photographed with a scanning electron microscope using 20 kV (HITACHI S-4700, Japan) in the Scanning Microscopy Laboratory at the Institute of Geological Sciences Jagiellonian University.
To measure the dark areas of the exopodites with pseudotracheae (validation of our proxy of lung size), we manually outlined the edges of the dark areas of the exopodites with ImageJ image analysis software (NIH, Bethesda, USA). The measurements of the areas of whole exopodites were performed with a multistep semiautomatic segmentation algorithm based on a Canny (1986) edge detector and morphological operations. More precisely, the successive steps of the algorithm were as follows: (i) load the RGB image; (ii) transform the image to grayscale image; (iii) downscale the image; (iv) blur the image; (v) apply the Canny edge detector; (vi) apply morphological dilation; (vii) fill the holes; (viii) apply morphological opening; (ix) apply morphological erosion; and (x) upscale the segmentation. The algorithm was implemented in MATLAB (The MathWorks, Inc., USA), which returned an image with coarse segmentation that was then corrected by the researcher. After correction, the software generated a CSV file with the area measurements in μm 2 as the final output. For each individual, we calculated the summed area of the four exopodites (exopodite area), our proxy of lung size, which served as a measure of the capacity for obtaining oxygen from the air. We excluded individuals for which we did not have a complete set of lung measures. Eventually, we obtained complete data from 334 growing individuals sampled over the course of the experiment and from 120 mature individuals sampled at the end of the experiment. Data obtained from the former group was used to study the effects of experimental treatment on the ontogenetic changes in the size of pleopods and body mass (growth curves); data from the latter group were used to study the effects of experimental treatments on the pleopod size and body mass of adult individuals, considering sex differences in these two traits.

Statistical analysis
Statistical analyses were conducted in R software (R Core Team, 2018) with use of the effects (Fox and Weisberg, 2018) and ggplot2 (Wickham, 2016) packages, for mature isopods of both sexes and for the growing females sampled in the course of their development, respectively. To validate our proxy of lung size, we correlated the total of the dark areas of exopodites with visible pseudotracheae and the area of the whole exopodites. The analysis was performed separately for each pair of exopodites, while considering sex and developing woodlice vs adult woodlice. The data on our proxy of lung size (exopodite area) were analyzed with the general linear model with the exopodite area as a dependent variable, oxygen and temperature as fixed predictors, and body mass as a covariate. A similar model was applied for the mature animals, but it included sex as an additional fixed factor. We fitted the full models with all possible interactions and then applied a step-wise procedure that used the Akaike information criterion (AIC) to obtain the best model. If an interaction between a covariate and categorical predictor was significant, we tested the significance of the main effects at the minimal, mean, and maximal values of the covariate. If the interaction of two (or more) categorical factors was significant, we analyzed differences between groups with a general linear hypothesis test using the multcomp package (Hothorn et al., 2008). Body mass and exopodite area were log transformed prior to the analysis.
To analyze links between growth and respiratory capacity, we used the non-linear least square procedure (nls) to fit von Bertalanffy growth curves (von Bertalanffy, 1957) to our body mass, M t , data and the corresponding age, t, of animals. This analysis involved only data for females that were sampled for pleopod dissections over the course of development (note that body mass and age for these animals were recorded only once in their lifetime). Although the mechanism proposed by Bertalanffy to explain the origin of indeterminate growth does not follow evolutionary principles, Bertalanffy's mathematical formula can be used as a tool to describe the indeterminate growth pattern phenomenologically (Czarnoleski and Kozłowski, 1998). We used a three-parameter formula: M t ¼ M A ð1 À e À kðtÀ t0Þ Þ 3 , where M A is an asymptotic size, K is a growth rate coefficient, and t 0 is the hypothetical age in which an animal would achieve a body size equal to 0.
Generally, our comparison of growth curves between experimental treatments involved two steps. The first step aimed at assessing overall differences in the shape of growth curves, and it involved the nls procedure for curve fitting, combined with random shuffling of data between experimental groups. This allowed us to calculate a Diff statistic that measured a degree of shape differences between growth curves (see below our detailed description of this method). The second step tested directly which of the three Bertalanffy's parameters was a main driver of the growth curve differences detected in the first step. The second-step analysis used the nls procedure for fitting Bertalanffy's model to data, additionally considering effects of our two treatments (temperature and oxygen) on values of Bertalanffy's parameters. Technically, we converted our two categorical predictors (temperature and oxygen treatment groups) to two numeric dummy variables, each taking values of either 0 or 1 (thermal treatment: value 0 assigned to 22 � C and value 1 assigned to 15 � C; oxygen treatment: value 0 assigned to 22% and value 1 assigned to 10%). Fitting Bertalanffy's growth formula with nls, the model considered an independent effect of each dummy variable on each growth parameter. For simplicity, the final model derived by our second-step analysis did not consider interactive effects of dummy variables on growth parameters because all interactions appeared to be highly nonsignificant.
Here, we explain in more detail methods applied in the first step of our growth curve analysis. This analysis used an approximate randomization method for growth curve comparisons described by Czarnoleski et al. (2005) and Weinberg and Helser (1996). To compare two curves, we fitted Bertalanffy's curves to two datasets (1 and 2), separately for each set, obtaining two curves with respective residual sums of squares, SS 1 and SS 2 . After pooling the data from the two datasets, we fitted Bertalanffy's curve to the pooled data (1 þ 2), obtaining a new residual sum of squares SS 1þ2 . This method considers whether the two curves being compared are identical; the total residual sum of squares does not change after fitting a single curve to the pooled data (SS 1 þ SS 2 ¼ SS 12 ), but it increases if the curves have different shapes (SS 1 þ SS 2 < SS 12 ). Therefore, the statistic Diff ¼ SS 12 -(SS 1 þ SS 2 ) measures the shape difference between two curves. If the curves are identical, then Diff is equal to 0, but the more the curves differ, the higher the value of Diff become. To assess the statistical significance of the shape difference in the curves, we compared the Diff statistics with the distribution of their randomly generated values. The distribution was produced via 10,000 random shufflings of data. Each time, data from the two groups under comparison were pooled and randomly reassigned to the groups, keeping the number of observations per group the same as in the original groups. After each shuffling, we calculated Diff* according to the same procedure used to calculate Diff. Note that the distribution shows how frequently a given difference between curves (Diff) arises entirely from random processes, so when there is only one general population, two groups of data on age and body size are drawn randomly. Statistical significance represented the probability (p* value) that the randomly generated Diff* statistic was as large as or larger than the observed Diff statistic. The randomization procedure was performed in R.
In our 2x2 experiment, thermal conditions (warm and cold) were crossed with oxygen conditions (normoxic and hypoxic). Consequently, to compare the growth curves between the warm and cold conditions, we had to simultaneously control the effects of oxygen conditions, and vice versa, i.e., the growth comparison between hypoxia and normoxia involved controlling the effects of thermal conditions. Therefore, to compare the curves between the warm and cold conditions, we calculated the Diff statistic separately for the data obtained under normoxic and hypoxic conditions. After summing these two Diffs, we obtained a Diff Temp statistic, which measured the shape difference of the temperature curves while simultaneously integrating information from the normoxic and hypoxic conditions. Note that when the thermal conditions do not affect the growth curve, i.e., the curve is consistent between hypoxic and normoxic conditions, then Diff calculated for normoxic conditions equals 0 and Diff calculated for hypoxic conditions also equals 0, which leads to Diff Temp ¼ 0. The more the growth curves differ between temperatures, and the more consistent this effect is in normoxic and hypoxic conditions, the higher the value of Diff Temp . If the temperature and oxygen have interactive effects on the growth curves, e.g., when Diff ¼ 0 for normoxic conditions (no thermal effects on growth curves in normoxic conditions) and Diff >0 for hypoxic conditions (thermal effects on growth curves in hypoxic conditions), Diff Temp has an intermediate value. To compare the curves between normoxic and hypoxic conditions, we calculated the Diff statistic separately with the data obtained from warm and cold conditions, and summing these values, we obtained Diff Oxy . This statistic measured the shape differences in the oxygen level curves, integrating information obtained under either warm or cold conditions. Finally, with reference to our randomly generated distributions (see previous paragraph), we evaluated the statistical significance of Diff Temp and Diff Oxy .

Exopodite area
When we analyzed the data on female woodlice in the course of their growth and development (Table 1; Fig. 3), we found a significant relationship between body mass and exopodite size (p < 0.001, log-log regression slope ¼ 0.74). There was also a significant body mass and oxygen interaction (p ¼ 0.01), and the regression slope was increased in low oxygen by 0.04. The temperature and log body mass interaction was close to significance (p ¼ 0.07), so we decided not to ignore it in our interpretation of model results (the regression slope was decreased in high temperature by 0.03). To improve the interpretability of our model with two interactions, we examined the results of the same model after centering the covariate (log body mass) at its two extreme and mean values (Fig. 3). Small animals had smaller exopodites in hypoxic conditions than in normoxic conditions (t ¼ À 2.04, p ¼ 0.04); the effect of temperature was not significant (t ¼ 1.22, p ¼ 0.22). In average-sized animals, we did not find effects of oxygen level on exopodite area (t ¼ 0.92, p ¼ 0.36), but increased temperature resulted in decreased exopodite area (t ¼ 1.94, p ¼ 0.05). Large animals had larger exopodites in hypoxic conditions than in normoxic conditions (t ¼ 2.66, p ¼ 0.01) and relatively smaller exopodites at increased temperatures (t ¼ 2.63, p ¼ 0.01).
In mature woodlice (Table 2; Fig. 4B), we found a significant increase in exopodite size associated with body mass (p < 0.001). Interestingly, there was a significant sex and body mass interaction (p ¼ 0.05); in females, this relationship was more prominent (the scaling exponent equaled 0.64) than that in males (the scaling exponent equaled 0.57). Regardless of sex, increased temperatures reduced the mean area of exopodites (p ¼ 0.004, Fig. 4A). There was a significant oxygen and sex interaction (p ¼ 0.002), with exopodites being larger under hypoxic conditions than normoxic conditions, but only in males (t ¼ 4.08, p < 0.001; Fig. 4B). In females, the size of the exopodites was not significantly different between oxygen groups (t ¼ 1.36, p ¼ 0.41, Fig. 4B). As the sex and covariate interaction was significant, we further analyzed the differences between the sexes, centering the covariate (log body mass) at three different values. Males had consistently larger exopodites than females at the minimal (t ¼ 8.76, p < 0.001), mean (t ¼ 16.81, p < 0.001) and maximal (t ¼ 7.59, p < 0.001) log body mass, but this sexual dimorphism was less visible among large woodlice.

Growth curves
In the first step of the analysis of Bertalanffy's growth curves, our randomization tests showed that oxygen conditions did not affect growth (p* ¼ 0.87, Fig. S1), but growth curves differed between our two thermal conditions (p*<0.0001, Fig. S1). In the second step of growth curve analysis (Fig. 5, Table 3) we explored which of Bertalanffy's parameters were main drivers of the detected thermal plasticity of growth curves. Our results show that cold-reared woodlice under either hypoxic or normoxic conditions were consistently characterized by a smaller growth rate coefficient (K) and larger asymptotic body mass (M A ), but note that only the result for K was close to statistical significance (t ¼ 1.79, p ¼ 0.07). In accordance with the results of our first-step analysis, we did not find significant effects of oxygen conditions during development on either of the three Bertalanffy's parameters (Table 3). Overall, we conclude that our two-step analysis of Bertalanffy's growth curves indicated that thermal but not oxygen conditions during development significantly affected the growth pattern of woodlice females, but none of the three Bertalanffy's parameters alone drove the detected thermal changes in growth pattern.

Discussion
Our study of P. scaber demonstrated that larger woodlice developed larger pleopodite-exopodites with lungs, which was evident among both growing subadult woodlice females and male and female adults with different body sizes. The nature of this mass-scaling suggests that the gas-exchange surface of pleopodal lungs might follow a negative allometry with body mass (scaling exponents smaller than 1). In accordance, studies of the aquatic species of isopods Asellus aquaticus and Saduria entomon (Babula, 1979;Babula et al., 1978;Bielawski and Babula, 1980) reported shallow mass-scaling of the area of respiratory surface and body mass (mass-scaling exponent 0.74 for A. aquaticus and 0.688 for S. entomon). Certainly, these conclusions should be treated with caution because these studies, including ours, did not measure directly the area of gas-exchange surface but instead used proxies that allow more for qualitative than quantitative conclusions. The positive relationship between body size and the size of the gas exchange organs is certainly a prevailing phenomenon in nature, and it was previously documented in many distantly related organisms such as turtles (intra and interspecifically) (Hochscheid et al., 2007), bats (interspecifically) (Maina et al., 1991), fish (intraspecifically), mammals (interspecifically), lizards (intraspecifically) and birds (interspecifically) (Hughes, 1984). Consistent with the principle of symmorphosis, organisms are expected to coordinate development to match the physiological capacities of different tissues and organs (Taylor and Weibel, 1981), so a change in the size of an isopod's lungs should be paralleled by adequate changes in other elements of the oxygen transportation system. However, in contrast with the expectations of symmorphy, mammals have been shown to develop excessive gas exchange capacity in the lungs in relation to the capacities of other elements of the oxygen transportation system (Weibel et al., 1991). Certainly, a better understanding of the capacity of the gas exchange system in isopods would require information on cuticular gas exchange (Edney and Spencer, 1955), the amount and binding efficiency of hemocyanin, the volume and flow rates of hemolymph, and the diffusion of oxygen through the tissue and cells.
The results of our experiment demonstrated that the size of the gas exchange organs in woodlice developmentally responded to environmental conditions, which, to our knowledge, might be the first report of this kind of phenotypic plasticity in isopods. The effect of hypoxia on tracheal architecture has been demonstrated in Blatella germanica Table 1 Results of the general linear model comparing the exopodite areas (proxy of lung size) in female Porcellio scaber during their postmarsupial development at two oxygen levels (10% and 22%) and two temperatures (15 and 22 � C). cockroaches (Vandenbrooks et al., 2012) and Drosophila melanogaster flies (Centanin et al., 2010;Harrison et al., 2018b). In line with the general idea that the gas exchange surface of isopod lungs is matched to the metabolic demand and oxygen supply, we found that the size of P. scaber lungs was sensitive to thermal and oxygen conditions during development. Nevertheless, the pattern of these responses suggests a compromise between the rates of oxygen uptake and water loss. It is especially notable that the size of lungs increased differently during ontogenetic growth of P. scaber females, depending on the developmental temperature. Small woodlice had similar sized lungs irrespective of thermal conditions, but among large woodlice, warm woodlice had smaller lungs than cold woodlice (Fig. 3). Consistently, when we analyzed our dataset for mature males and females, we found that lungs were smaller in adult woodlice originating from our warm treatment than from the cold treatment (Fig. 4). The development of relatively smaller lungs under warm conditions does not follow the expectation that a larger gas exchange surface helps ectotherms in warm conditions meet their increased demand for oxygen. However, the advantage of large lungs in warm environments can be reduced by the positive effect of temperature on oxygen diffusion rates combined with the increased rates of water evaporation via the gas exchange surface. This scenario agrees with the macroevolutionary trends of covering the lungs and limiting possible surface evaporation, which is associated with adaptation of species to more xeric environments (Hoese, 1982;Hsia et al., 2013;Schmidt and W€ agele, 2001). Consistent with the idea that an increased gas exchange surface facilitates oxygen uptake in hypoxic environments, we found that among mature adults, males reared under hypoxic conditions had larger lungs than those reared under normoxic conditions. The same response pattern also occurred in females, but only among growing individuals with a relatively larger size: as females grew in size, their lung areas  increased faster under hypoxic than under normoxic conditions. Intriguingly, the opposite effect of oxygen conditions on lung size characterized growing females if we considered smaller individuals (but notably, younger and thus smaller females were exposed to the oxygen treatment over a shorter time interval than the larger and older females). Note here that in addition to the lungs, isopods also use the outer surfaces of other body parts for gas exchange with the environment (Edney and Spencer, 1955). Therefore, it is worth considering that large woodlice with their low surface-to-volume body ratio should benefit more from large lungs than small woodlice characterized by relatively large body surface. This effect might explain why we observed large but not small woodlice responding to hypoxia via the development of larger exopodite lungs. Nevertheless, this mechanism cannot explain why males and females responded differently to oxygen conditions in their lung size; among adults, the lungs of males became larger under hypoxic conditions than those of males under normoxic conditions, but the lungs of females remained the same size under different oxygen conditions. Notably, P. scaber males developed relatively larger lungs than females, so it appears that the sex with the relatively larger exopodites (males) responded to hypoxia by enlarging the lungs, while the sex with the relatively smaller exopodites did not respond to oxygen conditions. We believe that this intersexual difference in exopodite/lung size should be interpreted at least partially with consideration of a reproductive function of exopodites (see also Fig. 2) -the second pair of male exopodites takes part not only in gas exchange but also in sperm transfer. However, this function still does not help to explain why males and females react differently to hypoxia in terms of their pleopod exopodite size. The sexual dimorphism of lung size and the sex differences in oxygen dependence of lung size suggest that males might have higher metabolic demand for oxygen or be less water constrained than females. Spiders have evolved intersexual metabolic differences, though which sex has an increased metabolic rate depends on the species (Kotiaho, 1998;Tanaka and Itô, 1982;Watson and Lighton, 1994). In contrast, mice and mosquitofish show little or no intersexual difference in metabolic rates (Cech et al., 1985;Schulz et al., 2002). Given the reproductive biology of isopods (Kight, 2008), increased metabolic demands might be expected in females rather than in males of P. scaber because females carry offspring in their marsupial pouches and provide offspring with resources and oxygen. Nevertheless, we consider that for females, the important ecologically limiting factor might be the water supply because they supply their offspring with aqueous conditions for a considerable part of their marsupial development (Horv� athov� a et al., 2017). In fact, Horv� athov� a et al. (2017) demonstrated that females of P. scaber increase the duration of the aqueous phase inside the marsupium in hypoxic environmental conditions, thereby delaying entry of offspring into gaseous surroundings. If so, we hypothesize that female terrestrial isopods might be under stronger selective pressure than males to compromise their capacity to obtain oxygen for the benefit of water conservation. Available information on thermal sensitivity of growth patterns shows that thermal conditions during development exert more pronounced effects on the growth of aquatic than terrestrial species (Horne et al., 2015). Our data on growth of woodlice under thermal and oxygen treatments indicated that the oxygen conditions in our experiment did not significantly modify the growth pattern, which is consistent with our previous study on the growth of woodlice during their first 13 weeks of postmarsupial development, where we did not find an effect of oxygen (Horv� athov� a et al., 2015a). However, we found that cold-reared woodlice grew slower and obtained larger asymptotic body mass than warm-reared woodlice, which agrees with the predictions of the so-called temperature-size rule (TSR), a thermal plasticity pattern that has evolved in many ectotherms (Atkinson, 1994). Thus far, the thermal dependence of growth patterns has rarely been studied in isopods. On a geographic scale, the woodlouse Porcellio laevis follows the TSR and exhibits an increased body size in Chilean populations that are located further from the equator (Lardies and Bozinovic, 2006). Moreover, the TSR was demonstrated experimentally in the freshwater isopod Asellus aquaticus under hypoxic conditions (Hoefnagel and Verberk, 2015). Our finding that oxygen conditions did not modify the growth curves of P. scaber does not follow the expectations that the TSR is influenced by changes in oxygen availability combined with oxygen demand Hoefnagel and Verberk, 2015;Walczy� nska et al., 2015). Certainly, the nonsignificant effects of oxygen availability on the growth curves in our experiment might be caused by relatively good access to oxygen even in our hypoxic treatment (10%), especially for litter dwelling animals such as P. scaber. It is also possible that the developmental plasticity of the lungs revealed in our experiment indicates compensatory changes that reduce the effect of oxygen availability on growth capacity. In nature, isopods are challenged with desiccation risk, which was minimized in our experiment, as we kept the animals under conditions with the highest possible humidity. The high humidity might have allowed them to develop lungs large enough to maintain an adequate oxygen supply (but note our sex-specific responses in pleopod size to oxygen conditions).
Our previous study (Antoł et al., 2019) demonstrated that the woodlouse P. scaber responded behaviorally to acute hypoxia by choosing cooler sites, which possibly helped survival by decreasing the Fig. 5. Von Bertalanffy's growth curves for female Porcellio scaber reared from eggs to adults under two oxygen levels (10% and 22%) combined with two temperatures (15 and 22 � C). The curves for the actual range of data used in the curve fitting are displayed, the parameters were calculated with use of nls procedure; Table 3 provides detailed results of the statistical model.

Table 3
Effect of temperature and oxygen on von Bertalanffy growth curve M t ¼ M A ð1 À e À kðtÀ t0Þ Þ 3 parameters for female Porcellio scaber that developed under two oxygen levels (10% and 22%) combined with two temperatures (15 and 22 � C). M A -asymptotic body mass, Kgrowth rate coefficient, t 0 -a hypothetical age at which zero size would be achieved. Parameters were estimated with use of non-linear least square (nls) model that considered effects of temperature and oxygen in a form of two dummy variables affecting each Bertalanffy's parameter (thermal treatment: 0 assigned to 22 � C and 1 assigned to 15 � C; oxygen treatment: 0 assigned to 22% and 1 assigned to 10%). Graphical visualization of curves predicted by the model is shown in Fig. 5 Here, we demonstrated that during development in different oxygen and thermal conditions, P. scaber underwent changes in the size of the gas exchange organs, which might be explained by a balance between the capacity to obtain oxygen and conserve water. Oxygen availability, thermal conditions and water conditions might be key environmental factors involved in the early evolution of terrestrial isopods (Horv� athov� a et al., 2017). Approximately 300 Mya, after isopods invaded land, oxygen availability drastically dropped (by 50%; Berner et al., 2007), and this shift coincided with the evolutionary diversification of land isopods. During that time, this group of crustaceans evolved numerous adaptations (such as lungs), which must have played a role in the colonization of various terrestrial habitats. Our results suggest that this evolution preserved genotypes that have the capacity to respond in a plastic way to developmental conditions, which would help explain why terrestrial isopods thrive among a wide range of environmental gradients.

Declaration of competing interest
The authors declare no conflicting interests.