Environmental factors determining the diversity of darkling beetles (Coleoptera: Tenebrionidae) in arid, high-altitude environments in Northwestern Argentina

Abstract Factors influencing the diversity of tenebrionid beetles in arid, high-altitude environments in the northwest of Argentina are little-known. Using pitfall traps and suction sampling in 30 sites, we collected these beetles in Altos Andes (AA) and Puna (PU) ecoregions and evaluated how local and regional factors influenced their assemblages. During each sampling date, we registered variables related to climate, vegetation, and soil in each sampling site. In total, we recorded 270 individuals from 21 species, richness of species was higher in PU than in AA, but without a statistically significant difference. Twelve species were present exclusively in PU, while two in AA. Tenebrionid assemblages in both ecoregions had different dominant species, in PU, Psectrascelis cariosicollis while Epipedonota interandina in AA. Beta diversity between ecoregions was moderate and we recorded a high beta diversity and high turnover of species in each ecoregion. The geographical distance between the sites influenced the similarity of the assemblages inside and between ecoregions, therefore, an increase in the geographical distance between sites determined a decrease in the similarity of the assemblages. The increase in elevation of the sampled sites did not produce significant changes in the abundance and richness of Tenebrionidae. There is no individual factor responsible for the darkling beetle communities structuration in these high-altitude environments. Different combinations of both local (soil and vegetation) and regional (climate and geographical distance) environmental factors together explain the ecoregional tenebrionid assemblages. We suggest including them among the focus groups in conservation studies to evaluate anthropogenic activities in this area.


INTRODUCTION
Cold-adapted, high-altitude species are very vulnerable to climate breakdown and particularly sensitive to global warming (Wilson et al. 2007, Dieker et al. 2013).All mountainous regions in the 20th century experienced aboveaverage warming compared to the world average (UNESCO 2014).This affects frost lines, dew point elevation, and other environmental factors that have a significantly greater impact on the functioning of ecosystems at high elevation than in the lowlands (Herzog et al. 2012).Many mountain species react to this climate breakdown changing their ranges to higher elevations to find appropriate living conditions (Dieker et al. 2013).They would eventually be unable to find sufficient elevations that would allow survival and thus face a greater risk of extinction (Lovejoy 2010).In turn, this displacement causes biotic wear in the lowlands as species that normally live in this region to ascend along the elevation gradient (Colwell et al. 2008).For this reason, knowledge concerning lowland-species conservation has recently driven the study of life at high elevation (Altshuler & Dudley 2006).
The Puna (PU) and Altos Andes (AA) are arid South American ecoregions that are adjacent and are distinguished by their high elevation (Morello et al. 2012).In Argentina, they extend from the border with Bolivia in the extreme northwest of the country to the South until the north of Neuquén province (AA, Reboratti 2005) and to the north of Mendoza province (PU, Martínez Carretero 1995).They are included in the Global 200 as "Vulnerable" and are considered priority conservation areas (Olson & Dinerstein 2002) because its fragile desert environments contain a great wealth of endemic species threatened mainly by the increase in temperature as a consequence of global warming, the increase in mining activity and ecotourism (Morello et al. 2012).In these arid ecosystems, arthropods play important roles as decomposers, herbivores, granivores and predators, controlling the cycle of nutrients and energy that flows through the different levels of the trophic chain (Polis 1991, Ayal 2007).However, until recently, knowledge about arthropod biodiversity in the Argentine PU and AA was incipient.In recent years, some authors studied arthropod diversity in the Salta province PU (Cava et al. 2012, 2013, González-Reyes et al. 2012, 2017, Cruz et al. 2017), but there were no similar studies for AA.
Among arthropods, darkling beetles (Coleoptera: Tenebrionidae), especially those from the subfamily Pimeliinae, represent dominant faunal elements in biomass and number of individuals in desert and semidesert environments (Watt 1974, Doyen 1994); they present morphological, physiological and ethological adaptations for life in these environments (Flores 1998, Cloudsley-Thompson 2001, Matthews et al. 2010).Their distribution and abundance are affected by factors such as vegetation type and soil topography (Ayal & Merkl 1994, Krasnov & Ayal 1995, Botes et al. 2007).Darkling beetles play a key role in the fragmentation of vegetal and animal matters, nutrient cycling and trophic networks, mainly as components of the diet of other consuming organisms (Flores 1998, Flores & Debandi 2004, Aballay et al. 2016).Their assemblages generally show spatial and temporal changes concerning environmental condition variations (Cartagena 2002, Botes et al. 2007, Lescano et al. 2017).Therefore, they are considered a good indicator group of the state of environments and climate breakdown (Cartagena 2002, Michaels 2007).
Spatial heterogeneity and seasonality (Polis 1991) are considered determinants of species diversity in arid environments.In different arid and semi-arid ecosystems around the world, heterogeneity of vegetation, soil type and climate has been documented to affect abundance and composition patterns of the darkling beetle communities (de los Santos et al. 2000, Sackmann & Flores 2009, Lescano et al. 2017).On the other hand, some studies carried out in Mongolia (Pfeiffer & Bayannasan 2012, Yu et al. 2016, Khurelpurev & Pfeiffer 2017) detected an influence of geographical distance on the structuring of mountain darkling beetle communities.At present, in Argentina, relatively little is known about the biodiversity of darkling beetles that inhabits the PU and AA, and nothing is known about geographical distance and environmental factors that could affect their diversity patterns.
The purpose of this study was to analyze the spatial variations of darkling beetles assemblages in PU and AA in Salta province and to identify the factors that determine them.Therefore: (1) we compared the darkling beetle alpha and beta diversities between ecoregions as well as their variation in terms of abundance and species richness; (2) we evaluated the relative contribution of alpha and beta diversities to regional diversity (gamma); (3) we identified the possible abiotic (elevation, climatic and soil heterogeneity) and biotic (vegetation heterogeneity) factors that can influence darkling beetle diversity patterns; and (4) we assessed whether geographic distance explains the community dissimilarities.We specifically test the following hypotheses: a) PU and AA darkling beetle communities differ in richness, abundance and species composition, b) Existence of a high beta diversity of darkling beetles between sites in the same ecoregion and between ecoregions, and c) Geographical distance and environmental factors (elevation, climate, soil, and vegetation) structure the patterns of diversity of darkling beetles from PU and AA.As each ecoregion is a unit of biodiversity on a regional scale (Dinerstein et al. 2000), with its own and homogeneous environmental characteristics (Morello et al. 2012) and containing different species communities (Olson et al. 2001), it is expected that there are particular assemblages of darkling beetles in each ecoregion.On the other hand, as the PU and AA ecoregions are characterized by presenting several ecosystem complexes with different plant communities in each of them (Morello et al 2012), which makes them heterogeneous environments; they would produce a marked turnover of darkling beetle species.Finally, as darkling beetle communities could be influenced by elevation, vegetation, soil, and environmental changes (Botes et al. 2007, Sackmann & Flores 2009, Lescano et al. 2017), it is predicted that changes in climate and vegetation with the altitude in short geographical distances may also be the main factors in determining the diversity patterns of these communities in these high-altitude environments.

