Challenging ecogeographical rules: Phenotypic variation in the Mountain Treeshrew (Tupaia montana) along tropical elevational gradients.

Bergmann's and Allen's rules were defined to describe macroecological patterns across latitudinal gradients. Bergmann observed a positive association between body size and latitude for endothermic species while Allen described shorter appendages as latitude increases. Almost two centuries later, there is still ongoing discussion about these patterns. Temperature, the common variable in these two rules, varies predictably across both latitude and elevation. Although these rules have been assessed extensively in mammals across latitude, particularly in regions with strong seasonality, studies on tropical montane mammals are scarce. We here test for these patterns and assess the variation of several other locomotory, diet-associated, body condition, and thermoregulatory traits across elevation in the Mountain Treeshrew (Tupaia montana) on tropical mountains in Borneo. Based on morphological measurements from both the field and scientific collections, we found a complex pattern: Bergmann's rule was not supported in our tropical mountain system, since skull length, body size, and weight decreased from the lowest elevations (<1000 m) to middle elevations (2000-2500 m), and then increased from middle elevations to highest elevations. Allen's rule was supported for relative tail length, which decreased with elevation, but not for ear and hindfoot length, with the former remaining constant and the latter increasing with elevation. This evidence together with changes in presumed diet-related traits (rostrum length, zygomatic breadth and upper tooth row length) along elevation suggest that selective pressures other than temperature, are playing a more important role shaping the morphological variation across the distribution of the Mountain Treeshrew. Diet, food acquisition, predation pressure, and/or intra- and inter-specific competition, are some of the potential factors driving the phenotypic variation of this study system. The lack of variation in body condition might suggest local adaptation of this species across its elevational range, perhaps due to generalist foraging strategies. Finally, a highly significant temporal effect was detected in several traits but not in others, representing the first phenotypic variation temporal trends described on treeshrews.

Introduction

Bergmann’s [1] and Allen’s [2] ecogeographical rules describe geographical variation of morphological traits primarily in endotherms. Bergmann’s Rule describes a pattern of individuals being larger in colder habitats. Bergmann [1] explained this trend through the heat conservation hypothesis: in colder climates (higher latitudes) larger animals, with a lower surface-area-to-volume ratio, increase their ability to conserve heat. Allen’s Rule [2] describes a pattern of individuals having shorter appendages in colder habitats (higher latitudes). The reduction of surface area to volume ratio through a decrease in appendage length increases thermal efficiency.

Although these rules were initially developed based on patterns across latitudinal gradients, they were later extensively tested across elevational gradients (e.g., [3-7]). Elevational gradients share certain properties with latitudinal gradients, including the key change in temperature [4]. Common characteristics of elevational gradients (e.g. decrease in temperature and oxygen levels as altitude increases) are consistent across large geographic scales, enabling the study of the influence of environmental variables in natural populations across numerous replicated systems. In contrast to latitudinal gradients, elevational gradients exhibit drastic environmental transitions at spatial scales that are small relative to the dispersal capability of many species, making these optimal to study whether adaptive divergence is possible when confronted with gene flow [4].

Many studies have tested these rules in mammals both within and between species, and in different habitats (e.g. [5, 8–13]). Some mammals have patterns consistent with these predictions (e.g. [3, 14–18]) while others do not (e.g. [3, 7, 15, 17, 19–24]). However, there are multiple biases in the literature. There is a tendency to test these rules in species that are known to have a high/moderate degree of morphological differentiation across their geographical distribution and ignore those that do not [10]. There is also a strong bias towards studies in habitats with a strong temperature seasonality, especially the Holarctic (e.g. [5, 7–9, 18, 24–26]). This bias is unfortunate because seasonality by itself likely affects the evolution of body size [27-30].

Despite the valuable insights obtained from elevational approaches, few studies have tested these rules in mammals along elevational gradients in tropical areas with little seasonality [31, 32]. Indeed, the few studies from tropical regions seem to refute Bergmann’s rule and provide mixed support for Allen’s rule. A negative effect of latitude on body size was found both in the Common Treeshrew in tropical Asia (Tupaia glis; [23]), and in peccaries in the American tropics (Tayassuidae; [33]) indicating the inverse of Bergmann’s rule. Bergmann’s rule was not supported either at the interspecific nor at the intraspecific level along elevational gradients in New Guinea for rats (Rattus sp.; [31]) or forest passerines [34]. At an interspecific level, Bergmann’s rule was not supported in Bornean birds, which decrease in mass with elevation, while Allen’s rule had mixed support, with high-elevation communities exhibiting narrower, shorter and relatively smaller bills but relatively longer tarsi [35]. Neotropical Torrent ducks (Merganetta armata) followed Bergmann’s rule across latitude, but not along elevation, while Allen’s rule was neither supported across latitude nor elevation, showing an opposite trend for wing length [36].

