Physiological, morphological and ecological traits drive desiccation resistance in north temperate dung beetles

Background Increasing temperatures and changes in precipitation patterns threaten the existence of many organisms. It is therefore informative to identify the functional traits that underlie differences in desiccation resistance to understand the response of different species to changes in water availability resulting from climate change. We used adult dung beetles as model species due to their importance to ecosystem services. We investigated: (i) the effect of physiological (water loss rate, water loss tolerance, body water content), morphological (body mass) and ecological (nesting behaviour) traits on desiccation resistance; (ii) the role of phylogenetic relatedness in the above associations; and, (iii) whether relatively large or small individuals within a species have similar desiccation resistance and whether these responses are consistent across species. Results Desiccation resistance decreased with increasing water loss rate and increased with increasing water loss tolerance (i.e. proportion of initial water content lost at the time of death). A lack of consistent correlation between these traits due to phylogenetic relatedness suggests that the relationship is not determined by a shared evolutionary history. The advantage of a large body size in favouring desiccation resistance depended on the nesting behaviour of the dung beetles. In rollers (one species), large body sizes increased desiccation resistance, while in tunnelers and dwellers, desiccation resistance seemed not to be dependent on body mass. The phylogenetic correlation between desiccation resistance and nesting strategies was significant. Within each species, large individuals showed greater resistance to desiccation, and these responses were consistent across species. Conclusions Resistance to desiccation was explained mainly by the dung beetles’ ability to reduce water loss rate (avoidance) and to tolerate water loss (tolerance). A reduction in water availability may impose a selection pressure on body size that varies based on nesting strategies, even though these responses may be phylogenetically constrained. Changes in water availability are more likely to affect dweller species, and hence the ecosystem services they provide. Supplementary Information The online version contains supplementary material available at 10.1186/s40850-021-00089-3.


Box S4 The phylogenetic analysis
For the phylogenetic analysis, 35 COI-5P sequences (ingroup and outgroup taxa) were obtained from the BOLD and GenBank databases (see Table S3 for the list of the accession numbers), and aligned using the muscle method by Mega v10 (UPGMB clustering method, Kumar et al., 2018). The dataset was analysed by Maximum Parsimony (MP) and Maximum Likelihood (ML) approaches. In both analyses the nonparametric bootstrapping was performed to assess statistical significance for support of internal nodes (Felsenstein, 1985).
In order to evaluate the best fit model of sequence evolution for the ML analysis, the fit statistics among 24 different nucleotide substitution models were calculated ( Table S5). The evolutionary divergence between sequences was also estimated, and the number of base substitutions per site with the standard error are given ( Table S6). The following options were used for the MP heuristic search: search method = TBR; initial trees = 10; max No. trees retained = 100; substitutions type = nucleotide level; bootstrap n reps = 1000. The ML analysis options were set according to the best fit value (Table S5) with the substitution model = GTR and rates = G+I; the other options were set as substitutions type = nucleotide level; heuristic method = NNI, with default initial tree; bootstrap n reps = 1000; branch swap filter = very strong. All the trees generated by the analyses showed the same phylogenetic relationships among the taxa. Here the ML tree is presented ( Figure S7).
For the Phylogenetic Signal (PhylSig) analysis, a reduced matrix was built, in which only a terminal node for each ingroup species was included (    The ML tree with the highest log likelihood (-3191.7887) is shown. The ML tree was obtained by the GTR (General Time Reversible) model. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.9309)). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 52.8325% sites). The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The tree was rooted by the outgroup method, using the Trypocopis species as outgroup.

Figure S8
The bootstrap tree from the ML tree with the highest log likelihood (-2731.88); on the tree the numbers of the internal nodes were added (in red). The evolutionary history was inferred by using the Maximum Likelihood method, GTR model. The percentage of trees in which the associated taxa clustered together is shown next to the branches (in black). Here, all the branches are showed as fully resolved. A discrete Gamma distribution was used to model evolutionary rate differences among sites, with 5 categories (+G, parameter = 1.62). The rate variation model allowed for some sites to be evolutionarily invariable ([+I], 55% sites).

Box S9 The Phylogenetic signal
The software PhyloCom v4.2 (Webb et al., 2008) was used to define the phylogenetic signal of the traits. A matrix of the mean values of each trait for each species and the reduced ML tree (Fig A4) were used to correlate the eight species and and six trais ( In the traits file the desiccation resistance, dry body mass, water loss rate, water loss tolerance, and initial water content were coded as continuous data, while the functional group was coded as ordered multistate. The AOT module was used to conduct the test of phylogenetic signal (= the tendency for close relatives to resemble each other) and traits correlations.
The variance of the standardized contrasts (VarContr, in Supplemental Material, Table S1) values were used as a measure of the phylogenetic signal, being the smaller values an index of stronger dependence of the trait, as is here the case for WRL and DR. If related species are similar to each other, a small variance of contrast values will results, since the magnitude of independent contrasts will be similar across the tree. In this framework, the VarContr values for BM and iWC vouch for a very great variance for both traits.
The evolutionary correlations between traits were calculated by using independent contrasts defining in turn each trait as the independent variable (traitX) and thence all the others as the dependent variable (traitY). The significance testing of the analysis was conducted by randomization (N = 9999 by default). The traits were analysed in pairs, and the significance of the resulting correlations (Table S2)