Study area
The study area was located between 24-25°S and 66-67°W (Table SI -Supplementary Material), which included the departments Cachi, La Poma and Los Andes in Western Salta province in Argentina (Fig. 1).This area is well represented by the Puna (PU) and Altos Andes (AA) ecoregions.The samples included protected natural environments, such as the Abra del Acay Natural Monument and the Reserva Natural de Uso Múltiple Los Andes (one of the largest reserves in Argentina).Additionally, were included unprotected natural environments, such as the Los Cardones-La Poma and Sierra del Cobre Longitudinal Corridor, which are of great interest for conservation and deserve protection (Chebez 2005).
The Puna (PU) is an extensive plateau located in west-central South America that occupies parts of Argentina, Bolivia, Chile, and Peru (Cabrera & Willink 1980).It is an undulating plain crossed by mountain ranges located above 3,000 m.a.s.l in extreme northwestern Argentina.The climate is dry, windy, and cold, with a seasonal thermal range between 2-16°C and marked daily temperature fluctuations (approximately 30°C).Most annual rains occur during the summer and vary between 100-800 mm.Like the AA ecoregion, the PU has the highest proportion of rocky outcrops with no soil, and its edaphic characteristics explain the low productivity level of the land (Morello et al. 2012).The vegetation is steppe and consists mainly of low tola bushes (Paraestrephia sp.), añagua (Adesmia horridiscula), and yareta (Azorella yareta).Grasses of the genus Distichlis appear only occasionally in places where accumulated water in spring creates differentiated microenvironments called "vegas".The arbustales and tolillares (Fabiana sp.) dominate in drier PU areas; rich-ricales (Acantholippia salsoloides) and fasciculate-type grasslands occur in less dry areas (Morello et al. 2012).Vegetation coverage is highly variable and depends on climatic factors and relief, especially water accumulation and permanence in the soil.Only in some places, plant communities cover less than 10% of the soil.
The Altos Andes (AA) occupies high mountain areas with deep valleys in west Argentina along the Andes mountains located above 4,000 m.a.s.l although towards the south of the country this limit drops to 2,500 m (Reboratti 2005).The climate is cold and dry with very strong winds, average monthly temperatures are below 0°C for more than half of the year.Annual rainfall ranges between 100-200 mm.Vegetation consists mainly of low and woody shrubs of the genera Azorella, Pycnophyllum, Nototriche, Werneria, and Xenophyllum, in addition to cushion-shaped plants and some grasses well adapted to local conditions, such as Festuca orthophylla, Festuca chrysophylla, Poa gymnantha and Stipa sp.Low temperature and precipitation prevent the formation of dense vegetal cover and complete soils.