Vegetation zonation can be expected to have an important impact on habitat structure and resource availability, both of which might pose strong selection pressures on local populations [37]. Populations might adapt to these different habitats through changes in locomotion, diet-associated traits, and/or insulation (changes in fat levels and/or fur density/length to alter heat loss; [6, 31, 35]). If populations are well adapted to these habitats, limited variation in body condition should be expected along elevational gradients [4]. Nevertheless, local adaptation must be constrained to a certain degree since species usually have restricted elevational distributions [4].

A likely explanation for the few intraspecific studies of morphological variation across elevational gradients in tropical habitats without a strong temperature seasonality is the lack of species with an appropriate distribution to test this, since species in the tropics tend to have narrow elevational ranges [38]. The Mountain Treeshrew (Tupaia montana) provides a rare opportunity to test Bergmann’s and Allen’s rules along wide tropical elevational gradients. These are one of the few Bornean mammals distributed along a wide elevational gradient (800–3400 m) on some of the tallest tropical mountains in tropical east Asia, including Mt. Kinabalu [39].

Here we assess phenotypic variation across an elevational gradient (in which elevation is treated as a surrogate of temperature) in the tropics in the Mountain Treeshrew, an endemic Bornean species. We combine field and museum data to test Bergmann’s rule (increasing size with elevation using weight, head-and-body length, and skull length); Allen’s rule (relatively smaller appendages with increasing elevation using ear, hindfoot, and tail length as a proportion of head-and-body length); changes in body condition as a surrogate of lipid content which is frequently assumed to be directly positively linked to fitness; trophic variation using rostrum length, upper toothrow length, and zygomatic breadth as a proportion of skull length; changes in locomotor behavior using hindfoot, and tail length as a proportion of head-and-body length; and insulation using guard hair length.

Materials and methods

Model species

The Mountain Treeshrew (Tupaia montana) is a terrestrial/ scansorial member of the Order Scandentia that has an omnivorous diet [40, 41]. It is an endemic Bornean species and montane specialist, distributed above 800 m in the highlands of northern Borneo. It has been recorded in all four vegetation zones described on Mt. Kinabalu [42, 43], but it is rare in the lowland forest (< 1,200 m) and becomes much more common through the lower montane (1,200–2,000 m) and upper montane forest (around 2,000–2,800 m) where it reaches highest densities, and in the subalpine forest (2,800–3,400 m) where its population density possibly decreases again [39].

Data collection

Mountain Treeshrews were sampled in the field along elevational gradients on Mt. Kinabalu (4,095 m; 6.07° N 116.56° E), Mt. Tambuyukon (2,579 m; 6.20° N 116.66° E), Mt. Alab, Crocker Range (2,050 m; 5.83° N, 116.34° E), and Mt. Trus Madi (2,643 m; 5.55° N, 116.52° E) (, [39]). We used Tomahawk and Sherman live-traps baited with banana, dry fish, shrimp paste with rice flour, coconut, and/ or palm nuts. Trapping was performed following ethical standards according to the guidelines of the American Society of Mammalogists [44]. Animal care and use committees approved the protocols (Smithsonian Institution, National Museum of Natural History, proposal number 2012–04 and Estación Biológica de Doñana proposal number CGL2010-21424). Field research was approved by Sabah Parks (Refs: TS/PTD/5/4 Jld. 45 (33), TS/PTD/5/4 Jld. 47 (25), TK/PP:5/78 Jld 5), the Economic Planning Unit (Ref: 100-24/1/299), Sabah Wildlife Department (Ref: JHL: (HQ):100-42/ 1 JLD.27), Sabah Forestry Department [Refs: JPHTN/TP(FSP) 100-14/18/2/KLT.31(31) & JPHTN/TP(FSP) 100-14/18/2/KLT.31(17)] and the Sabah Biodiversity Council/Centre (Ref:TK/PP:8/8Jld.2). Two hundred thirty-three Mountain Treeshrews were captured from 836 to 3,382 m during three field seasons (2012, 2013 and 2016). We sampled a total of 81 unique individuals (ear clipping allowed identification of recaptured individuals) in 1737 trap-nights on Mt. Kinabalu (0.047 individuals per trap-night), 103 in 4369 trap-nights on Mt. Tambuyukon (0.023 individuals per trap-night), 9 in 960 trap-nights on Mt. Alab (0.009 individuals per trap-night), and 40 in 1600 trap-nights on Mt. Trus Madi (0.025 individuals per trap-night) (, individual BOR numbers as in [39]). These individuals constituted the Estación Biológica de Doñana collection (BOR and EBD-CSIC numbers, for released and collected specimens). We weighed animals with Pesola® 300 g with ± 0.3% precision scale. We measured body length as the distance from the tip of the snout to the base of the tail, with the ruler held along the dorsum, to the nearest 1 mm (head-and-body length, HBL); tail length (TL), from posterior margin of the anus to tail tip, excluding hairs; hindfoot length (HFL), from back of the heel to tip of the longest digit, including the claw; and ear length (EL), from notch at the base of the ear to distal edge. These measures were taken immediately after specimens were collected (or from live specimens in analyses). Sex was determined from external sexual organs, and in the case of collected animals, also through dissection and discrimination between ovaries, testes and immature gonads during specimen preparation. After specimen preparation, the following craniodental measurements were collected on the cleaned and dried skulls with a Fowler High Precision electronic digital caliper to the nearest 0.01 mm by A. H. and as described in [23]: condyloincisive length (CIL), as the greatest distance between anterior‐most surface of I1 and caudal surface of occipital condyle; zygomatic breadth (ZB), as the greatest distance between lateral surfaces of zygomatic arch; upper toothrow length (UTL), as the greatest distance between anterior‐most surface of I1 and posterior‐most surface of M3; and rostrum length (RL), as the greatest distance between anterior most surface of premaxilla and anterior most surface of lacrimal foramen (Fig 1 in ). We measured the right side for all bilateral measurements.

