Predictions and test of multiple climate-species richness hypotheses to explain the spatial distribution of tenebrionid beetles in mountain environments

Abstract

Few studies have evaluated how climate is mechanistically related to species richness in mountain environments. We used path analysis to evaluate predictions of several mechanistic hypotheses based on their hypothesized mechanism relating climate with richness of darkling beetles (Coleoptera: Tenebrionidae). We modeled the influence of spatial covariation on climatic variables and tenebrionid richness. Results showed that richness peaks at mid elevations, chiefly influenced by precipitation and temperature, both directly and indirectly through geographic range sizes. The best fitting model explains 84% of the variance of tenebrionid richness. We suggest this pattern is induced by a water-energy balance along the altitudinal gradient. At low elevations, energy availability is high but water deficit may limit species richness; in contrast, at high elevations water availability is high, but energy deficit may limit species richness. These results suggest high susceptibility of the study region to future global climate change.

Key words
darkling beetles; species diversity; elevation; water-energy balance

INTRODUCTION

Several hypotheses relate climatic variables with species richness. One that has received much attention is the “species–energy hypothesis” (Wright 1983, Ballesteros-Mejía et al. 2017). Derived from species–area theory, this hypothesis posits that energy availability generates and maintains species richness, resulting in a positive, saturating relationship between energy availability and species richness. Another widely accepted hypothesis is the so-called species–climate theory (Rohde 1992), which posits that climatic conditions limit species richness because organisms can live only within a range of physiologically tolerable climatic conditions. The are three main versions of the species-climate hypothesis: (a) the “climate-variability hypothesis” (Janzen 1967), which asserts that a monotonic increase in climate variability leads to a negative relationship with species range sizes, which in turn may influence species richness (see Stevens 1989, Suggitt et al. 2015), (b) the “climate extreme hypothesis” (von Humboldt 1808), which asserts that the pool of species available for colonization declines with increasing harshness of environmental conditions (Moser et al. 2005), and (c) the “optimal-climate hypothesis” (Pither 2003, Bozinovic & Naya 2015), which asserts that species may experience some optimum level of climatic conditions that lead to a positive relationship between range size and richness. Although these climate-based hypotheses share many similarities, they differ importantly in the underlying mechanisms of climatic regulation of species richness.

Unfortunately, most previous studies evaluating these hypotheses have rarely evaluated the underlying causal mechanisms responsible for them. In contrast, these analyses are based only on correlative procedures that merely suggest causal mechanisms among highly correlated variables (e.g., Currie 1991, Hawkins et al. 2003, Werenkraut & Ruggiero 2014, Werenkraut et al. 2015, Yu et al. 2016). These multivariate procedures present two additional limitations: the collinearity between variables may hamper the identification of the effect of explanatory variables when ecological inferences are the target of study (i.e., different hypotheses; see, Legendre & Legendre 1998) and, most importantly, this methodology does not provide an explicit test of alternative explanations (Rahbek et al. 2007). Therefore, an accurate understanding of richness-climate relationships implies the simultaneous evaluation of alternative models that explicitly represent hypothesized mechanisms (Currie et al. 2004). The simplest way to achieve this sort of evaluation is to represent mechanistic hypotheses as causal models. Causal models provide estimates of the magnitude and significance of hypothesized causal connections between sets of variables and allow the comparison of competing models.

When applying causal models, care must be taken to avoid falsely attributing causality to environmental variables, which may be in fact correlated because of a common spatial structure present in both the predictor and response data (Legendre 1993). Thus, an empirical test of these causal models requires the evaluation of the spatial structure in which they are embedded. Partial regression analysis allows the examination of the dependence of species richness on different groups of variables by decomposing total explained variance into environmental variance, spatial variance, variance shared by environment and space and unexplained variance (Legendre & Legendre 1998).

Here, we compare multiple climate-based hypotheses to evaluate alternative mechanisms explaining spatial variation in darkling beetles (Coleoptera: Tenebrionidae) in montane environments in Argentina. Tenebrionids are the most representative group of insects in these environments. They are mostly wingless (low dispersal capability), and are characterized by narrow distributional ranges and high sensitivity to climatic variability (Carrara et al. 2011). Using data on geographic ranges and abundance, we pursue two objetives: (1) to evaluate mechanistic hypotheses relating environmental variables with darkling beetle richness throughout altitudinal environments, and (2) to assess the spatial structure of darkling beetle richness.

Hypotheses, mechanisms, predictions and climatic variables

Species–energy hypothesis (SEH)

This hypothesis has also been called the “more individuals hypothesis” (Srivastava & Lawton 1998). It proposes two different underlying causal mechanisms linking energy availability with species richness, depending on how energy availability is measured (Evans et al. 2005). First, solar energy metrics (i.e., temperature) influence species richness by enhancing organisms metabolic rates, especially in ectotherms, promoting large populations with low extinction vulnerability (Allen et al. 2002). Second, productive energy metrics (i.e., the amount of resources available for consumers to turn into biomass) influence species richness by allowing a greater number of individuals in more productive sites, again promoting large populations with low extinction vulnerability (Wright 1983). We use the Enhanced Vegetation Index (EVI) as a measure of productive energy in an area, which is a satellite derived measure of the greenness of vegetation used in previous studies (see Waring et al. 2006, Qian et al. 2007). As tenebrionid beetles are mostly phytophagous insects, we assume that EVI will represent a reliable measure of plant productivity to detect relationships between tenebrionid population sizes and richness.

Modified species–energy hypothesis (SEH_M)

Recently, Carrara & Vázquez (2010) modified the original species–energy hypothesis to include energy variation. They argued that areas not only with low quantity of resources but also with high temporal variation in resources will have low species richness (see also Ruggiero & Kitzberger 2004). The mechanism proposed is that widely-fluctuating population abundance will result in high extinction risk. We will evaluate this modified version of SEH (hereafter SEHM), using measures of solar and productive energy corrected by their annual coefficient of variation.

Under SEHM we predict that solar-energy, measured as average annual temperature corrected by its annual coefficient of variation (T_TV), has a positive effect on N, which in turn has a positive effect on S. Finally, we predict that productive energy in an area, measured as mean annual EVI corrected by its annual coefficient of variation (ETV), has a positive effect on N, which in turn has a positive effect on S (Figure 1b).