Sampling
The samplings were conducted in December 2014 and June 2015 because these are the only times of the year where the roads are drivable in both ecoregions.Thirty geo-referenced sites (12 AA and 18 PU) were considered, separated from each other by at least 4 km.Each site covered an area of approximately 2 ha and in all of them, the composition of the surrounding landscape (600 meters around the collected area) was homogeneous.Darkling beetles were collected using complementary sampling techniques as suggested by Halffter et al. (2001): we used pitfall traps to capture epigean darkling beetles and a G-Vac suction device (Garden-Vaccum) to capture species that live above the vegetation.During each sampling date and in each site, five darkling beetle samples using the suction method and five samples with pitfall traps were taken simultaneously, making up a total of 600 samples.Suction samples were taken at random with a Stihl aspirator with a 1.10 m long tube with a 12 cm diameter.Each suction sample was defined as the suction of the vegetation in an area of one square meter for one minute.At each sampling point a linear transect was drawn at random on which the pitfall traps were placed with a separation of 20 meters.

Local variables
During each sampling date, a set of variables related to vegetation and soil at each site was taken in the field to characterise the heterogeneity of the sampled habitats.To analyse the plant communities at each site, we measured species richness and variables related to horizontal and vertical vegetation structures.For horizontal structure, the vegetation coverage percentage was considered using five quadrants 5 x 5 meters taken at random.For vertical structure, we evaluated vertical stratification (number of individuals per shrub and herbaceous stratum) and plant architecture: density of individuals with central stems without branching (cardons), with ramification from the base (bushes), and grass density.Thus, the percentage of vertical stratum was analysed with the Vertical Vegetational Structure Analysis (VESTA) photographic method (Zehm et al. 2003).For this last method, at each site, three random points were selected where eight photographs were taken (two for each cardinal point up to 2 m high); we placed a 1 x 1 meter white panel behind the vegetation for contrast.Subsequently, using Adobe Photoshop CS6, we quantified the percentage of each stratum.Productivity for each study site was calculated using the Normalised Vegetation Index (NDVI) because it allows one to determine vegetation distribution and dynamics (Pettorelli et al. 2005).This index was calculated from the combination of spectral bands (NDVI = (NIR-RED)/(NIR + RED)) of images from the OLIS sensor of the Landsat satellite obtained from the website https://earthexplorer.usgs.gov/with Path = 232 and Row = 77 for the study dates.The images were processed with the program QGis 2.18-Las Palmas (QGIS Development Team 2016).
To analyse soil heterogeneity, we considered five squares 0.5 x 0.5 m randomly distributed over two ha.They were photographed perpendicularly and later analysed with Adobe Photoshop CS6 following the method proposed by Gilbert & Butt (2009) to obtain the percentage of sand, rock, gravel, soil cover and percentage of dead organic matter.A differential colour was assigned for each category, and the values for each colour were obtained as the percentage of the total number of pixels in the photograph.To perform the analysis we used average values of different soil covers per site.