Physical map of study area.

Map of northern Borneo showing mountains and potential distribution of Mountain Treeshrews based on its elevation range (areas above 800 m). Areas between 800 and 1600 m are highlighted in light gray, those between 1600 and 2400 m in dark gray and those over 2400 m in black.

To complement our field data, specimen-label-associated measurements were incorporated from the collections of the Field Museum of Natural History (N = 33, from 1937 and 1965), Harvard Museum of Comparative Zoology (N = 72, collected in 1937), and Sabah Museum (N = 88, from 1971 to 2007), while skull and specimen-label-associated measurements were taken from United States National Museum of Natural History (N = 72, collected in 1951, 1953 and 1961) and Natural History Museum of London (N = 25, collected in 1952 and 1961). An additional mountain, Maligan Range (1953 m; 4.84° N, 115.76° E) was added to this study based only on museum data (). In total, data from 523 individuals collected between 1937 and 2016 were evaluated ().

We aged (adult/ non-adult) individuals from museum specimens based on gonad examination in necropsies of wet specimens or on the permanent dentition eruption and lack of deciduous teeth [45]. Individuals that were released in the field were aged based on HBL, body weight and the presence of external sexual organs.

We generated five datasets. The first dataset (D1) included weight measurements of museum specimen adults aged based on gonad examination or on the confirmation of their permanent dentition (N = 84). The second dataset (D2) included HBL and TL measurements of the subset of D1 from the EBD-CSIC collection (N = 51 for HBL, N = 52 for TL). The third dataset (D3) was a subset of D2 that included HFL and EL measurements from adults aged based on their dentition or gonads measured by a single collector (MTRH, N = 35 for HFL, N = 34 for EL). With the creation of D2 and D3 we tried to reduce potential observer bias, live/postmortem stage variation, and ontogeny-associated errors [46, 47]. The fourth dataset (D4) included external measurements (weight, HBL, TL, and EL) of adults aged by any of the methods previously described plus all individuals from Sabah Museum collection, for which we did not have information regarding age (N = 374, pregnant females were excluded). The fifth dataset (D5) included HFL measurements taken including nails (following the U.S.A. procedure: only those coming from EBD-CSIC and U.S.A. collection specimens) from adults aged by any of the methods previously described (N = 167). The analyses of dataset D4 and D5 are included in the as complementary evidence to the main text findings. The sixth dataset (D6) was composed of skull measurements of 136 adult individuals aged based on dentition from several collections and taken by A. H.

Finally, we measured hair length on the dorsum at the rump (HRL) and scapula (HSL) in 17 adult dry skins from Mt. Tambuyukon in order to examine potential pelage differences along elevation. We followed [48] and measured under a magnifier lamp overfur length on the dorsum near the rump and scapula levels by placing a ruler at a right angle to the skin surface and recording the approximate mark where ends of the bunched hairs rested. As [48] pointed out, “the technique is unsophisticated and the results imprecise, but still provide a descriptive estimate of lengths for those pelage constituents”. Pictures of the pelage were taken of several specimens with a Zeiss SteREO Discovery.V8 at 5X to exhibit graphically the conspicuous changes found in hair density and hair thickness along elevation ().

Inclusivity in global research

Additional information regarding the ethical, cultural, and scientific considerations specific to inclusivity in global research is included in the Supporting Information ()

Statistical analyses