Figure 1
Path diagrams representing the relationship expected by each climate-based hypothesis (a to f). Full arrows represent positive relationships, dotted arrows negative relationships, and double arrows covariance. S, species richness; N, individual abundance; RS, species range size; T_A, average annual temperature; P_A, average annual precipitation; E_A, average annual normalized difference vegetation index; T_TV, average annual temperature corrected by its annual coefficient of variation fitted to N; E_TV, average annual EVI corrected by its annual coefficient of variation fitted to N; T_R, annual temperature range; P_R, standard deviation of annual precipitation; T_M, annual minimum temperature; P_M, annual minimum precipitation.

Climate variability hypothesis (CVH)

This hypothesis posits that climate variability influences species richness because relatively few species should be physiologically equipped to tolerate climatic variations (Moser et al. 2005). Species that tolerate large magnitudes of climatic variability should occupy larger geographical areas than species with narrow climatic tolerances (Pither 2003). Conversely, low climatic variability promotes niche differentiation (i.e., evolutionary adaptation and speciation; Currie 1991), inducing higher species richness in stable sites.

According to CVH, we predict that large variability in temperature, measured as annual temperature range (T_R), and in precipitation, measured as annual precipitation range (P_R), will have a positive effect on species range size (i.e., areas of occupancy, hereafter RS), which in turn has a negative effect on S because not all species can tolerate climatic variability. We also predict that T_R and P_R will have a negative direct effect over species richness because variability retards evolutionary adaptation (Figure 1c).

Climatic extremes hypothesis (CEH)

This hypothesis suggests that the pool of species available for colonization declines with increasing climatic extremes (Moser et al. 2005). Thus, climatic extremes within the range of a species’ distribution may relate to its range size in a way analogous to that proposed by CVH: species that are able to tolerate extreme climatic conditions at any locality are climate-generalist, which allows them to occupy larger geographical areas than species with narrower climatic tolerances (Pither 2003). Also, extreme climatic conditions promote the development of physiological adaptation to survive; thus, because these adaptations are very costly, sites with strong climatic extremes have lower species richness than sites with weak climatic extremes (Turner et al. 1987).

Under CEH, we predict that extreme temperature, measured as minimum temperature of the coldest month (T_M), and extreme precipitations, measured as minimum precipitations of the driest month (P_M), have a negative effect on RS, which in turn has a negative effect on S, because not all species can tolerate climatic extremes. We also predict that T_M and P_M have a positive direct effect on species richness because climatic extremes retard physiological adaptation (Figure 1d).

Optimal climate hypothesis (OCH)

This hypothesis suggests that species richness tends to be greatest where the potential amount and duration of biological activity is greatest (i.e., optimum), decreasing as the climatic potential for biological activity decreases (Mourelle & Ezcurra 1996). Species that inhabit in this climatic optimum are unable to tolerate climatic conditions outside of this optimum, which allows them to occupy smaller geographical areas than species with broader climatic tolerances (Pither 2003, Bozinovic & Naya 2015). However, optimal climatic conditions may influence species richness directly by physiological adaptations, because it is easier to adapt under optimal habitat condition (Turner et al. 1987, Mourelle & Ezcurra 1996).

According to this hypothesis, we predict that T_A, E_A and annual precipitation (P_A) will have a negative effect on RS, which in turn has a negative effect over S, because most of species that can tolerate climatic optimum have small range sizes. Also, we predict that T_A, E_A and P_A have a direct positive effect over S (Figure 1e).

This hypothesis can also be represented by a quadratic function (O´Brien et al. 1998), where the optimum is modeled as a unimodal curve. Thus, we can modify OCH to incorporate these quadratic effects (OCHQ). Under OCHQ, we predict that values of T_A, E_A and P_A have positive effects over RS and that T_A², E_A² and P_A² have negative effects over RS, which in turn has a negative effect over S. Also, we predict a direct positive effect of T_A, E_A and P_A on S and a direct negative effect of T_A², E_A² and P_A² on S. Because quadratically transformed variables are functions of the original (untransformed) variables, we expect covariation between them (see, Scheiner et al. 2000; Figure 1f).

MATERIALS AND METHODS

Study site

The study area covers both Andean and extra-Andean mountain environments in central Argentina (32° to 35° S). The Andes is the highest mountain chain in the Americas, with elevations that surpass 6000 m. The extra-Andean mountain ranges are physically separated from the Andes, but their uplift was influenced by the Andes orogeny. The climate in this region is arid-semiarid, changing from temperate to cold with increasing elevation, with mean minimum temperatures of -15 °C and mean maximum temperatures of 30 °C. We worked at 34 sampling locations spread over four altitudinal transects and covering the major biogeographic types. The geographic location of sample sites are represented in Figure 2, and the main characteristics of biogeographic types are described below:

Figure 2
Distribution of sampling sites across Andean and extra-Andean mountain environments in central Argentina around parallel 32, 34 and 35.

Andean transects

Aconcagua Provincial Park. This transect rans through the Main Andean Range (Cordillera Principal), following the international road to Chile (national road 7). It included three sampling locations between 2400 m and 3500 m a.s.l, with xerophilous vegetation closely related to elevation and slope exposition. The lowest sampling location, Punta de Vacas (32.86 S, 69.76 W, 2400 m), lies at the first floor of High-Andean biogeographic province. The vegetation is dominated by Adesmia pinifolia, A. remyana and Mulinum spinosum (Ambrosetti et al. 1986). The second location, Morenas de Horcones (Aconcagua Provincial Park, 32.81 S, 69.94 W, 2900 m), is at the second floor of High-Andean biogeographic province. The vegetation is characterized by the presence of tall shrubs with Adesmia pinifolia, Chuquiraga oppositifolia and Tropelum polyphyllum (Martínez Carretero & Méndez 1992). The highest site, Las Cuevas (32.80 S, 70.07 W, 3500 m), encompasses the upper vegetation limit, with pulvinar vegetation dominated by Adesmia subterranea and A. hemisferica.