Regional variables
Due to the influence of climate in the selected ecoregions, for each site, we considered 19 climatic variables related to temperature and precipitation, obtained from Worldclim (2018) USGS-WIST (NASA; http://www.worldclim.org) with a resolution of 30 sec (approximately 90 m).Additionally, soil salinity was considered as a variable because the excess salt in arid-and semi-arid-zone soil is a global phenomenon (Asfaw et al. 2016, Elhag 2016).The Salinity Index (SI) was calculated for each site (Dehni & Lounis 2012), which was measured from the combination of spectral bands of images from Landsat satellite, Path = 232 and Row = 77 for the study dates, with corresponding radiometric corrections.The images were processed with QGis 2.18-Las Palmas (QGIS Development Team 2016).

Data analysis
Several sites did not contain individuals and, for this reason, these were not considered in the analyses that required the presence of species.Originally 30 sites were sampled but only in 22 of them we recorded darkling beetles (Fig. 1).
Normality tests (goodness-of-fit test) and homogeneity of variances (Bartlett) were performed on all obtained data groups.No data group met the requirements, therefore nonparametric statistical tests were applied.

Inventory and alpha diversity
To determine if the sample obtained is representative of the measured attribute, we evaluated the quality of the inventories obtained in each ecoregion using non-parametric estimators.These use the information on the number of rare or infrequent species in the inventory to estimate the number of undetected species (Chao et al. 2009).We calculated the integrity of inventories, such as the relationship between observed species richness and true diversity (Jost 2006(Jost , 2007) ) estimated by ACE (abundance-based coverage estimator), obtained using the SPADE software (Chao & Shen 2010).The ACE estimator was chosen following the proposal of Chao & Shen (2010), who recommend its use when the communities are highly heterogeneous, such as in this study.The completeness values of the inventories obtained for each ecoregion were evaluated following the criteria of Cardoso (2009), which establishes three categories: "reasonable" (50-70% completeness), "complete" (70-80%), and "exhaustive" (90-100%).
To compare species richness among ecoregions and sites, we used extrapolation analysis of rarefaction of the Hill numbers for q = 0 (species richness), q = 1 (Shannon diversity) and q = 2 (Simpson diversity) (see calculation formulas in Jost 2006), based on individuals with the same level of sample coverage, with 95% confidence intervals and 100 permutations, using the iNEXT program (Hsieh et al. 2013).This method ensures that samples are compared with equal integrity, regardless of sample size, a feature that allows for more robust inferences about community species richness patterns (Chao & Jost 2012).The proportion of estimated richness detected in each site, as well as the deficit of coverage of the sample (1-coverage of the sample) (Chao & Jost 2012) were considered to evaluate the representativeness of the samples obtained on each site.Additionally, to determine and compare the structure of the darkling beetle communities registered in different ecoregions, we used abundance range curves (Whittaker curves), where the range was considered as species listed in decreasing abundance order.

Effect of geographical distance and environmental variables on darkling beetle beta diversity
The contribution of each species to the dissimilarity observed between the samples was evaluated with a SIMPER analysis using the Bray-Curtis similarity measure with the PAST 3.16 program (Hammer et al. 2001).The relative importance of geographical distance in the obtained beta diversity patterns was verified by a simple Mantel test that considered a faunistic matrix (generated using Jaccard distance data as a measure of the similarity between darkling beetle communities between sites) and a geographic distance matrix (generated using site coordinates in WGS84 format and geographic distance as a measure of the distance between pairs of sites).The Mantel test was performed at a regional level and in each ecoregion with the PAST 3.16 program (Hammer et al. 2001) using the Pearson correlation coefficient and 10,000 permutations to evaluate the statistical significance between matrices.
To evaluate beta diversity between ecoregions and among multiple-sites of the same ecoregion, we used the approach of Baselga (2010a) that measures total beta diversity using the Sørensen dissimilarity index (β SOR ).Because beta diversity (β SOR ) can reflect two underlying phenomena: turnover (β SIM ) and nestedness (β SNE ) and these are generated by different processes (replacement of species and loss of species, respectively) (Baselga 2010a), the relative contribution of these components was quantified and their values were expressed as percentages.This analysis was performed using the Betapart package in R (Baselga et al. 2013), where β Sørensen = β Simpson + β Nestedness-resultant .We follow the same notation as Baselga (2010a) and Baselga et al. (2013), using upper-case letters for the multiple-site metrics (β SOR , β SIM , β SNE ), and lower-case letters for pairwise dissimilarities (β sor , β sim , β sne ).
To test whether there are significant differences in the abundance and richness of species depending on the elevation of the sampled sites we used the Kruskal-Wallis test.This analysis was carried out with the PAST 3.16 program (Hammer et al. 2001) with a significance level of 95% and 21 degrees of freedom.
We used generalized linear models (GLM) to evaluate the influence of geographical distance and environmental variables (climate, soil, and vegetation) on the darkling beetles' distribution pattern at a regional scale and by each ecoregion.Without considering singletons and doubletons, species richness and abundance of darkling beetles was used as the response variable, while the environmental factors (climate, soil, vegetation) and geographic distance as explanatory variables.Spearman's rank correlation was used to identify collinearity between the independent variables, as it is a non-parametric measure of statistical dependence (Zar 1999).The variables were excluded when the coefficient r> 0.7.The variance inflation factor (VIF) was evaluated for any remaining collinearity in the full models and the variables excluded with VIF> 5 since they indicate collinearity between predictors (Heiberger & Holland 2004).For each model, a backward elimination procedure was performed to debug insignificant variables without losing relevant information (significance level p <0.05), and in this way to be able to obtain the adequate minimum model.Finally, a pseudo-R2 was calculated (Zuur et al. 2009), with deviation values for the best models: (null deviance -residual deviance) / null deviance * 100.These statistical analyses were performed using R (R Core Team 2019).Before the GLM analyses, all the variables considered in this study (vegetation, soil, climate, and geographical distance) were subjected to a principal components analysis (PCA) with the PAST 3.16 Program (Hammer et al. 2001) to evaluate and overcome the problems of correlations between variables.In this way, the variables that were correlated and according to their influence on the first two axes of PCA were eliminated.This analysis has been widely used as a tool in exploratory analyses (Zarechahouki 2011).The selected variables were relativized through the standard deviation (mean = 0, variance = 1) because it is a very useful transformation for environmental variables (McCune et al. 2002).
Importance of beta diversity at the regional scale Species diversity distribution across the landscape, and the contribution of each level to regional diversity (gamma diversity), was analysed and compared by a multiplicative partition of the species richness with the PARTITION 3.0 program (Veech & Crist 2009).Different authors (Jost et al. 2010, Baselga 2010b, Tuomisto 2010) recommend their use to evaluate the relative differentiation between communities because under this method beta (β) is independent of alpha (α) and gamma (γ), a feature that allows one to study the contributions of each component to gamma diversity without any type of distortion and within a rigorous and consistent mathematical framework (Baselga 2010b, Jost 2010).Multiplicative beta diversity measures how much more diverse the region is compared with the average of its zones (Pereyra & Moreno 2013).Thus, we considered three hierarchical levels (samples, sites and ecoregions), where γ = α1 (within the sample) x β1 (between samples) x β2 (between sites) x β3 (between ecoregions), considering that sample is equal to microhabitats.These observed α and β values were compared with the expected values in a distribution of the random samples (Crist et al. 2003) using the same program.If p = 1 or close to this value, the species are randomly distributed at the considered level, and therefore, the samples are independent.Conversely, if p = 0 or close, the species are not randomly distributed at that level.

Inventory and alpha diversity
A total of 270 darkling beetle specimens were collected belonging to 21 species, 13 genera, 8 tribes, and 2 subfamilies (Table SII).Of the 21 registered species, four were completely novel species: Emmallodera sp., Physogaster sp., Physogasterini gen.nov, sp.nov., and Praocis (Postpraocis) sp.Nineteen darkling beetle species were reported in the PU and nine in the AA.The ecoregions studied shared seven species, while twelve species were present exclusively in PU sites and two were exclusive to AA.
The inventory achieved in PU detected 62.5% of the species estimated by ACE (30.4 species) and 78.6% in AA (11.4 species).These values indicate that the inventory obtained for PU was "reasonable", while for AA it was "complete".The sample coverage reached 96% in both ecoregions denoting an almost complete record of the species present.There were no significant differences between them, and thus we inferred that the samples were comparable.The interpolation (rarefaction)/extrapolation curves revealed that, when comparing darkling beetle species richness at the same level of sampling coverage, there were no statistically significant differences between the PU and AA assemblages (Fig. 2a) (bootstrap with significance level of 95%, iNEXT program).The average number of individuals collected per site was low (~ 12 individuals per site).However, it was observed that in 20 sampled sites (91%) it was possible to detect between 70-100% of the estimated species, while in the remaining sites (PU1 and PU8) only half of the expected species were detected (Fig. 2b).The deficit of coverage of the sample was equal to zero in 14 sampled sites (64%) (Table SII), that is to say, that in them there are no species to detect.In six sites (27%), the sample coverage deficit was low (between 0.1 and 0.3) and, conversely, only in two sites was its value greater than 0.6 (Table SII).Therefore, the sample size obtained provides representative data of the studied sites.
The Whittaker curves (Fig. 3) showed that the PU and AA darkling beetle communities were structured differently; each had some endemic species of high-elevation environments.The PU presented a greater number of species with intermediate and rare abundance, with Psectrascelis cariosicollis Fairmaire (42.5%) being markedly dominant, followed by Hylithus forsteri meridionalis Kaszab (13.5%) and Praocis magnoi Molinari (10.4%).In contrast, the AA community presented lower species richness and, it was dominated by Epipedonota interandina Vidal and Flores (40.3%), followed by H. forsteri meridionalis (33.8%).
Beta diversity: the effect of geographic distance and environmental variables on darkling beetle beta diversity SIMPER analysis demonstrated that 80% of the dissimilarity between the PU and AA communities was explained by 7 species.Hylithus forsteri meridionalis (24.64%) contributed most to the difference between the ecoregional communities, followed by Epipedonota interandina (17.78%) and Psectrascelis cariosicollis (11.39%).
Mantel tests carried out regionally and for each ecoregion showed that the geographical distance between the sampled sites influenced the similarity of the dark beetle assemblages, a result that confirms that the communities were more different as the distance between them increased (regional level R = 0.3064, p = 0.001; PU R = 0.3832, p = 0.0016; AA R = 0.5723, p = 0.0005).Beta diversity between ecoregions (β sor = 0.5) was moderate, with 60% of its value due to nestedness and the other 40% due to species turnover (β sne = 0.3, β sim = 0.2).In contrast, the beta diversity among the multiple sites within each ecoregion was high and very similar for the PU (β SOR = 0.95) and AA (β SOR = 0.94), with spatial species turnover the main factor that was responsible for most of the ecoregional beta diversity (PU: β SIM = 0.88, β SNE = 0.06; AA: β SIM = 0.79, β SNE = 0.15).
The best model at the regional scale (GLM) explained 79.68% (pseudo-R 2 , Table I) of the total variance in the abundance of darkling beetles.On a regional scale, the variables that had a significant and positive relationship with the abundance of these beetles were the density of cardons (Z = 2.60, p = 0.009), vertical stratification of the vegetation (Z = 6.20, p = 5.35e-10), percentage of organic matter in the soil (Z = 10.19,p = <2e-16) and geographical distance (Z = 4.08, p = 4.45e-05).In addition, the annual average temperature (Z = -5.58,p = 2.31e-08) and the annual precipitation (Z = -12.85,p = <2e-16) had a significant and negative relationship with abundance (Table I).The best regional scale model for species richness (GLM) explained only 37.7% (pseudo-R 2 ) of the total variance.Soil salinity (Z = -2.55,p = 0.01) showed a significant and negative relationship with species richness, while geographic distance (Z = 1.99, p = 0.04) was positively related (Table I).
The best abundance model for PU (GLM) explained 92.36% (pseudo-R 2 ) of the variance, the variables that showed a significant and positive relationship with the abundance of darkling beetles in this ecoregion were: vertical stratification of vegetation (Z = 2.06, p = 0.03), percentage of rock in the soil (Z = 2.01, p = 0.04), and annual average temperature (Z = 7.12, p = 1.03e-12).The percentage of gravel in the soil (Z = -6.86,p = 6.72e-12), the rainfall in the warmest quarter of the year (Z = -7.01,p = 2.23e-12), and the geographical distance (Z = -4.68,p = 2.85e-06) showed a significant and negative relationship with abundance (Table I).In all the PU darkling beetles species richness models constructed, no variable considered in this study had a significant effect on richness.
The best abundance model for AA (GLM) explained 66.35% (pseudo-R 2 ) of the variance, the variables with significant and negative relationship with the abundance of darkling beetles in AA were the percentage of sand in the soil (Z = -3.05,p = 0.002) and geographical distance (Z = -4.13,p = 3.52e-05), while the annual average temperature (Z = 4.44, p = 8.64e-06) showed a positive relationship (Table I).Finally, the best model of species richness for AA (GLM) explained 62.11% (pseudo-R 2 ) of the variance and a single variable, the geographical distance (Z = -2.04,p = 0.04) showed a significant and negative relationship with the richness of AA (Table I).

Importance of beta diversity at the regional scale
Multiplicative partitioning of gamma diversity revealed that beta diversity among sites and between ecoregions was significantly higher than expected by chance (p = 0).Conversely, at the sample level, alpha and beta diversity did not differ significantly (p = 1) indicating that species are randomly distributed within and between samples.In other words, all species had the same probability of being present in the samples.At the landscape scale, the greatest contribution was represented by beta diversity between sites, representing 55.5% of the total regional darkling beetle diversity (gamma), followed by beta diversity among samples that contributed 17.5%.Beta diversity between ecoregions represented 14.3% and finally, the alpha diversity of the samples provided 12.6%.

DISCUSSION
In this study we demonstrated that the abundance and species composition of Tenebrionidae was different between ecoregions (PU and AA) but we did not find significant differences in species richness.Consistent with our second prediction, the heterogeneity of the PU and AA environments produced high beta diversity and marked species turnover between sites.The variation in elevation of the sites in this study did not generate significant changes in the abundance and species richness of darkling beetles.On a regional scale, the abundance and species richness patterns of these beetles were determined by the density cardones, the vertical stratification of the vegetation, Table I. Results of generalized linear models calculated on a regional scale and for each ecoregion, showing the variables that influenced the abundance and species richness of darkling beetles from Puna and Altos Andes.

pseudo R²
Predictor variables Z value p-value percentage of organic matter in the soil, soil salinity, annual average temperature, annual precipitation, and geographical distance.None of the environmental variables considered in this study showed a significant relationship with the species richness of PU.However, more than 90% of the pattern observed in the abundance of PU was determined by a set of environmental variables such as: vertical stratification of vegetation, percentage of rock in the soil, percentage of gravel in the soil, annual average temperature, rainfall in the warmest quarter of the year and the geographical distance.The abundance and species richness patterns of darkling beetles in AA were determined by three environment variables: the percentage of sand in the soil, the geographical distance and the annual average temperature.This is the first study to analyse darkling beetle diversity patterns in Argentine highelevation ecoregions with different analysis scales (habitat and landscape).The results indicate that these ecoregions in Salta province harbour several darkling beetle species (S= 21) and some of them were cited as endemic to high-altitude environments (Flores 2000, Roig-Juñent et al. 2003, Flores & Pizarro-Araya 2004); of these 21, four are novel species that will be described in future works.This data shows that the Reserva Natural de Uso Múltiple Los Andes and the western Los Cardones-La Poma Longitudinal Corridor are important protected natural areas to support darkling beetle species richness.The species richness recorded in this study surpassed what is known for the Chilean transitional coastal desert (Alfaro et al. 2016) and western Abu Dhabi, United Arab Emirates (Saji & Al Dhaheri 2011), although it is lower than that reported for Monte (Roig-Juñent et al. 2001, Flores et al. 2004) and the Argentine Patagonian steppe (Sackmann & Flores 2009).
These darkling beetle species belong to subfamilies Pimeliinae (90%) and Tenebrioninae (10%; Table SII); members of Pimelinae exhibit good adaptations to arid conditions, with a generally restricted geographical distribution due to the fact that they are usually wingless and have undergone adaptive changes for ambulatory life (Flores 1998, Matthews et al. 2010).In Argentina, Pimeliinae species inhabit xeric environments of the PU, Prepuna, AA, Monte and south of Espinal and Patagonia (Cabrera & Willink 1980), and they are most representative of tenebrionids in these biogeographic provinces (Flores 1998).The two species of Tenebrioninae belong to Scotobiini, an endemic tribe to South America abundant in arid and semiarid environments (Doyen 1994, Matthews et al. 2010) with similar adaptations to living in such habitats (Silvestro et al. 2015).These beetles lay eggs in the soil and their larval development is hypogeous, meaning that practically all their biological cycle takes place above and below the ground (Flores 1998, Matthews et al. 2010, Silvestro & Michat 2016).
In recent years, mining activity has increased in the PU and AA, mainly for the extraction of lithium in areas close to the sites of this study.Lithium extraction generates important environmental changes because it involves consumption and pollution of water in high-altitude arid zones, construction of pools of flowering and decantation in sensitive ecosystems, for which the excavations remove tons of rocks and generate solid and chemical waste (Aguilar & Zeller 2012).All these activities alter the soils and could affect darkling beetles as they are very vulnerable to environmental changes.For this reason, these beetles could be used as a focus group to monitor the effects of mining activity in future studies of the region.
Of the twelve species present exclusively in PU in this study Omopheres difficilis, O. laevicollis, Achanius bicolor, A. rufonitens and other members of the genera Epipedonota Solier, Physogaster Guérin-Méneville, Platyholmus Solier, and Emmallodera Blanchard inhabit also the ecoregion of Monte (Flores & Vidal 2001, Silvestro & Flores 2012, IADIZA, collection data), indicating a biogeographic relationship with lands lower than PU (Roig-Juñent et al. 2003).The two species present exclusively in AA are distributed also in Chilean high elevation places (AA): Antofagapraocis subnudus Flores (Flores 2000) and Caenocrypticoides translucidus Kaszab (Flores & Pizarro-Araya 2004) as well as in other sites of Argentinean and Bolivian AA (IADIZA, collection data), showing a distribution within AA but not in other ecoregions.
Although species richness in both ecoregions was not markedly different, the PU and AA darkling beetle communities had different dominant species.In the PU environments, Psectrascelis cariosicollis dominated, while Epipedonota interandina dominated in AA.Species of the genus Psectrascelis Solier are nocturnal; they hide during the day under creeping plants, rocks or branches and are preyed upon mainly by lizards, birds, foxes and small mammals (Peña 1985).Epipedonota species are diurnal; they hide during the night or hottest hours under rocks, shrubs, and vegetal debris (Flores & Vidal 2001).Both species belong to Nycteliini, the most diverse darkling beetle tribe in South America (Flores 1997) and are the largest tenebrionids among the captured species (Table SII).The predominance of species of big size in our study is likely due to the fact that this morphological adaptation to the arid conditions allowed them to increase control of their homeostasis and longevity (Carrara et al. 2011).
The beta diversity between PU and AA was moderate and, when dividing it into its two components, it is clear that both species turnover (40%) and nestedness (60%) contribute markedly and in a similar way (see Baselga 2010a for interpretation details).In nature, species turnover may reflect species sorting by the environment or dispersal processes, whereas nestedness related to ordered extinctioncolonization dynamics (Si et al. 2016).Nestedness between sites occurs due to changes in species richness -sites with fewer species are strict subsets of richer sites (Baselga 2010a) and reflects an orderly loss of darkling beetle species to AA.This loss of species occurs because in AA there is a greater number of local and regional extinctions of Tenebrionidae species related to their lower productivity, historical factors and harsh environmental conditions (Baselga 2010a).In contrast, beta diversity in each ecoregion was high and caused almost entirely by species turnover between sites.This contrast in the results is explained because beta diversity and turnover are affected by changes at the scale (Soininen et al. 2018, Antão et al. 2019).For organisms with limited dispersal (such as darkling beetles), if the sampled area is small, more different will be the sites compared (Nekola & White 1999, Mac Nally et al. 2004, Barton et al. 2013).If the scale level increases, more species are shared between the sampled areas and microenvironmental differences are attenuated (Antão et al. 2019).Therefore, in searching to estimate general heterogeneity in communities, multiple-site metrics can be more reliable than average pair-wise metrics (Baselga 2013).
Other studies on insects performed in different mountainous environments coincide with our results because they also recorded patterns of high beta diversity and high turnover in ants (Liu et al. 2018, Castro et al. 2020), butterfly (Beirão et al. 2020), dung beetles (Moctezuma et al. 2018), bees andwasps (Perillo et al. 2017).Also, our findings have important implications that must be taken into account when making conservation decisions for the region.The high turnover within each ecoregion indicates that it would be necessary to protect multiple areas to ensure the conservation of the species in the region (Calderón-Patrón et al. 2012, Moctezuma et al. 2018).Due to the consequences of global warming, the two ecoregions are valuable for conservation regarding Tenebrionidae.We encourage AA conservation, despite having less alpha diversity, as this would preserve higher elevation habitats where PU species could move, in search of better conditions for survival.
Geographic distance influences the PU and AA darkling beetle communities, specifically by decreasing their assemblage similarity as the distance between the sites increases.Pfeiffer & Bayannasan (2012), Yu et al. (2016) and, Khurelpurev & Pfeiffer (2017) obtained similar results when studying darkling beetles along an ecological gradient in arid Mongolia.This pattern of diminishing similarity with geographical distance could be controlled jointly by processes based on niche, spatial configuration and neutrality (Tuomisto et al. 2003, Cottenie 2005, Soininen et al. 2007).For the group under study, the structural and environmental properties of the mountainous landscapes would represent a limitation to the dispersion of these beetles, which are flightless and predominantly endemic.Thus, there is greater species replacement between sites.
Our results show that geographical distance is a strong driver of the structuring of darkling beetle communities at high elevation.Climate affects these mountain beetle communities, but also there are other important local factors to consider, such as vegetation and soil heterogeneity.This fact is consistent with the idea that, in arid and semi-arid systems, both the abundance and composition of darkling beetle communities may be limited and influenced by vegetation and soil type (de los Santos et al. 2000, Fattorini 2009, Pan et al. 2015, Alfaro et al. 2016, Lescano et al. 2017).
In studied ecoregions, the vegetation is discontinuous in space, with different dominant plant species in different communities (Morello et al. 2012).This phenomenon is due to soil heterogeneity, as well as marked seasonal and climatic conditions (marked rainfall seasonality, large diurnal and seasonal thermal amplitudes and strong daily climatic fluctuations, among others) that could affect the presence of the species.These biotic and abiotic factors exercise control over the activities of epigean species subject to severe environmental conditions (Heatwole 1996, Whitford 2002, Lescano et al. 2017).The habitat preference of darkling beetles is related to edaphic factors.Tenebrionid beetles prefer sandy soils for laying of eggs because most of them have soil-inhabiting larvae which feed on all available organic matter including plant roots and seedlings (Matthews et al. 2010); in contrast, these insects had less preference for stony, clay, compact, and poorly-permeable soils (Parmenter & MacMahon 1984).Additionally, soil characteristics, which are strongly interrelated with climatic factors (Whittaker 1975) -mainly precipitation and temperature-can determine available humidity and thus limit organisms whose life cycles are mostly developed on the soil in arid mountain systems.
Overall, our results corroborate the idea that darkling beetles have different demands for multidimensional ecological resources in desert and semi-desert ecosystems (Yang et al. 2010, Pfeiffer & Bayannsan 2012).In other words, there is no single factor that structures their communities.In addition, due to the high percentage of recorded endemism, it is also important to consider the effect of geographic isolation and the biogeographic history of the different environments present in these mountain ranges as explanations of the current diversity reported in the darkling beetle communities.

CONCLUSION
Our study presented the first exhaustive list of darkling beetle diversity from ecoregions above 3,000 meters above sea level in Northwestern Argentina.We show that both the Puna and the Altos Andes ecoregions are home to distinctive fauna with significant endemisms.Furthermore, we show that on a regional scale and in both ecoregions geographic distance, soil heterogeneity, and climate influence the patterns of abundance and species richness of darkling beetles.The variation in elevation of the sites did not affect the abundance and richness patterns.Finally, we show that on a regional scale and in the Puna, the heterogeneity of montane grasslands and shrubs influences the abundance of darkling beetles.We suggest including darkling beetles among the focus groups in conservation studies as part of evaluating the impact of anthropogenic activities in these areas, such as mining and ecotourism, which threaten fragile high-elevation environments.

Figure 1 .
Figure 1.Map of the sampling area located in Northwestern Argentina showing the location of the sites considered for each ecoregion (circles = Puna (PU), rhombuses = Altos Andes (AA), RNUMLA= Los Andes multipurpose nature reserve área, MNAA= Abra del Acay natural monument area).

Figure 2
Figure 2. Interpolation (rarefaction) / extrapolation curves of Tenebrionidae (darkling beetles) diversity recorded in each ecoregion (a) and in each site (b), based on individuals with the same level of sample coverage, with 95% confidence intervals and 100 permutations.

Figure 3 .
Figure 3. Range / abundance curves of Tenebrionidae (darkling beetles) species collected in Puna and Altos Andes in Northwestern Argentina.