To test Bergmann’s rule, we tested the effect of elevation over the following variables, as proxies of body size: weight (W), head-and-body length (HBL), and condyle-incisive length (CIL). We fit a linear model for each one of the cited variables. We standardized ‘elevation’ to include it as a quadratic term in the models of W and CIL given that these variables followed a U-shaped distribution along elevation in a preliminary plot. Standardization was done by subtracting the mean of ‘elevation’ and by dividing by its standard deviation (we used the function scale() to standardize ‘elevation’). ‘Elevation’ was added as a linear and non-standardized term in the HBL model after showing HBL followed a linear relationship with elevation. ‘Sex’ was included as a factor of interest in all the models to test for sexual dimorphism. All models initially also included ‘mountain’ and ‘year’ of collection to control for their effect, but these two variables were excluded from the final models when their effect over the response variable was not significant (see final model formulation in ). Homoscedasticity and normality of residuals were checked by visual exam of scatterplots. Significance was evaluated with an F-test with the function Anova() available in the R car package [49]. Pairwise comparisons between sexes were done with the Tukey post hoc test using the function emmeans() from the emmeans package [50]. We analyzed D4 and D5 in the same manner (). All statistical analyses and graphs were done in R 4.0.3 [51], using RStudio Version 1.3.959 [52].

To test Allen’s rule, we tested the effect of elevation over relative tail length (TL), hindfoot length (HFL), and ear length (EL). We also investigated variation of diet-associated traits across elevation by analyzing zygomatic breadth (ZB), rostrum length (RL), and upper tooth-row length (UTL). We fit Linear Models for each one of the cited variables. ‘Elevation’ was standardized and included as a quadratic term in the models of TL and UTL, and as a linear and non-standardized term in the models of HFL, EL, RL, and ZB as explained above. ‘Sex’ was included in all the models to test for sexual dimorphism, as well as ‘mountain’ and ‘year’ to control for their effect. These two last variables were excluded from the final models when their effect over the response variable was not significant (see final model formulation in ). To account for body size and control for allometric covariance, we also included HBL in the models of TL, HFL, and EL, and CIL in the models of RL, ZB, and UTL. Model quality checking, evaluation of significance, pairwise comparisons between sexes, and analyses over D4 and D5 were performed as in Bergmann’s rule analyses (S1 Table; S1 File).

We tested for changes in individual body condition along elevation with the Scaled Mass Index (SMI), following [53, 54]. SMI adjusts the mass of all individuals to the mass they would have if they had the same body size. The SMI was computed for each individual of D2 as follows: SMI = Wi*(HB0/HBi)b; where Wi and HBi are the weight and the head-and-body length of the individual, respectively, HB0 is the arithmetic mean value of head-and-body for the whole study population, and b is the slope estimate of a standardized major axis (SMA) regression of log-transformed weight on log-transformed head-and-body (b = 1.8399). We fit a Linear Model with SMI as response variable, ‘elevation’ and ‘sex’ as explanatory variables of interest, and ‘mountain’ and ‘year’ as factors to control for. Model quality checking and evaluation of significance were performed as previously described.

Results

Bergmann’s rule

W and CIL followed a significant U-shaped distribution along elevation (D1; W: F1,79 = 4.6719, P = 0.0337; D6; CIL: F1,102 = 8.7565, P = 0.0038), (adjusted R2 (R) and F-statistic for the models of W and CIL: W R = 0.2498, F4,79 = 7.909; CIL R = 0.2363, F3,102 = 11.83). This trend of W was also observed when analyzing D4 (). The smallest average values of W per vegetation zone were observed in the upper montane zone (mean = 123 g, standard error (SE) = 2.14), and in the lower montane zone (mean = 129 g, SE = 3.55) while W values were higher on average in the lowland zone (mean = 136, SE = 5.02) and in the subalpine zone (mean = 139, SE = 8.24). HBL did not show a significant trend over elevation in the small dataset analysis (D2; F1,48 = 0.1862, P = 0.6680), but did follow a significant U-shaped distribution along elevation in D4 analysis () (HBL model R = -0.0326, F4,79 = 7.909) (Figs and ).

Elevational distribution of body size.

Distribution of body size, approximated by (A) weight (W, in g) for D1 individuals and (B) head-and-body length (HBL, in mm) for D2 individuals across the studied elevational gradient. Vegetation zonation is indicated with vertical dotted lines. (A) Dot relative shading corresponds to the year of collection. Weight (W) followed a significant U-shaped distribution along elevation (P = 0.0337), with no significant differences between sexes (P = 0.8292), but a significant trend to increase over the years (P<0.0001) (W model R = 0.2498, F4,79 = 7.909). (B) Black dots represent females, gray dots represent males. Head-and-body length (HBL) was not significantly affected by any of the response variables studied (HBL model R = -0.0326, F4,79 = 7.909).

Elevational distribution of craniodental variation.