Laguna del Diamante Provincial Park. This transect also rans through the Main Andean Range following the provincial road 98, and included five sampling points at four elevations. At the lowest elevation, two sampling points were located on north and south slopes of one mountain (34.24 S, 69.38 W, 2400 m). The vegetation is characteristic of the first floor of the High-Andean biogeographic province, and is dominated by shrubs Adesmia pinifolia, A. scneiderii and A. ovata and the grass Poa halciformis. The third sampling point (34.23 S, 69.50 W, 3700 m), lies at the third floor of the High-Andean biogeographic province. The vegetation is characterized by grasses dominated by Poa halciformis, Junellia uniflora and Stipa chrysophylla (Martínez Carretero et al. 1999). The fourth site (34.19 S, 69.53 W, 3300 m), within the second floor of High-Andean biogeographic province, is characterized by dwarf shrubs and dominated by Adesmia subterranea and A. hemisferica. The last site, Vegas del Maipo (34.21 S, 69.70 W, 2700 m), is characterized by peat bogs and dominated by Festuca desvauxii, F. argentina, Eleocharis albibracteata, Oxychloe andina, Festuca argentina and Werneria pygmaea.

Extra-Andean transects

Villavicencio Nature Reserve. This transect followed provincial road 52 through Villavicencio Nature Reserve. The first site (32.62 S, 68.94 W, 1100 m), within the Monte biogeographic province, is characterized by desert vegetation dominated by shrubs, particularly Larrea. The second site (32.51 S, 69.01 W, 1900 m), within Prepuna biogeographic province, is characterized by vegetation closely related to the Monte but differing from it by the absence of Larrea, the high abundance of characteristic cacti and bromeliads and the presence of ferns (Haumman 1947). The third site (32.49 S, 69.06 W, 2800 m), within the Puna biogeographic province, is characterized by a tall grassland and shrub arid steppe dominated by Plazia daphnoides, Lycium decipiens, L. fuscum, Ephedra multiflora, Verbena diversifolia, Baccharis thymifolia, B. polifolia and B. tola (Roig-Juñent et al. 2003). The fourth site (32.72 S, 69.15 W, 3000 m), lies in a flat plain with a grassland dominated by Stipa tenuissima. The final site (32.37 S, 69.44 W, 2400 m), within the High-Monte biogeographic province in Uspallata valley, is characterized by desert vegetation, dominated by shurbs of the genus Larrea.

Cerro Nevado. This transect rans through the extra-Andean mountain Cerro Nevado (3800 m), and included four sampling points. The first point in this transect (35.62 S, 68.54 W, 2400 m) falls within the Patagonian biogeographic province. It is characterized by a shrubby steppe on sandy and basaltic soils, dominated by Neosparton aphyllum and Sporobolus rigens. The second point (35.61 S, 68.52 W, 2600 m) lies at the first floor of the High-Andean biogeographic province. The vegetation is dominated by the shrubs Adesmia pinifolia, A. scneiderii and A. ovata and the grass Poa halciformis. The third point (35.60 S, 68.51 W, 2800 m) falls within the second floor of the High-Andean biogeographic province. It is characterized by dwarf shrubs and dominated by Adesmia pinifolia and Anarthrophyllum rigidum. Finally, the fourth point (35.59 S, 68.50 W, 3000 m), within the third floor of High-Andean biogeographic province, is characterized by low vegetation dominated by Panthacantha ameghinoi.

Sampling methods

Darkling beetles were sampled between December and February 2004-2005 and 2005-2006. At each site we placed 32 pitfall traps for ten days, arranged in eight groups of four traps separated by 30 m from each other. Traps were cylindric plastic containers of 10 cm diameter, with 100 cm³ of a 30% ethylene glycol-water solution to preserve the material. Once traps had been collected, samples were washed and placed in 90% ethanol. Beetles were identified to species by using keys and revisions from Kulzer (1955, 1958, 1962, 1963), Flores (1999, 2004) and Flores & Chani-Posse (2005). Voucher specimens are deposited at the Entomological Collection of the Argentine Institute of Dryland Research (IADIZA, CCT CONICET Mendoza, Argentina).

Climatic variables