Distribution of skull measurements taken from D6 individuals across the studied elevational gradient. Vegetation zonation is indicated with vertical dotted lines. (A) Condyle-incisive length (CIL), (B) rostrum length (RL), (C) upper tooth-row length (UTL), and (D) zygomatic breadth (ZB) across elevation. The limits of the vegetation zones are marked with vertical dotted lines. Black dots: females; gray dots: males. Lines show the relationship of the dependent variable with elevation for both sexes (only one line in the plot), for females (black lines) and for males (gray lines). (A) CIL followed a significant U-shaped distribution along elevation (P = 0.0038). Females showed significantly smaller CIL values than males (P = 0.0096) (CIL model R = 0.2363, F3,102 = 11.83). (B) RL significantly decreased linearly with elevation (P = 0.0148). No significant differences were found between sexes (RL model R = 0.8127, F3,94 = 141.3). (C) UTL followed a significant U-shaped distribution along elevation (P = 0.0378). No significant differences were found between sexes, but differences among mountains were significantly supported (P = 0.0167) (UTL model R = 0.818, F6,98 = 78.93). (D) ZB significantly increased linearly with elevation (P = 0.0004). Females showed significantly lower ZB values than males (P = 0.0079). ZB significantly increased over the years (P = 0.0177) (ZB model R = 0.3857, F4,101 = 17.48). CIL had a significant effect over RL (P< 0.00001), UTL (P< 0.00001) and ZB (P< 0.00001).

Females showed significantly smaller CIL values than males (F3,102 = 6.9790, P = 0.0096; female mean CIL = 44.1 mm, SE = 0.20, sample size N = 45; male mean CIL = 44.6 mm, SE = 0.18, N = 61), while no significant differences between sexes were found for W or HBL (W: F1,79 = 0.0469, P = 0.8292; HBL: F1,48 = 0.1496, P = 0.7006). Similar results were obtained regarding sexual dimorphism when analyzing D4 ().

W followed a significant trend across years, with the most recently collected specimens increasing in W (F1,79 = 17.3116, P<0.0001); this result was also observed in the larger D4 dataset. HBL also significantly increased across years in the larger D4 dataset (). CIL did not change significantly over the years.

Differences among mountains in body size proxies were unsupported for the smaller datasets (D1 and D2), but these differences were statistically significant for W and HBL in the larger D4 analyses ().

Allen’s rule

The relative size of the appendages followed different patterns along elevation (). TL followed a significant U-shaped distribution (F1,46 = 6.2188, P = 0.0163, TL model R = 0.353, F4.46 = 7.821), although males followed a linear decrease (). This was confirmed by the larger dataset D4, where a significant decrease in TL over elevation was observed (). EL did not follow a significant trend with elevation (, F1,30 = 0.0806, P = 0.7784, EL model R = -0.0895, F3,30 = 0.0966; D4, ). HFL globally significantly increased linearly with elevation (, F1,30 = 12.954, P = 0.0011, HFL model R = 0.5665, F4,30 = 12.11), but showed significant differences between Mt. Kinabalu and Mt. Tambuyukon (F1,30 = 11.261, P = 0.0022): Mt. Kinabalu individuals had longer HFL (mean HFL = 44.5 mm, SE = 0.33, N = 22) than Mt. Tambuyukon ones (mean HFL = 42.6 mm, SE = 0.44, N = 13; HFL means averaged over sexes), and followed opposite trends along elevation.

Elevational distribution of body shape.

Distribution of (A) tail length (TL) of D2 individuals, (B and D) hindfoot length (HFL) of D3 individuals, and (C) ear length (EL) of D3 individuals (in mm) over the elevational gradient studied. Vegetation zonation indicated with vertical dotted lines. Black dots: females; gray dots: males. Lines show the relationship of the dependent variable with elevation for females (black lines) and for males (gray lines). (A) Tail followed a U-shaped distribution over elevation (P = 0.0163). Females had significantly shorter tails than males (P = 0.0321) (TL model R = 0.353, F4.46 = 7.821). (B and D) Although HFL globally significantly increased linearly over elevation (P = 0.0011), significant differences were observed between Mt. Kinabalu (B, N = 22) and Mt. Tambuyukon (D, N = 13, P = 0.0022). Females had significantly shorter hindfeet than males (HFL model R = 0.5665, F4,30 = 12.11). (C) No significant trends were observed for EL. HBL had a significant effect only over HFL (P = 0.0018).

Sexual dimorphism was statistically supported in TL (F1,46 = 4.8827, P = 0.0321) and HFL (F1,30 = 7.639, P = 0.0097). Females had shorter tails and hindfeet than males (female mean TL = 136 mm, SE = 2.49, N = 21; male mean TL = 142 mm, SE = 2.12, N = 30; female mean HFL = 42.8 mm, SE = 0.40, N = 15; male mean HFL = 44.3 mm, SE = 0.35, N = 20; HFL was averaged across mountains). No significant differences between sexes were observed for EL (F1,30 = 0.0670, P = 0.7975). No significant differences in EL and TL were found between mountains, although differences in TL among mountains were shown in D4 analysis ().

HBL showed a significant effect over HFL (F1,30 = 11.669, P = 0.0018) but not over EL (F1,30 = 0.1985, P = 0.6591) or TL (F1,46 = 3.1997, P = 0.0802). When analyzing the more extensive datasets D4 and D5, HBL had a significant effect over HFL, EL, and TL ().

Diet-associated and insulation traits

Diet-related craniodental variables were distributed differently along elevation (, ). While RL decreased linearly with elevation (F1,94 = 6.1665, P = 0.0148), UTL followed a U-shaped pattern (F1,98 = 4.4334, P = 0.0378), and ZB increased linearly (F1,101 = 13.6604, P = 0.0004) (RL model R = 0.8127, F3,94 = 141.3; UTL model R = 0.818, F6,98 = 78.93; ZB model R = 0.3857, F4,101 = 17.48). The three diet-related craniodental variables were significantly affected by CIL (RL F1,94 = 314.0722, P< 0.00001; UTL F1,98 = 271.1278, P< 0.00001; ZB F1,94 = 33.5914, P< 0.00001).

Significant sexual dimorphism was observed in ZB (F1,101 = 7.3389, P = 0.0079; female mean ZB = 25.6 mm, SE = 0.09, N = 45; male mean ZB = 26.0 mm, SE = 0.08, N = 61). As with HFL, the small difference among sexes is unlikely to be biologically significant. Significant differences were observed among mountains only in UTL (F2,98 = 4.2664, P = 0.0167), with a significantly smaller UTL on Mt. Kinabalu than on Mt. Tambuyukon (posthoc Tukey test estimate = -0.23, df = 98, P = 0.0275; Mt. Kinabalu, mean UTL = 24.85 mm, SE = 0.06, N = 77; Mt. Tambuyukon, mean UTL = 25.08 mm, SE = 0.07, N = 23; Mt. Trus Madi, mean UTL = 25.13 mm, SE = 0.16, N = 5, UTL means averaged over sexes). A significant trend across years was observed only for ZB, with higher values in the more recent individuals (F1,101 = 5.8175, P = 0.0177).

SMI did not significantly change over elevation (F1,45 = 1.2575, P = 0.2681) or between sexes (F1,45 = 0.0856, P = 0.7712) (Fig 3 in ). No significant differences were observed among mountains (F2,45 = 0.8265, P = 0.4441) or years (F1,45 = 1.1595, P = 0.2873) (model R = -0.0190; model F5,45 = 0.8132) ().

Hair length (HRL, HSL) significantly increased with elevation (R = 0.5618, F1,15 = 21.52, P = 0.0003; R = 0.5033, F1,15 = 17.21, P = 0.0009, ).

Elevational distribution of dorsum hair length.

Distribution of dorsum hair length (in mm) at the rump (HRL) and scapula (HSL) of Mountain Treeshrews over an elevational gradient (in m) on Mt. Tambuyukon. The limits of the vegetation zones are marked with vertical dotted lines. Hair length (HRL, HSL) significantly increased with elevation (HRL model R = 0.5618, F1,15 = 21.52, P = 0.0003; HSL model R = 0.5033, F1,15 = 17.21, P = 0.0009).

Discussion

Bergmann’s rule was not supported in our tropical, montane (elevation) system. Skull length, head-and-body length, and weight decreased from the lowest elevations (<1000 m) to middle elevations (2000–2500 m), then increased from middle to higher elevations, showing a pattern different from that predicted by this rule. This is in contrast with sympatric Bornean shrews which increase with size along elevation [55], but mirrors other studies along tropical elevational gradients such as the New Guinean rats [31], Neotropical Soft Grass mice [32], Wood-wrens [56], and Torrent ducks [36] or subtropical small mammals [6] which were also inconsistent with this rule, with individuals generally becoming smaller with elevation. Neotropical Howler monkeys, Water and Woolly opossums, peccaries, and Crab-eating Fox also followed an inverse Bergmann’s rule pattern across latitude at an interspecific or intraspecific level, although the Crab-eating Fox also showed a strong Bergmannian pattern south of the Equator [33, 57–60].

In the Mountain Treeshrew, individuals from the upper montane forest (2000–2500 m elevation range) exhibited the smallest size and weight. A similar trend was shown for the Yunnan Field Mouse (Apodemus ilex) in Southern China, which decreased in size with elevation up to 3200 m and then exhibited a strong positive relationship with elevation [6]. The previously documented negative correlation between temperature and primary productivity with elevation found in Mt. Kinabalu [42], suggests alternative selective forces may drive these phenotypic changes. Different environmental or biotic factors such as increasing rainfall and water surplus, reduced visibility generated from thick mist, lowered nutrient availability, decreasing radiation, changes in food quality/abundance, and increasing intraspecific competition (the highest population density is found in the upper montane forest) might be interacting and shaping this mid-elevation body size decrease (see [61] for a scheme on the complex interactions that can govern body size in island mammals). Similarly, the unexpected increase in size in the open stunted forest and subalpine meadows of Mt. Kinabalu could be due to: a predator release effect, a reduction of interspecific competition (since Mountain Treeshrews share that habitat with just five other small mammal species [39]), or an increase in relevant food resources. In any case, testing these alternative mechanistic hypotheses is beyond the scope of this study and will require a combination of additional ecological, physiological, and diet data.