Climatic data were extracted from the mean monthly climate database published by Hijmans et al. (2005), which was compiled at 1 × 1 km latitude-longitude resolution for the period of 1961-1990. This source represents one of the most accurate published databases on contemporary climate. The productivity measure used was the enhanced vegetation index (EVI). This measure was extracted from satellite images at 1 × 1 km latitude-longitude resolution, freely available at the NASA web site (http://modis.gsfc.nasa.gov/data). Monthly measures between January of 2003 and January of 2006 were extracted and processed to obtain mean annual measures (E_A) and measures of E_A corrected by its variation (E_V), as explained above. To evaluate the SEH_M model derived in Carrara & Vázquez (2010), we first estimated parameters conducting non-linear regression analyses by using maximum likelihood optimization with the mle2 function of the bbmle package (Bolker 2008) of R statistical software (R Development core team 2009).

Data analyses

Species richness, range sizes and abundances

Because some pitfall traps were destroyed by introduced and native ungulates, we conducted sample-based rarefaction curves for all sites to obtain an accurate estimation of species richness. To construct these curves, we used the freely available software EstimateS V7.5.2 (Colwell 2000). Rarefaction curves were rarefied to the lowest number of samples (14) recorded at a site to ensure valid comparisons of species richness (S) between different sites (Gotelli & Colwell 2001). Although sample-based rarefaction curves are useful to compare species density (i.e., the number of species per unit area) and not species richness, we performed sample-based and not individual-based rarefaction because the latter implies a loss of abundance information, which is necessary to evaluate the mechanisms proposed by SEH and SEH_M.

We used the areas of occupancy of tenebrionid species in different trapping sites as a measure of species range size. To this end, we defined each sampling site in montane environments as a particular “habitat,” according to their biogeographic characteristics. Then, we recorded the presence of each species in each site, and estimated mean range size for each sampling site. Mean range size was defined by sample sites following Gaston & Chown (1999) as the sum of areas occupied by each species that occurred in a site divided by total number of species that occurred in the site. Abundance (N) was measured as total individuals in a site.

Statistical analyses

Path analysis, a type of structural equation modeling (SEM), was used to evaluate the fit of several climate-based hypotheses to the data. Path coefficients, variances and covariances and overall goodness of fit to the data were estimated using maximum likelihood from the variance-covariance matrix among the set of variables included in a particular model (for more specifications see Shipley 2000). These statistics allowed us to estimate the direct and indirect strengths of the effects of climatic variables on species richness and to compare and select between alternative models. The null hypothesis implicit in the model fitting procedure is that there are no differences between observed and predicted covariance matrices; thus χ² p-values greater than 0.05 indicate that it is not possible to reject the model (see Shipley 2000). Because there is no consensus as to the appropriateness of χ² as the best goodness-of-fit index, we also used other alternative indices to assess the fit between the models and the data: the comparative fit index (CFI), the non-normed fit index (NNFI), and Schwartz-Bayes information criterion (BIC). CFI and NNFI are known to perform well, especially when sample sizes are small (Shipley 1997). BIC is based on the model’s χ² value, degrees of freedom, number of free parameters and sample size. Two other goodness of fit indices, GFI (Goodness of Fit Index) and NFI (Normed Fit Index), were not appropriate because of non-normality of the data (GFI) and the small size (NFI).

A basic assumption of path analysis is that data follow a multinormal distribution. This assumption was not met by the data used here, suggesting that ML χ² estimates are biased; to solve this problem we apply the Satorra-Bentler χ² statistic using the package sem.additions (Byrnes 2009) of R statistical software, which accounts for deviations from normality. To circumvent problems resulting from small sample size, we corrected goodness of fit statistics using the swain function (Boomsma & Herzog 2008) of R statistical software, which provides robust fit of NNFI and CFI when samples are small.

In SEM, the strength of the direct causal effect between two variables is represented by a path coefficient. Indirect effects are those of a given variable on another one through a third or more mediator variable, estimated by multiplying the path coefficients along the entire path. The coefficient of determination, calculated as , represents the proportion of variance unexplained by predictor variables (de Lafontaine & Houle 2007). Once the best fitting model had been selected, we evaluated if the model could be improved by eliminating non-significant paths and re-assessing fit. We used the Satorra-Bentler χ² difference (Satorra & Bentler 2001) to select between the original and the reduced models. Significant χ² values favor complex models, while non-significant χ² favor simple models.

Spatial analysis

The variation in a given variable (e.g., richness) may result from spatial autocorrelation of the variable itself, from relationships between the variable of interest and another variable that is spatially structured or from the space-independent relationships between the two variables (Legendre & Legendre 1998). To evaluate the influence of spatial structure on the SEM models, we conducted a trend surface analysis, widely used in ecology to account for spatial autocorrelation among continuous variables (van Rensburg et al. 2002). Trend surface analysis was applied to determine the best combination of spatial variables that explained variation in response variable. The spatial component of the variation in tenebrionid beetles was modeled using a third-order polynomial that extracts linear and more complex spatial features from data:

$$f(x,y)=b_0+b_{1}x+b_{2}y+b_3{x}^2+b_{4}xy+b_{3}x{y}^2+b_6{x}^3+b_7{x}^{2}y+b_{8}x{y}^2+b_9{y}^3$$ (eq. 1)

where x and y are longitude and latitude, respectively (Legendre 1993). As in Legendre & Legendre (1998), after standardizing the variables we performed ordinary least square regression (OLS) with the x and y third degree polynomial term of eq. 1 fitted on tenebrionid richness. Nonsignificant terms were removed from the model, after which fit of the simplified model was assessed again. The resulting trend surface equation was highly significant (R² = 0.59, p < 0.001).

Thereafter, the trend surface was introduced into modeling to interpret the influence of spatial heterogeneity on climatic variables and species richness. To this end, we used partial regression analyses to partition the variation of tenebrionid beetles into four components: (a) the non-environmental spatial fraction (the component of the spatial variation in species richness that is not shared with the environmental variables); (b) the spatially structured environmental fraction (the spatial structuring in the species richness data that is shared with the environmental variable data); (c) the non-spatial environmental fraction (the component of the spatial variation in species richness that can be explained by the environmental variables independent of any spatial structure); and (d) the unexplained (residual) variation (Legendre & Legendre 1998). Notice that the richness coefficient of determination in causal models and the trend suface analysis are both obtained by multiple regression analyses. Thus, variance partition using multiple regressions allowed me to disentangle the real contribution of the climatic variables and spatial components on species richness. These analyses were performed with the freely available statistical software SAM (http://www.ecoevol.ufg.br/sam; Rangel et al. 2010).

RESULTS

A total of 34 darkling beetle species and 2307 individuals were collected during the two years of study (taxonomic position and species abundance are listed in Table I). The fitted causal models for each climate-based hypothesis are shown in Figure 3. SEM results indicate that three models fit the data reasonably well (CEH, OCH and OCHQ; Table II), indicating that the covariance structure of the data did not deviate significantly from the covariance structure of the model. Model selection identified OCHQ as the best fitting model: OCHQ shows the highest values of CFI and NNFI and the lowest value of BIC. OCHQ was also the causal model explaining the greatest proportion of total variance in geographic range size (R² = 61 %) and species richness (R² = 62 %). Removal of two non significant paths from OCHQ, E_A and E_A ², improved model fit significantly (Satorra-Bentler χ² difference p>0.05; Figure 4).

Figure 3
Path analysis results showing standardized path coefficients (letters a to f refer to the same hypotheses specified in Figure 1). Solid arrows represent positive relationships and hatched arrows negative relationships, double arrows covariance and vertical arrows unexplained variance. S, species richness; N, individual abundance; RS, species range size; T_A, average annual temperature; P_A, average annual precipitation; E_A, average annual normalized difference vegetation index; T_VN, average annual temperature corrected by its annual coefficient of variation fitted to N; T_VS, average annual temperature corrected by its annual coefficient of variation fitted to S; E_VN, average annual EVI corrected by its annual coefficient of variation fitted to N; T_R, annual temperature range; P_R, standard deviation of annual precipitation; T_M, annual minimum temperature; P_M, annual minimum precipitation.

Figure 4
Best fitting model (OCHQ, Figure 1f) after elimination of non-significant paths. Arrows and variables as in Figure 1.

Table I
Taxonomic position of total species recorded in study sites: Abundance represents total number of individuals collected throughout sampling.

Subfamily	Tribe	Species	Abundance
Alleculinae	Alleculini	Lobopoda breyeri	1
Lagriinae		Gen. sp. undetermined	1
Pimeliinae	Caenocrypticini	Caenocrypticoides triplehorni	26
	Edrotini	Arthroconus sp. 1	31
		Arthroconus sp. 2	19
		Hylithus kovacsi	4
		Hylithus tentyroides tentyroides	33
	Epitragini	Epitragus sp.	3
		Nyctopetus sp.	3
		Nyctopetus tenebrioides aconcaguensis	88
		Omopheres sp.	3
	Nycteliini	Epipedonota cristallisata	543
		Epipedonota senex	138
		Mitragenius araneiformis	27
		Nyctelia laevis	148
		Nyctelia nevadoensis	367
		Nyctelia plicatipennis	80
		Nyctelia setipennis	12
		Psectracelis semistrigosa	14
		Psectrascelis mamillionea	29
		Scelidospecta lobata	1
	Physogasterini	Gen. sp. undetermined	2
	Praocini	Calymmophorus cucullatus	2
		Falsopraocis australis	336
		Platyholmus uspallatensis	5
		Praocis concinna	13
		Praocis rotundata	27
		Praocis sp. 1	91
	Stenosini	Discopleurus sp.	2
		Ecnomoderes bruchi	1
Tenebrioninae	Opatrini	Blapstinus punctulatus	35
	Scotobiini	Scotobius andrassy	114
		Scotobius wittmeri	14
		Scotobius punctatus	94

Table II
Goodnes of fit used for model selection and analysis of the spatial structure of causal models; SEH, species-energy hypothesis; SEH_M, species energy hypothesis modified; CVH, climatic variability hypothesis; CEH, climatic extreme hypothesis; OCH, climatic optimal hypothesis; and OCHQ, quadratic version of the climatic optimal hypothesis. SB^{- χ}², Satorra-Bentler-^χ²; NNFI, Non-Normed Fit Index (high values indicate good fit); CFI, Comparative Fit Index (high values indicate good fit); BIC, Schwartz-Bayes Information Criterion (low values indicate good fit); R², proportion of variance explained (high values indicate good fit); a, variance explained by nonenvironmental spatial fraction; b, variance explained by spatially structured environmental fraction; c, variance explained by nonspatial environmental fraction; and d, unexplained variation. AIC, Akaike’s Information Criterion used for spatially structured model selection; dAIC, relative model fit (the best fitting model has dAIC = 0).

Model	SB-χ ²	NNFI	CFI	BIC	R2	a	b	c	d	AIC	dAIC
SEH	0.01	-0.18	0.41	0.04	<0.01	0.56	0.03	0.02	0.39	7.78	14.50
SEHM	<0.01	-0.89	0.06	40.79	<0.01	0.56	0.03	0.02	0.39	7.78	14.50
CVH	<0.01	-3.33	0.28	16.13	0.24	0.38	0.21	0.02	0.39	11.04	17.76
CEH	0.05	0.64	0.94	0.18	0.32	0.37	0.22	0.09	0.32	7.87	14.59
OCH	0.12	-0.13	0.66	-4.89	0.21	0.53	0.06	0.15	0.26	4.77	11.49
OCHQ	0.99	0.95	0.98	-42.28	0.62	0.25	0.34	0.25	0.16	-6.72	0.00

When spatial autocorrelation is taken into account, the variance explained by environmental variables decreases substantially. This change in the explanatory power of models stems from the spatial structure of environmental variables (Table II). In the majority of models the percentage of variance explained by spatial structure of environmental variables was always less than half of total variance. However, pure spatial variance explained a high proportion of the total variance while pure environmental variance accounts only for a small proportion of tenebrionid richness. The sole exception to these findings was the OCHQ model (see below). As in path analysis, the OCHQ spatially structured model shows the best fit to the data (see AIC values in Table II). As a whole, this model explains 84% of tenebrionid richness, while the spatially structured environmental variables explain 34% of variance and pure spatial and environmental components each explain 25%. These findings indicate that climatic variables have great predictive explanatory power for species richness, even after accounting for spatial structure.

DISCUSSION

Our study shows an approximation of how climate is mechanistically related to tenebrionid richness in mountain environments. From the results showed in this work it was possible to recognize whether the predictions of different hypotheses that relate climate with richness through some biological mechanisms were supported by data. These findings are particularly interesting given that most of previous researches on beetles diversity (e.g., Botes et al. 2007, González-Megías et al. 2008, Yu et al. 2013, Corcos et al. 2018, Şenyüz et al. 2019) were limited to assess the role of multiple climatic variables directly on richness without testing the hypothesized causal relationships between population abundances or species range sizes with species richness (but see Gebert et al. 2020, Njovu et al. 2021). Therefore, the possibility to recognize the underlying mechanism that operate behind different climate-species richness relationships offer a better understanding of the ways in which climate structures species richness.

Tenebrionid richness showed a striking unimodal relationship with elevation and was strongly influenced by climate (OCHQ model). Thus, species richness exhibits a unimodal relationship with both increasing temperature and precipitation, reaching a peak at intermediate levels of both climatic variables. Overall, quadratic temperature and precipitation accounted for nearly 59% of explained variance of total richness. Previous studies have found that species richness increases with increasing temperature (i.e., a linear, positive relationship; McIntyre 2000, de los Santos et al. 2002, Fattorini 2006), although weak unimodal relationships were found by Botes et al. (2007). In contrast, the relationship between precipitation and tenebrionid richness is contradictory, with studies that reported no relationship (Cepeda-Pizarro et al. 2005), or positive (Yu et al. 2016) and negative covariation between them (McIntyre 2000). Our results help clarify this relationship between tenebrionid richness and precipitation.

The quadratic relationship found between richness and climate suggests a balance between temperature and precipitation (i.e., water-energy dynamics; sensu O’Brien 1998), driving tenebrionid richness across altitudinal environments. Thus, tenebrionid richness appears to be greatest where the potential amount and duration of biological activity is greatest, decreasing as the climatic potential for biological activity decreases. Interestingly, O’Brien (1993) found that species richness could be described by the relationship S= W + (E - E²), where S is species richness, W a measure of water availability (e.g., precipitation), and E a measure of energy availability (e.g., temperature), with a peak of richness at intermediate levels of water and energy availability. Similarly, we found that tenebrionid richness is best described by the OCHQ model, which indicates that the greatest potential for biological activity is found at intermediate elevations, with intermediate levels of precipitation and temperature. However, contrary to O’Brien’s assumptions, in our model species richness is not a positive linear function of water; instead, tenebrionid richness is a quadratic function of water, with tenebrionid richness peaking at intermediate precipitation levels. In fact, tenebrionid beetles appear to be rare or absent in areas at the lower and upper extremes of the precipitation gradient (G. Flores, unpublished results).

Consideration of regional climatic dynamics is useful to understand the climatic basis of our results. The study region includes some of the highest peaks of the Andes, which prevent moist air masses from the west (Pacific winds) from reaching eastern (Argentine) slopes. Highlands receive a greater amount of precipitation, mostly from the west as snow during winter, while lowlands receive precipitation from eastern winds during the summer. So, in this case total precipitation is positively correlated with elevation. In contrast, because air is less dense at high elevations, temperature decreases with increasing elevation (Brown & Lomolino 1998). Lowlands are thus characterized by low precipitation and high temperatures, highlands by high precipitation and low temperatures, and middle elevations by moderate temperatures and precipitation (Figure 5).

Figure 5
Annual temperature (red line) and precipitation (black points and dashed line) throughout the altitudinal gradients analyzed.

Water and energy regimes have been demonstrated to be essential to the physiology of plants and directly influence distribution and species richness (O´Brien 1993, McCain 2007); however, for animal richness patterns the effects of water and temperature may be direct (physiological) or indirect (via resources; see Hawkins et al. 2003). Results of path analyses suggest that the water-energy balance has a direct effect rather than an indirect effect on tenebrionid richness. This direct effect of climate on tenebrionid species is likely determined by their drought, cold and freezing tolerances (see Cloudsley-Thompson 2001). In contrast, path analysis results indicate that the effect of water-energy balance on tenebrionid richness does not operate indirectly through resources.

The analysis including spatial structure demonstrated that the pure climatic component was high (25%), giving strong support to climatic models. Nonetheless, the spatially structured component of the variation in the environmental variables accounted for greater proportion of variation in species richness (34%), indicating similar patterns of spatial autocorrelation between richness and climatic variables (van Rensburg et al. 2002). Thus, spatial structure of climatic variables is closely related to complex spatial heterogeneity of mountainous environments. The pure spatial component (25%) indicates that tenebrionid richness is spatially autocorrelated, independently of the spatial structure of the explanatory variables analyzed. This spatial component could reflect other environmental variables not considered in this study (e.g., population- and community-level processes, interaction with other parts of the biotic community, historic processes; Legendre 1993). Finally, the low proportion of unexplained variance (16%) indicates high predictability of tenebrionid richness (He et al. 1996).

Overall, these findings have important implications to efforts to predict the effect of future climate change on species richness. As is widely known, global climate change poses a number of potential risks to biodiversity by inducing transformation of habitats, which in turn may induce changes in geographic distributions of species (Walther et al. 2005). Here, tenebrionids showed a strong relationship with water-energy balance, and because it is predicted that global warming influences both water and energy regimes, tenebrionids should be strongly affected by anthropogenic climate change. Thus, more research is needed to predict the response of tenebrionids to global change. Also, because tenebrionid richness is highly predictable based on environmental and spatial variables, they could be good candidates as bio-indicators of climate change.

ACKNOWLEDGMENTS

We thank A.H. Hurlbert for exchange of ideas and comments, and L.A. Puppo Bryant for help with the design of figures. This research was supported by Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET, Argentina) and the Agencia Nacional de Promoción Científica y Técnica, Argentina (ANPCyT) through the following projects: PICT#2013-3128; PICT#2013-1539, PIP# 112-201101-00987, PI-IADIZA.

REFERENCES

• ALLEN AP, BROWN JH & GILLOOLY JF. 2002. Global biodiversity, biochemical kinetics, and the energetic-equivalence rule. Science 297: 1545-1548.
• AMBROSETTI JA, DEL VITTO LA & ROIG FA. 1986. La vegetación del paso de Uspallata, Provincia de Mendoza, Argentina. Veröff Geobot 91: 141-180.
• BALLESTEROS-MEJÍA L, KITCHING IJ, JETZ W & BECK J. 2017. Putting insects on the map: near-global variation in sphingid moth richness along spatial and environmental gradients. Ecography 40: 698-708.
• BOLKER BM. 2008. Ecological models and data in R. Princeton Univ. Press, Oxford.
• BOOMSMA A & HERZOG W. 2008. R function swain: Correcting structural equation model fit statistics and indexes under small-sample and/or large-model conditions (http://www.ppsw.rug.nl/~boomsma/swain.pdf).
• BOTES A, MCGEOCH MA & CHOWN SL. 2007. Ground-dwelling beetle assemblages in the northern Cape Floristic Region: Patterns, correlates and implications. Austral Ecol 32: 210-224.
• BOZINOVIC F & NAYA DE. 2015. Linking physiology, climate, and species distributional ranges. In Martin LB, Ghalambor CK & Woods HA (Eds) Integrative organismal biology, 1st edn. Wiley, New York, p. 277-290.
• BROWN JH & LOMOLINO MV. 1998. Biogeography. Sinauer, Oxford.
• BYRNES J. 2009. sem.additions: additional methods for structural equation modelling. R package version 0.1-04/r3. http://R-Forge.R-project.org/projects/sem-additions/
  » http://R-Forge.R-project.org/projects/sem-additions/
• CARRARA R & VÁZQUEZ DP. 2010 The especies-energy theory: a role for energy variability. Ecography 33: 942-948.
• CARRARA R, VÁZQUEZ DP & FLORES GE. 2011. Habitat specificity can blur the predictions of species–energy theory: A case study of tenebrionid beetles adapted to aridity. J A Environ 75: 703-710.
• CEPEDA-PIZARRO J, PIZZAROO-ARAYA J & VÁSQUEZ H. 2005. Composición y abundancia de artrópodos epígeos del Parque Nacional Llanos de Challe: impactos del ENOS de 1997 y efectos del hábitat pedológico. Rev Chil de Hist Nat 78: 635-650.
• CLOUDSLEY-THOMPSON J. 2001. Thermal and water relations of desert beetles. Naturwissenschaften 88: 447-460.
• COLWELL R. 2000. EstimateS: Statistical Estimation of Species Richness and Shared Species from Samples, version 7.5.2. (http://viceroy.eeb.uconn.edu/estimates).
• CORCOS D ET AL. 2018. Predator and parasitoid insects along elevational gradients: role of temperature and habitat diversity. Oecologia 188: 193-202.
• CURRIE DJ. 1991. Energy and large-scale patterns of animal- and plant-species richness. The Am Nat 137: 27-49.
• CURRIE DJ ET AL. 2004. Predictions and tests of climate-based hypotheses of broad-scale variation in taxonomic richness. Ecol Lett 7: 1121-1134.
• DE LAFONTAINE G & HOULE G. 2007. Species richness along a production gradient: a multivariate approach. Am J Bot 94: 79-88.
• DE LOS SANTOS A, DE NICOLÁS JP & FERRER F. 2002. Habitat selection and assemblage structure of darkling beetles (Col. Tenebrionidae) along environmental gradients on the Island of Tenerife (Canary Islands). J A Environ 52: 63-85.
• EVANS KL, GREENWOOD JJD & GASTON KJ. 2005. Dissecting the species–energy relationship. P Roy Soc B-Biol Sci 272: 2155-2163.
• FATTORINI S. 2006. Detecting biodiversity hotspots by species–area relationships: a case study of Mediterranean beetles. Conserv Biol 20: 1169-1180.
• FLORES GE. 1999. Systematic revision and cladistic analysis of the Neotropical genera Mitragenius Solier, Auladera Solier and Patagonogenius gen. n. (Coleoptera: Tenebrionidae). Entomol Scand 30: 361-396.
• FLORES GE. 2004. Systematic revision and cladistic analysis of the Patagonian genus Platesthes (Coleoptera: Tenebrionidae). Eur J Entomol 101: 591-608.
• FLORES GE & CHANI-POSSE M. 2005. Patagonopraocis, a new genus of Praocini from Patagonia (Coleoptera: Tenebrionidae). Ann Zool 55: 575-581.
• GASTON KJ & CHOWN SL. 1999. Elevation and climatic tolerance: a test using dung beetles. Oikos 86: 584-590.
• GEBERT F, STEFFAN-DEWENTER I, MORETTO P & PETERS MK. 2020. Climate rather than dung resources predict dung beetle abundance and diversity along elevational and land use gradients on Mt. Kilimanjaro. J Biogeogr 47: 371-381.
• GONZÁLEZ-MEGÍAS A, GÓMEZ JM & SÁNCHEZ-PIÑERO F. 2008. Factors determining beetle richness and composition along an altitudinal gradient in the high mountains of the Sierra Nevada National Park (Spain). Ecoscience 15: 429-441.
• GOTELLI NJ & COLWELL RK. 2001. Quantifying biodiversity: procedures and pitfalls in the measurement and comparison of species richness. Ecol Lett 4: 379-391.
• HAUMMAN L. 1947. La vegetación de la Argentina. In: Sociedad Argentina de estudios geográficos (GAEA) (ed) Geografía de la República Argentina. Casa Coni, Argentina, p. 5-346.
• HAWKINS BA ET AL. 2003. Energy, water, and broad-scale geographic patterns of species richness. Ecology 84: 3105-3117.
• HE F, LEGENDRE P & LAFRANKIE J. 1996. Spatial pattern of diversity in a tropical rain forest in Malaysia. J Biogeogr 23: 57-74.
• HIJMANS RJ, CAMERON SE, PARRA JL, JONES PG & JARVIS A. 2005. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol 25: 1965-1978.
• JANZEN DH. 1967. Why mountain passes are higher in the tropics. Am Nat 101: 233-249.
• KULZER H. 1955. Monographie der Scotobiini (Zehnter Beitrag zur Kenntnis der Tenebrioniden). Entomologische Arbeiten aus dem Museum Georg Frey 6: 383-478.
• KULZER H. 1958. Monographie der südamerikanischen Tribus Praocini (Col.) (16 Beitrag zur Kenntnis der Tenebrioniden). Entomologische Arbeiten aus dem Museum Georg Frey 9: 1-105.
• KULZER H. 1962. Neue Tenebrioniden sus Südamerika (Col.) (23 Beitrag sur Kenntnis der Tbnebrioniden). Entomologische Arbeiten aus dem Museum Georg Frey 13: 79-100.
• KULZER H. 1963. Revision der südamerikanischen Gattung Nyctelia Latr. (Col. Teneb.) (24 Beitrag zur Kenntnis der Tenebrioniden). Entomologische Arbeiten aus dem Museum George Frey 14: 1-71.
• LEGENDRE P. 1993. Spatial autocorrelation: trouble or new paradigm? Ecology 74: 1659-1673.
• LEGENDRE P & LEGENDRE L. 1998. Numerical ecology. 2nd English Edition, Elsevier, Amsterdam.
• MARTÍNEZ CARRETERO E, BORGUI CE & DALMASSO AD. 1999. Evaluación de impacto ambiental. Propuesta metodológica expeditiva y estudio de caso en la reserva altoandina Laguna del Diamante, Argentina. Multequina 8: 111-120.
• MARTÍNEZ CARRETERO E & MÉNDEZ E. 1992. La vegetación de la vertiente oriental de la cordillera real, Mendoza-Argentina. Multequina 1: 99-106.
• MCCAIN CM. 2007. Could temperature and water availability drive elevational species richness patterns? A global case study for bats. Global Ecol Biogeogr 16: 1-13.
• MCINTYRE NE. 2000. Community structure of Eleodes beetles (Coleoptera: Tenebrionidae) in the shortgrass steppe: scale-dependent uses of heterogeneity. West N Am Naturalist 60: 1-15.
• MOSER D, DULLINGER S, ENGLISCH T, NIKLFELD H, PLUTZAR C, SAUBERER N, ZECHMEISTER HG & GRABHERR G. 2005. Environmental determinants of vascular plant species richness in the Austrian Alps. J Biogeogr 32: 1117-1127.
• MOURELLE C & EZCURRA E. 1996. Species richness of Argentine cacti: A test of biogeographic hypotheses. J Veg Sci 7: 667-680.
• NJOVU HK, STEFFAN-DEWENTER I, GEBERT F, SCHELLENBERGER COSTA D, KLEYER M, WAGNER T & PETERS MK. 2021. Plant traits madiate the effects of climate on phytophagous beetle diversity on Mt. Kilimanjaro. Ecology 102: e03521.
• O’BRIEN EM. 1993. Climatic gradients in woody plant species richness: Towards an explanation based on an analysis of Southern Africa’s woody flora. J Biogeogr 20: 181-198.
• O’BRIEN EM, WHITTAKER RJ & FIELD R. 1998. Climate and woody plant diversity in southern Africa: relationships at species, genus and family levels. Ecography 21: 495-509.
• PITHER J. 2003. Climate tolerance and interspecific variation in geographic range size. P Roy Soc B-Biol Sci 270: 475-481.
• QIAN H ET, FRIEDLEY JD & PALMER MW. 2007. The latitudinal gradient of species–area relationships for vascular plants of North America. Am Nat 170: 690-701.
• R DEVELOPMENT CORE TEAM. 2009. R: a language and environment for statistical computing, version 2.8.1. R Foundation for Statistical Computing, Vienna, Austria.
• RAHBEK C, GOTELLI NJ, COLWELL RK, ENTSMINGER GL, RANGEL TF & GRAVES GR. 2007. Predicting continental-scale patterns of bird species richness with spatially explicit models. P Roy Soc B-Biol Sci 274: 165-174.
• RANGEL TF, DINIZ-FILHO JA & BINI LM. 2010. SAM: a comprehensive application for spatial analysis in macroecology. Ecography 33: 1-5.
• ROHDE K. 1992. Latitudinal gradients in species diversity: the search for the primary cause. Oikos 65: 514-527.
• ROIG-JUÑENT SA, FLORES GE & MATTONI C. 2003. Consideraciones biogeográficas de la Precordillera (Argentina), con base en artrópodos epígeos. In: Morrone J & Llorente-Bousquets J (Eds) Una perspectiva latinoamericana de la Biogeografía. Las Prensas de Ciencias, México DF, p. 275-288.
• RUGGIERO A & KITZBERGER T. 2004. Environmental correlates of mammal species richness in South America: effects of spatial structure, taxonomy and geographic range. Ecography 27: 401-417.
• SATORRA A & BENTLER PM. 2001. A scaled difference chi-square test statistic for moment structure analysis. Psychometrika 66: 507-514.
• SCHEINER SM, MITCHELL RJ & CALLAHAN HS. 2000. Using path analysis to measure natural selection. J EvolutionBiol 13: 423-433.
• ŞENYÜZ Y, LOBO JM & DINDAR K. 2019. Altitudinal gradient in species richness and composition of dung beetles (Coleoptera: Scarabaeidae) in an eastern Euro-Mediterranean locality: Functional, seasonal and habitat influences. Eur J Entomol 116: 309-319.
• SHIPLEY B. 1997. Exploratory path analysis with applications in ecology and evolution. Am Nat 149: 1113-1138.
• SHIPLEY B. 2000. Cause and correlation in biology. Cambridge University Press, Cambridge.
• SRIVASTAVA DS & LAWTON JH. 1998. Why more productive sites have more species: an experimental test of theory using tree-hole communities. Am Nat 152: 510-529.
• STEVENS GC. 1989. The latitudinal gradient in geographical range: how so many species coexist in the tropics. Am Nat 133: 240-256.
• SUGGITT AJ, WILSON RJ, AUGUST TA, FOX R, ISAAC NJB, MACGREGOR N, MORECROFT MD & MACLEAN IMD. 2015. Microclimate effects landscape level persistence in the British Lepidoptera. J Insect Conserv 19: 237-253.
• TURNER RG, GATEHOUSE CM & COREY CA. 1987. Does solar energy control organic diversity? Butterflies, moths and the British climate. Oikos 48: 195-205.
• VAN RENSBURG BJ, CHOWN SL & GASTON KJ. 2002. Species richness, environmental correlates, and spatial scale: A test using South African birds. Am Nat 159: 566-577.
• VON HUMBOLDT A. 1808. Ansichten der Natur mit wissenschaftlichen. Erläuterungen, Stuttgart.
• WALTHER G-R, BERGER S & SYKES M. 2005. An ecological ‘footprint’ of climate change. P Roy Soc B-Biol Sci 272: 1427-1432.
• WARING RH, COOPS NC, FAN W & NIGHTINGALE JM. 2006. MODIS enhanced vegetation index predicts tree species richness across forested ecoregions in the contiguous U.S.A. Remote Sens of Environ 103: 218-226.
• WERENKRAUT V, FERGNANI PN & RUGGIERO A. 2015. Ants at the edge: a sharp forest-steppe boundary influences the taxonomic and functional organization of ant species assemblages along elevational gradients in northwestern Patagonia (Argentina). Biodivers Conserv 24: 287-308.
• WERENKRAUT V & RUGGIERO A. 2014. The richness and abundance of epigaeic moutain beetles in north-western Patagonia, Argentina: assessment of patterns and environmental correlates. J Biogeogr 41: 561-573.
• WRIGHT DH. 1983. Species–energy theory: an extension of species–area theory. Oikos 41: 496-506.
• YU X-D, LÜ L, LUO T-H & ZHOU H-Z. 2013. Elevational Gradient in Species Richness Pattern of Epigaeic Beetles and Underlying Mechanisms at East Slope of Balang Mountain in Southwestern China. PLoS ONE 8: e69177.
• YU X-D, LU L, WANG F-Y, LUO T-H, ZOU S-S, WANG C-B, SONG T-T & ZHOU H-Z. 2016. The relative importance of spatial and local environmental factors in determining beetle assemblages in the inner Mongolia grassland. PLoS ONE 11: e0154659.

Publication Dates

• Publication in this collection
  28 Aug 2023
• Date of issue
  2023

History

• Received
  23 Mar 2021
• Accepted
  28 Oct 2022