In this study system, only relative tail length was consistent with Allen’s rule, which decreased with elevation (in the larger D4 and females of D2), as the rule predicts. Ear and hindfoot length did not follow this rule, since the former remained constant, and the latter increased with elevation (although differences were found among mountains). Allen’s rule has previously been supported for relative tail length in Xeric Four-striped Grass Rat (Rhabdomys pumilio) across elevation; and Virginia opossums (Didelphis virginiana), Herb Field mice (Apodemus uralensis), and Long-tailed macaques (Macaca fascicularis) across latitude; but not for South African Mouse shrews (Myosorex varius) or Soft-furred Tree mice (Typhlomys cinereus) across elevation, or fossorial Southern African Pouched mice (Saccostomus campestris) across latitude [3, 7, 26, 62–64]. Interestingly, tail length seems to exhibit a stronger relationship than ear or hindfoot length with elevation [3] and latitude [65, 66]. The factors driving this stronger relationship might be biological (e.g. temperature) or non-biological ones such as those derived from lower measurement errors from the larger appendages (i.e. tail) compared to the smallest (i.e. ear and hindfoot) when instruments with the same precision are used. In fact, [65] found that temperature poorly predicted tail length (or other extremities) and hypothesized that hindfoot length likely relates to arboreality in the tropics. Longer tails enhance aerodynamic performance and greater maneuverability, which are important to escape predation by leaping, and to improve foraging [65, 67]. Tail length has been shown to be related to scansoriality in squirrels [67] and cricetids [68], but it has only been hypothesized (but not tested) in treeshrews [69]. In the oak-dominated (Fagaceae) upper montane forest of Borneo the main resources provided in the canopy might be hard shell nuts, which treeshrews cannot feed on [40, 41]. This could perhaps encourage more terrestrial foraging at higher elevations, mirroring the pattern shown along latitude in New World small mammals [65]. Thus, the trend observed in treeshrew tail length might be the result of two converging selective agents (temperature and locomotion).

The rule reversal pattern shown in hindfoot length along elevation could also be consistent with more terrestrial locomotion and foraging at higher elevations (longer hindfoot) and a more scansorial fruit foraging behavior at lower elevations (short hindfoot and longer balancing tail). However, most of this study’s Mountain Treeshrews were sampled in the ground, perhaps suggesting a predominantly terrestrial foraging with different degrees of scansoriality across elevation [39]. Similarly, highland bird communities of our study area were also shown to contain species with relatively longer tarsus than lowland relatives [35]. Relative rostrum length decreased significantly with elevation, but it visually followed the same U-shaped pattern with elevation that CIL () and was highly correlated with it, casting doubts on whether the trend is biologically significant or if it might be driven by allometric covariation. Similarly, UTR length followed the same pattern with elevation as CIL, so this trait might also be driven by allometric covariation. UTR is a developmentally constrained trait [70], so stasis is expected, particularly in treeshrews which are not considered dietary specialists [71]. Phenotypic stasis in rostrum length would contrast with other studies [17].

Zygomatic arches are where masseter muscles attach and are integral to the biomechanics of mastication. In Mountain Treeshrews the zygomatic breadth increases along elevation, which could indicate adaptation to a diet with harder food items [72-74]. The observed increase in zygomatic breadth parallels the decrease in availability of soft fleshy fruits and some invertebrates with elevation [41, 75]. Further diet studies will be required to address this hypothesis.

While body size changed along elevation, body condition, an important trait that can be considered a surrogate of lipid content and is characteristic of environmental stress [54], did not. This lack of body condition variation might suggest successful local adaptation throughout its elevational range, perhaps due to generalist foraging strategies [4].

Sexual size dimorphism in the Mountain Treeshrew was statistically significant for some traits (skull length, zygomatic breadth, tail and hindfoot length), which contrasts with other studies on different treeshrew species [41, 45]. The observation of morphological differentiation on a backdrop of high geneflow [76] further highlights the possibility that this pattern is caused by phenotypic plasticity or strong divergent selection, which are some of the mechanisms that could maintain local adaptation despite gene flow [77-81]. Similar patterns of phenotypic differentiation with shallow genetic divergence have been observed in mountain chickadees [82] and Eurasian tree sparrows [83]. Additional novel, functional genomic, transcriptomic and/ or epigenetic approaches will be needed to address the specific molecular basis of these phenotypic traits.

Overall body size may be related to compliance with Bergman’s and Allen’s rules. Several studies [8, 10] have found a significantly lower tendency to conform to Bergmann’s rule within smaller mammals, more specifically, those with a body mass under 500 g, which is the case of the Mountain Treeshrew. Small mammals show alternative strategies to overcome changes in temperature and resource availability than larger ones. Small mammals are better at exploiting torpor, food caching, and microclimatic refugia through burrowing or nesting [84]. Mountain Treeshrews are diurnal and can behaviorally thermoregulate with the use of burrows when necessary, which could reduce the selective pressure of temperature on them as compared to nocturnal species or larger species that are more exposed to rainy and cold weather. An additional mechanism that may mitigate the effect of temperature on body size in this system is conspicuous, possibly adaptive changes in insulation across elevation [85]. High-elevation individuals have softer, more wooly underfur, with a conspicuously higher hair density and longer hairs than lower elevation Mountain Treeshrews (; ), similar to that found in sympatric Crocidura foetida sensu lato-C. baluensis [55], Peromyscus maniculatus [19], and New Guinean Rattus sp. [31].

A highly significant temporal effect was detected in several traits (in D1, D4 and D5 datasets). A temporal trend of increasing weight and head-and-body length was supported, but the magnitude of such changes is difficult to quantify due to the effect of other interacting variables such as elevation. Similarly, ear length decreased and hindfoot length increased with time, while tail length and body condition did not. Given that external measurements are prone to observer bias we remain highly skeptical of any temporal trends since these have not been verified with a less error-prone dataset as in the elevation analyses [46, 47]. Furthermore, skull length, which is exempt of such associated errors, did not increase across years. This lack of body size change across the last century represents the first temporal body size trend described on treeshrews. It mirrors that of other species such as Cape hares and Golden jackals in Israel [86], California Ground squirrels [71], Smith’s Red-backed voles [87] and all but 6 of the 52 populations of 22 species of carnivores analyzed in [88], but is in contrast with other species such as Northern chamois [89] and African Wild dogs [90] which have been shown to decrease in size or Belding’s Ground squirrels and Common Golden-mantled Ground squirrels [68], Large Japanese Field mice [87], Masked shrews [91] and Spain and Israel’s Red foxes [86, 92] which have been shown to increase in size. In the light of the average temperature increase of the last century, this lack of temporal body size change is consistent with our elevational findings, which suggest that temperature is not an important force shaping the morphological variation of the Mountain Treeshrew. It should also be noted that most studies assessing temporal body size variation are biased toward higher latitudes, where there has been a greater temperature increase than in the tropics during the last century [93]. Additional studies will be required to address whether temporal body size changes are less common in the tropics. The absence of an observed change in body size for CIL may also be the result of low statistical power due to the smaller sample size for this trait than HBL and W rather than indicating tolerance of climatic changes [94]. Relative zygomatic breadth increased through the years, perhaps driven by vegetation zonation shifts and associated changes in resource availability. Future studies including newly collected and additional museum specimens will be required to address any of these temporal trends with higher confidence.

These data from a tropical system both support and refute these long-standing biogeographical rules in ways that might shed light on conflicting observations recorded in the literature. Small mammals may experience their environment in a fundamentally different way through their ability to modify their own microenvironment. There also may be essential differences in the environmental clines generated by latitude vs. elevation. One obvious difference is scale, which impacts population size and gene flow and thus the efficiency of selection. Variability in the environment, such as seasonality, also has a major impact on selection. Small mammals in more seasonal mountains might comply more with Bergmann’s rule since resource-shortage periods at higher elevations are longer in colder areas. An increase in size (and fat reserves) in these colder habitats might allow individuals to cope better with these periods. The trends shown in tropical New Guinean rats, Mountain Treeshrews, and subtropical Chinese small mammals suggest a lack of support for Bergmann’s rule in tropical/ subtropical mountain small mammals. Understanding morphological patterns across elevational gradients can provide important insights on how species adapt to changes in different environmental factors and can thereby shed light on how montane specialists could react under a climate change scenario. Overall, our study suggests that selective pressures other than temperature, perhaps driven by diet, food acquisition, predation pressure and/or intra- and inter-specific competition, are playing an important role shaping the morphological variation of the Mountain Treeshrew.

Table of specimens examined in this study, its associated information and morphological data.

(XLSX)

Click here for additional data file.

Appendix.

Methods, results and figures of dataset D2 statistical analyses, figure of body condition, and pictures highlighting differences in fur of Mountain Treeshrews (Tupaia montana) along elevation.

(DOCX)

Click here for additional data file.

Checklist.

Additional information regarding the ethical, cultural, and scientific considerations specific to inclusivity in global research.

(PDF)

Click here for additional data file.