Introduction

Biodiversity studies across elevational gradients, with their long tradition in ecology and faunistics (McCoy 1990; Lomolino 2001), have evaluated characteristic patterns of diversity, community structure, and endemism and related them to climatic, spatial, historic, and/or biotic hypotheses (Körner 2007; McCain and Grytnes 2010; Grytnes and McCain 2013). Mountain ranges have also been studied to determine ecological speciation, colonization, and environmental filtering to improve our understanding of the processes that have formed mountainous communities (Graham et al. 2014; Merckx et al. 2015). Tropical mountains have been revealed as important reservoirs of biodiversity, providing a huge range of neighboring climatic and structural niches within a small geographical area (Graham et al. 2014). Moreover, elevational diversity gradients are currently of growing interest in the context of ecosystem responses to climate change (e.g., Sundqvist et al. 2013).

Although many investigations of elevational gradients of tropical arthropod diversity have been undertaken, a generally accepted pattern for all investigated taxa is lacking (cf. Körner 2007; Corcos et al. 2018). Instead, these studies have revealed various patterns ranging from monotonically decreasing diversity with increasing elevation (Wolda 1987; Rahbek 1995; Fisher 1996; Kumar et al. 2009; Lessard et al. 2011; Corcos et al. 2018), increasing diversity with increasing elevation (Sanders et al. 2003; Chatelain et al. 2018; Corcos et al. 2018), or mid-elevational peaks (Rahbek 1995; Fisher 1998; Brehm et al. 2007; Lessard et al. 2011; Plant et al. 2012; Yu et al. 2013; Bishop et al. 2014; Szewczyk and McCain 2019) suggesting that elevational patterns of species diversity should be considered on a case-by-case basis, and that various mechanisms drive these taxon-specific patterns with differing degrees of intensity. Potential explanatory factors include available elevational area (Rahbek 1995; Körner 2000; Lomolino 2001), the extension of Rapoport’s latitudinal rule to elevation (Stevens 1992), spatial isolation, temperature and moisture (both connected to productivity; Hutchinson 1959; Wright 1983, Yu et al. 2013; Szewczyk and McCain 2019), seasonality (Plant et al. 2018), vegetation, overall zoogeographical distribution range (Stevens 1992), predation and parasitism pressure, and habitat heterogeneity (Szewczyk and McCain 2019). Whereas several studies have addressed single factors or combinations thereof (Yu et al. 2013; Corcos et al. 2018; Szewczyk and McCain 2019), few have used different niche metrics to evaluate the influence of both evolutionary and ecological determinants on local diversity patterns. Perhaps such an approach has rarely been attempted because it requires knowledge of the full taxonomy at the species level, which is often unknown for arthropods in tropical regions. The advantage of calculated niche dimensions, however, is that they integrate a number of single explanatory evolutionary and ecological factors that usually co-vary along elevational gradients and are thus difficult to disentangle (Szewczyk and McCain 2019).

We present the results of a quantitative field study on staphylinid beetles of the subfamily Steninae across an elevational gradient of a tropical mountain area in northern Thailand. Megadiverse Steninae are predators in moist habitats and forage in various zones; they are (1) inhabitants of moist humus and plant debris near the ground, (2) plant climbers, or (3) surface runners on bare ground. Moreover, many “Dianous” species seem to be closely associated with (4) waterfalls and mountain rivers. In addition, certain ecomorphs can be distinguished because they correlate with these habitats in terms of their morphology (reviewed in Betz et al. 2018).

Our study is one of the few on tropical insects that focuses on a clearly defined lineage of closely related species (subfamily consisting of two genera) in combination with full taxonomic species resolution. We analyze species richness, heterogeneity, niche occupancy and community composition/similarity of the Steninae fauna occurring along an elevational gradient from approximately 400–2500 m. Aiming to maintain equal sampling efforts at each spot, we collected Steninae from a variety of habitat types, identified them to species level, analyzed the alpha and beta diversity of the sampled communities and calculated niche parameters of the occurring species. We addressed the following aims and hypotheses.

Aim 1—Species richness and species heterogeneity in relation to elevation: Previous work on elevational gradients has often revealed overall declining numbers of species and species diversity with increasing elevation, whereby possible hump-shaped declines might occur that are partly attributable to species overlap among zonal communities (cf. Lomolino 2001; Colwell et al. 2004; Dunn et al. 2007). Our aim was to determine whether such decline in species richness and heterogeneity also exists in the investigated Steninae, and whether it better fits a linear or a quadratic function showing a mid-elevational peak. In addition to documenting species richness and diversity patterns, we were interested in the changes of actual species composition along the elevational gradient, since the changing environmental conditions (e.g., temperature) along elevational gradients may have created defined beetle communities. Hypotheses: According to the prevailing influence of the mid-domain effect inherent in the zonal structuring of a mountain (Colwell et al. 2004; Dunn et al. 2007), species numbers are highest at mid-elevations and so can be best described via quadratic instead of linear functions. However, we expected that factors such as the elevation-dependent availability of certain habitats also influence the final pattern; therefore, we tested this hypothesis for single pre-defined habitats such as bare ground, litter, vegetation, and sprinkled rocks associated with streams. For beta diversity, at higher elevations, dissimilarity should be mainly determined by species loss from lower to higher collection sites, whereas at lower elevations species turnover (i.e., replacements of “lowland species” by “highland species”) should prevail.

Aim 2—Elevation and habitat niche partitioning: To evaluate potential causes that sustain numbers of species across the elevational gradient, niche partitioning among members of the assemblages was correlated with elevation and habitat. Hypotheses: For niche structuring, we expected to find both elevational and habitat specialists indicative of effective environmental filtering determining the composition of local communities. The existence of elevational specialists might indicate a prevailing influence of temperature in these species as recently established for ants considering species-level responses (Szewczyk and McCain 2019).

Aim 3—Climatic-niche conservatism: A prevailing influence of elevation as a proxy for temperature might be attributable to climatic-niche conservatism of certain clades, i.e., their phylogenetic tendency to retain ancestral niche attributes (Boucher et al. 2014). We tested whether differences occurred between Paleartic and Oriental species in terms of their elevational ranges. Hypothesis: Having colonized these mountain ranges from the more northern regions of the Himalayas, Palearctic species show a wider distribution across higher elevations compared with Oriental species that predominate at lower elevations.

Aim 4—Niche breadth as a predictor for range size: Niche breadth is predicted to be positively correlated to the geographical range of a species, with certain species being common and others rare (Slatyer et al. 2013). For mountain ranges, we can assume that elevational generalists have a broader physiological tolerance toward various abiotic conditions such as temperature and can thus be considered candidates that more easily colonize entire mountain ranges. Hypothesis: There is a correlation between niche breadth and species frequency (i.e., the number of locations of a species divided by the total number of collection sites), i.e. generalists are expected to show higher species frequencies than specialists.

We finally collate and discuss our results to widen our understanding of the influence of factors such as orography, colonization, environmental filtering, and ecological radiation on the formation of such species assemblages.

Materials and methods

Studied organisms

The Steninae is a subfamily within the Staphylinidae (rove beetles) and contains only two genera, i.e., Stenus Latreille 1797 and “Dianous”Footnote 1 Leach 1819. The biology of Steninae has recently been reviewed in Betz et al. (2018). With about 3000 species, Stenus is one of the largest beetle genera on earth and is distributed worldwide, except for New Zealand (Puthz 1971). The genus “Dianous” comprises > 300 species (Shi and Zhou 2011) and is present in Oriental, Palaearctic, and Nearctic regions (but not in the Ethiopian region) with its main distribution in Asia (India, China, and southeast (SE) Asia). Its center of distribution lies in the mountain areas between the Palaearctic and Oriental regions south of 31° northern latitude (Indochina Peninsula and southern slopes of the Himalaya) (Puthz 1981a). From here, it has probably dispersed in (1) a south-eastern direction, (2) a north-western direction, and (3) via the Bering Straits to North America (Puthz 1981a). We chose this subfamily for our study because of the extensive taxonomic and ecological knowledge of this clade and its high diversity, especially in the SE Asian tropics.

Study area

The study was conducted during the dry season from Dec 2013 to Feb 2014 in the highest and second highest mountain areas of Thailand, i.e., the Doi Inthanon (ITH) (18° 35′ 0″ N, 98° 29′ 0″ E) and the Doi Pha Hom Pok (PHK) (19° 59′ 16″ N, 99° 08′ 47″ E) National Parks in Chiang Mai Province, North Thailand. The mountains of northern Thailand belong to the Indo-Burma biodiversity hotspot identified by Myers et al. (2000) as being of global significance and showing a mixture of Palearctic and Indo-Malayan influences (cf. Plant 2009, 2014; Plant et al. 2011). The mountain area of Doi Inthanon (2565 m) (main site of our field studies) arose as a metamorphic core complex bounded by low-angle faults; it was uplifted 50 MYA and is now higher than the surrounding mountains (cf. MacDonald et al. 1993). Consequently, it shows a larger elevational extent supporting a greater range of elevational habitats than elsewhere in northern Thailand (Plant 2009). The Doi Pha Hom Pok mountain area (2285 m) (our second collection area) is the second highest mountain in Thailand; it belongs to the Dan Lao Mountain range. According to Chayamarit and Puff (2007), the climate in both mountain areas is monsoonal, i.e., it is divided into a rainy season between May and September, and a dry season between October and April. Humidity varies between 60% (dry season) and 80% (rainy season), whereby low clouds and fog occur almost daily, especially at higher elevations.

Several collection sites (mostly associated with mountain streams and waterfalls: see Online Resource 1) were selected based on elevation, i.e., (1) 400–900 m (lowland), (2) 900–1400 m (lowland–lower montane), (3) 1400–1900 m (lower montane), (4) 1900–2500 m (lower montane–upper montane) (Online Resource 2). These elevational zones approximately reflect major vegetation transitions and cover various primary vegetation types comprising the (Dry) Deciduous Dipterocarp Forest, the Evergreen Lower Montane Rain forest (including Pine stocks), the Evergreen Upper Montane Cloud Forest, and the Upper Montane Peat Bog (cf. Chayamarit and Puff 2007; Betz et al. 2014). Whereas the 46 collection sites mostly represented primary natural sites, four sites comprised secondary/anthropogenic habitats such as rice fields, artificial ponds, and secondary grassland (Online Resource 2).

Sampling and identification

Each collection site listed in Online Resource 2 was hand-sampled. The numbers of sampled individuals and species were recorded at each sampling site and collecting date. To ensure equal sampling among sites, each sampling event lasted for 2.5 h and involved four experienced insect collectors who focused equally on three active collecting techniques and habitat types: (1) Winkler leaf litter sieving, (2) sweeping a linen net in vegetation, and (3) actively collecting Steninae directly from sprinkled rocks, logs, and stream banks with a mouth-pooter. The species collected by the different methods were placed in separate vials to allow the various species and associated number of collected specimens to be assigned to the four pre-defined habitat types; this was necessary for calculating the niche parameters. All sampled Steninae specimens were killed with ethyl acetate and preserved in 70% ethanol for later identification. The pre-identification of the specimens was performed by using the keys of Puthz (1980, 1981a, b, 2000, 2013) and de Rougemont (1981, 1983, 1987). The final verification of these identifications and description of new species was undertaken by Puthz (2015). Over the three months at the Doi Inthanon and Doi Pha Hom Pok National Parks, we performed a total of 33 and 13 sampling events, respectively.

Data analysis

Aim 1—Species richness and species heterogeneity in relation to elevation: In addition to (1) species richness, i.e., the number of species per sample, (2) maximum number of species, i.e., the maximum number of species per collection site at each elevational zone, (3) species density, i.e., total number of species per elevational zone divided by the number of collection sites at this zone, and (4) species frequency, i.e., the number of collection sites where a species was collected divided by the total number of collection sites, we calculated (5) Simpson’s index, i.e., a measure of species diversity (heterogeneity), with the software Ecological Methodology v. 7.2 (Krebs 2011; see also Krebs 1999). In order to avoid Simpson’s index being affected by too low sample sizes, we calculated it for a sampling unit only when it met the following two criteria: number of species > 2 and total number of individuals per sampling unit > 10. Finally, to compare species similarity across pre-defined elevational zones, (6) various aspects of the Sørensen index were calculated with R software (see below). In particular, the following analyses were performed:

  1. 1.

    Evaluation of habitat- and elevation-specific species assemblages: A discriminant function analysis (DA) was employed separately for “Dianous” and Stenus by using SPSS 25 to analyze whether our pre-defined four elevational classes and habitat types could be distinguished by the potentially discriminating species composition. To reduce the influence of rare species, only species that occurred in more than two samples were considered, leaving 44 species (16 “Dianous”, 28 Stenus) in the analysis (cf. Online Resource 3). To reduce the number of discriminating variables (i.e., species) in relation to the number of cases (i.e., collection sites), a principal components analysis (PCA) (varimax rotation) of the log10-transformed number of collected specimens per species was conducted prior to the DA. Only the loadings of single species > 0.5 and < −0.5 were considered for the interpretation of the PCs. The extracted PCs with eigenvalues > 1 (elevation: six for “Dianous”, 11 for Stenus; habitat: five for “Dianous”, 11 for Stenus) were then fed into the DA.

  2. 2.

    Decrease of species richness versus decrease of species density with elevation: In order to test the hypothesis that declining species richness at higher elevations simply reflects the reduced land area available, we compared the regression slopes between the species richness and species density with elevation. A shallower slope in species density with elevation compared with the absolute number of species is indicative of such an area effect. We therefore sub-divided the entire elevational range into eight equal zones (i.e., 400–650, 650–900, 900–1150, 1150–1400, 1400–1650, 1650–1900, 1900–2150, > 2150). Response variables were log10(x+1) transformed to improve homoscedasticity (because values between 0 and 1 occurred in the case of the species densities, we added the value 1 to each measurement of both variables beforehand), and conducted a Generalized Linear Model (GLM) regression by using SPSS 25 with three variables, i.e., (1) a dummy variable “data set” coding for the variables “number of species” or “species density”, (2) the response variable containing concrete values of number of species and species density (assigned to the codings of the first variable) and (3) “elevation” as a covariate. The model included “data set” and “elevation” as main effects and “data set” * “elevation” as an interaction term, i.e., y ~ elevation + ”data set” + elevation * “data set”, whereby y represents number of species and species density. In such a model, a significant interaction term indicates that the slopes of both the data sets differ from each other. Since neither species richness nor species density were on the same scale (species density has by definition much lower values compared with number of species), we scaled them prior to analysis by dividing each value of the dependent variable by the mean of this variable.

  3. 3.

    Diversity-elevation and richness-elevation relationships: To test whether a linear relationship exists between species diversity and elevation, a linear regression analysis was conducted on the double logarithmic data set (homoscedasticity confirmed by z standardized plots of the predictors against the residuals). For species richness (comprising count data including zero values), GLMs were calculated in RStudio (R version 3.6.1), assuming a negative binomial distribution. The adequacy of this distribution was confirmed by evaluating the dispersion parameter (usually close to 1) and by model simulations (whose predicted zero portions corresponded well to those observed). In the case of significant results of the linear and quadratic model components, a Likelihood Ratio Chi2 test was used to test whether the quadratic model fits the data better than the linear one. These analyses were performed for both the genera (“Dianous” and Stenus) separately and combined.

  4. 4.

    Similarity measures between elevations and mountain areas: Dissimilarity was measured with the Sørensen index. This means that we compared the similarity between (1) both the national parks in general (without consideration of elevational ranges) and (2) different elevations. These analyses were performed for both the genera (“Dianous” and Stenus) separately and combined. We used the RStudio package betapart (Baselga et al. 2018) employing the function ‘beta-pairwise.R’ to compute total dissimilarity with the Sørensen index and its respective turnover (replacement) and nestedness components (Baselga 2010, 2012; Baselga and Orme 2012).

Aim 2—Elevation and habitat niche partitioning: To compare the various species with respect to their elevational and habitat preference, the Shannon–Wiener Evenness was calculated with the software Ecological Methodology v. 7.2 (Krebs 2011; see also Krebs 1999) for each species as an indicator for its niche breadth. In terms of their elevational occurrence, niche breadth was referred to the four elevation classes defined above. With respect to their habitat preference, niche breadth was applied to the following habitat types typical for Steninae ecomorphs as determined in earlier studies (reviewed in Betz et al. 2018): (1) bare or sparsely grown over (sandy or muddy) ground on stream banks, (2) litter layer and debris, (3) shrub-like and grassy vegetation, and (4) sprinkled (mossy) rocks and dead logs in streams and waterfalls. Whereas habitat types (2) and (3) were abundant across the entire elevation gradient, types (1) and (4) continuously decreased in abundance with higher elevations. To minimize the effect of overall low beetle abundances on the analysis, niche breadth was calculated only for species for which > 10 specimens were collected (cf. Aguiar et al. 2013).

The habitat preference of “Dianous” species for sprinkled rocks in streams and waterfalls (type 4 above) was tested by a Chi-squared test (conducted in SPSS 25) assuming equal distribution across the four habitats as our null hypothesis. This test was based on the patterns of presence-absence of “Dianous” species in each habitat type (see Fig. 4a).


Aim 3—Climatic-niche conservatism: Climatic-niche conservatism was tested via a fourfold Chi squared test (conducted in SPSS 25), considering the zoogeographic realm of each species in relation to its elevational occurrence. The mountain range was subdivided in such a way that both the lower pre-defined elevational zones formed one elevational category (350–1400 m) and both the upper zones the second category (1400–2500 m). We tested whether more than expected Palearctic species were found in the upper category, or more Oriental species occurred in the lower one. The zoogeographic affiliations were obtained from Herman (2001) and Löbl and Löbl (2015). This analysis was restricted to the members of the genus Stenus, since all the collected “Dianous” species belong to the Paleartic realm.


Aim 4—Niche breadth as a predictor for range size: The potential correlation of niche breath with species frequency was tested by bivariate nonparametric Spearman Rank correlation analysis. For the habitat niche, this analysis was performed for the entire mountain range and for each of four pre-defined elevational range separately.

Results

Aim 1: Species richness and species heterogeneity in relation to elevation

Overall numbers of species

In ITH and PHK, we found 22 “Dianous” and 60 Stenus species (Online Resource 3). At ITH, we collected 17 “Dianous” (1785 specimens) and 46 Stenus (1004 specimens) species (three collection trips), and at PHK, 14 “Dianous” (141 specimens) and 41 Stenus (218 specimens) species (one collection trip). Among the “Dianous” species at ITH, the highest species frequencies were attained by D. srivichaii (45%), D. obliquenotatus (30%), and D. karen (27%) (Online Resource 3). At PHK, the highest species frequencies were 23% for D. biscobriculifrons, D. karen, D. luteolunatus, D. nitidicollis, and D. srivichaii (Online Resource 3). Regarding Stenus, the highest frequencies at ITH were detected for S. rorellus cursorius (27%), S. srisukai (21%), and S. liangtangi and S. thanonensis (18%). At PHK, the highest frequencies were observed for S. explanipennis, S. pustulatus, and S. schwendingeri (31%) (Online Resource 3).

Among all species collected, six “Dianous” and one Stenus species were new species discoveries (Puthz 2015).

Evaluation of habitat- and elevation-specific species assemblages

Elevational ranges: In “Dianous”, none of the six extracted PCs significantly contributed to group separation in the DA, i.e., the “Dianous” species did not form characteristic assemblages distinguishable at the various elevational ranges. In Stenus, two (PCs 3 and 9) of the 11 extracted PCs showed significant (P < 0.001) χ2 values for group separation; their correlations (loadings) with the respective species are shown in Online Resource 4. The DA extracted three discriminant functions that significantly contributed to the separation of elevational ranges, explaining 57% (function 1), 26.6% (function 2), and 16.4% (function 3) of total variation (Fig. 1). Whereas discriminant function 1 mainly separated the highest elevation range (1900–2500 m) from the lower ones, function 2 mainly segregated the second highest range (1400–1900 m) (Fig. 1).

Fig. 1
figure 1

Discriminant analytical separation of the four predefined elevational ranges in the collected Stenus species. Each dot represents a different collection site. The numbers of collected specimens per species were log10-transformed prior to the analysis and fed into a PCA to reduce the number of variables. In the following DA, two of the 11 principal components (PC) extracted in this way significantly contributed to group separation and could thus be used as predictors for the elevational ranges. Only the first two functions are plotted, showing the following standardized canonical discriminant coefficients: Discriminant function 1 (eigenvalue = 2.06; canonic correlation: 0.82): PC3: -0.38; PC9: 0.91; Discriminant function 2 (eigenvalue = 0.96; canonic correlation: 0.70): PC3: 0.87; PC9: 0.45. The species assemblages associated with the PCs (cf. Online Resource 4) are provided in the diagram according to their association with the elevational ranges. Meaning of symbols: open triangle: 400–900 m, open square: 900–1400 m, open circle: 1400–1900 m, filled circle: 1900–2500 m. Large bold squares represent the group centroids

Full size image

Habitats: In “Dianous”, two of the five extracted PCs significantly (P < 0.05) contributed to group separation in the DA (Online Resource 5). Only the first discriminant function significantly contributed, although weakly, to the separation of the habitat containing waterfall- and stream-associated rocks, explaining 99.7% (eigenvalue: 0.36; canonical correlation: 0.51) of total variation. Both the extracted PCs (analysis given in Online Resource 5) showed “Dianous” species assemblages characteristic for this habitat type, the assemblage associated with PC3 more strongly contributing to this separation (cf. details in Online Resource 5). In Stenus, four of the extracted 11 PCs significantly (P < 0.05) contributed to group separation in the DA (Online Resource 6). Only the first discriminant function significantly contributed, but weakly, to the separation of the litter habitat, explaining 79.5% (eigenvalue: 0.34; canonical correlation: 0.50) of total variation. The four PCs extracted in Online Resource 6 show Stenus assemblages characteristic for the litter layer, with all assemblages equally contributing to this separation (cf. standardized canonical discriminant coefficients in Online Resource 6).

Decrease of species richness versus decrease of species density with elevation

The number of species per collection site varied from 1 to 20 and 1 to 22 in ITH and PHK, respectively (cf. Online Resource 2 for localities). Regarding the pre-defined elevational ranges, an inspection of (1) total numbers of species, (2) maximum number of species per sample, and (3) species densities confirmed decreasing species richness with increasing elevation, especially from 1400 m upward (Fig. 2). GLM analysis revealed a significant main effect of elevation on number of species and a significant interaction between elevation and both number of species and species density (denoted by the dummy variable “data set”) (Table 1). This indicates that species density declines to a lesser extent than number of species with elevation.

Fig. 2
figure 2

Total numbers of species (white column), maximum number of species per sample (gray column), and species densities (i.e., total numbers of species divided by the number of collection sites at each elevational zone) (black column). The shallower decline of species density compared with the number of species was statistically confirmed by a GLM regression (cf. Table 1). All numbers refer to both “Dianous” and Stenus and both the national parks together. The small numbers above the columns represent the exact values

Full size image
Table 1 Results of Generalized Linear Model (GLM) regression testing for differences in the decline of number of species and species density (i.e., the total number of species divided by the number of collection sites at the respective elevational zone) with elevation
Full size table

Diversity-elevation and richness-elevation relationships

Whereas species diversity (Simpson’s Index) showed no significant relationship with elevation, the numbers of species per collection site (not differentiated by habitat) declined with elevation in a linear fashion (Table 2, Fig. 3a). The differentiation of numbers of species according to microhabitats revealed a linear decline with elevation for the vegetation-inhabiting species (P < 0.05) (Fig. 3b), whereas the litter-inhabiting species (Fig. 3c) (P < 0.01) and the species associated with sprinkled rocks at streams and waterfalls (Fig. 3d) (P < 0.001) followed a quadratic function showing a mid-elevational peak (Table 2). Our testing of these relationships separately for “Dianous” and Stenus revealed that these relationships occurred for “Dianous” (all microhabitats together, waterfalls) and Stenus (litter, vegetation, waterfalls) (Table 2).

Table 2 GLM results for testing linear and polynomial quadratic models on the relationship between species richness and elevation (cf. Fig. 3)
Full size table
Fig. 3
figure 3

Linear (a, b) and quadratic (c, d) GLM regressions of number of species as related to elevation, both “Dianous” and Stenus species combined, and both the National Parks considered. a All microhabitats considered together, b–d Separate analyses for the various habitats. For statistical GLM results see Table 2. Each dot refers to a certain collection site

Full size image

Similarity measures between elevations and mountain areas

Although the sampling intensity at ITH was 2.6 times higher than that at PHK, the total number of “Dianous” and Stenus species was only slightly lower at PHK (see above). The overall Sørensen dissimilarity (βsor) of the Steninae/”Dianous”/Stenus species composition of the two National Parks was 0. 366/0.333/0.379. The respective dissimilarity components were 0.321/0.267/0.341 (spatial turnover βsim) and 0.0452/0.0667/0.0378 (nestedness βnes).

Across the elevational zones, the similarity in species composition decreased with increasing elevation (Table 3). In the transition from the lowest (400–900 m) to the next highest elevational zone (900–1400 m), dissimilarity was mainly explained by species turnover. However, across higher elevations, the nestedness component (i.e., the omission of certain or, in the case of “Dianous”, almost all species) became steadily more important as an additional (Stenus) or exclusive (“Dianous”) explanatory factor (Table 3).

Table 3 Compositional Sørensen dissimilarity of the Steninae communities between the four elevation ranges of both ITH and PHK
Full size table

Aim 2: Elevation and habitat niche partitioning

Overall habitat and elevational preferences

In terms of the preferred habitats of “Dianous” beetles, they showed strong hygrobiont preferences (petri- and bryomadicolous including dead logs and rocks immersed in water) associated with the spray zone of running water and cascades. Whereas the representatives of most species were usually observed to run actively along the waterline of sprinkled rocks (e.g., D. betzi, D. obliquenotatus) or to sit in water-sprayed mosses (e.g., D. meo, D. tubericollis), others (D. chetri) were encountered in clusters of several individuals at the underside of logs sprinkled by water. Less frequently, “Dianous” individuals were collected in wet litter close to streams and cascades. Although Stenus species were rarely collected from the spray-zones of rocks (e.g., S. aestivalis), most individuals showed strong associations with moist litter around streams. Here, they were found either among debris and single leaves (e.g., S. notaculipennis) or on sieving the thick humus layer at higher elevations (e.g., S. luteomaculatus). Individuals of other species were swept in higher numbers from brush-like (Elatostema) vegetation (e.g., S. angusticollis). Whereas riparian surface runners on bare ground were almost absent from our collections, we frequently encountered higher numbers of specimens on relatively dry sandy soil between grassy vegetation (e.g., S. diversiventris, S. rorellus cursorius).

Although several Stenus species could be collected at higher elevations around 2000 m (two species exclusively above 1400 m), most “Dianous” species were only found up to 1400 m (except for three species that were sieved from litter at higher elevations) (Fig. 4a).

Fig. 4
figure 4

Niche occupation of collected Steninae. Only species with > 10 collected individuals are considered. a Distribution of sampled “Dianous” and Stenus species across elevations [m a. s. l] (left column) and microhabitats (right column) in both Doi Inthanon and Doi Pha Hom Pok National Parks. The habitat specifications refer to (1) litter layer and debris, (2) sprinkled (mossy) rocks and logs in streams and waterfalls, (3) shrub-like and grassy vegetation, and (4) bare or sparsely overgrown (sandy or muddy) ground of stream bank. Numbers in brackets behind species names refer to the labels used in b. b Graphic depiction of niche breadth parameters (habitat versus elevation) of the species listed in a. The Shannon–Wiener niche breadths were taken from Online Resource 3

Full size image

Shannon–Wiener niche widths

With respect to elevation, “Dianous” preferred the lower elevations (up to 1400 m a.s.l.), as reflected in relatively low niche breadths (mostly considerably smaller than 0.5; Online Resource 3). In a few species, only single individuals were collected in higher regions. Among the collected Stenus species, we observed a greater differentiation between species, i.e., some species were restricted to lower (e.g., S. humeralis) or to higher (e.g., S. luteomaculatus) elevations (resulting in low niche breadths), whereas others (e.g., S. virgula) occurred across the entire mountain range from the lowest elevations to the summit. Correspondingly, these species showed considerably higher niche breadths (Online Resource 3).

With respect to habitat, “Dianous” clearly preferred the sprinkled (mossy) rocks and logs in streams and waterfalls (mostly resulting in low niche breadths) (Fig. 4a) (Chi-squared test, χ2 = 13.47, df = 3, P < 0.005; N = 30, standardized residual for rocks: 3.1). Hence, in graphic representations of niche breadth parameters, most “Dianous” species occurred at the extreme left (Fig. 4b), corresponding to low Shannon–Wiener evenness values. In Stenus beetles, we observed various habitat preferences and niche widths. Some species were collected by litter sieving only (e.g., S. luteomaculatus). Consequently, they showed the minimum possible habitat niche width, namely zero. Various species (with higher niche widths) were also found in other microhabitats, some of them additionally being swept in higher numbers from vegetation (e.g., S. thanonensis). The highest niche widths were attained by species collected from all four defined microhabitats (e.g., S. amoenus). Stenus species restricted to sprinkled rocks or sparsely vegetated ground (cf. the European S. comma) could not be detected among the collected specimens. Only S. beesoni were found under dry leaves lying upon rocks close to the waterside.

Aim 3: Climatic-niche conservatism

In Stenus species, our 2 × 2 table test revealed a significant connection between the relatedness of zoogeographic realm (Oriental versus Palearctic) and elevation (Pearson’s χ2 = 3.94, df = 1, P < 0.05; N = 76). According to the standardized residuals, this can be assigned to more Palearctic (standardized residual = 1.0) and fewer Oriental species (standardized residual = −1.4) than expected in the upper elevational zones between 1400 and 2500 m.

Aim 4: Niche breadth as a predictor for range size

The Rank correlation analyses between Shannon niche width and species frequency revealed low to moderately positive correlations for elevation (Spearman correlation coefficient: 0. 56 (P < 0.001), n = 35; Fig. 5) and habitat (Spearman correlation coefficient: 0.37 (P < 0.05), n = 35). These results also held true when “Dianous” (n = 16) and Stenus (n = 19) were analyzed separately (Spearman correlation coefficients for elevation: “Dianous”: 0.53, P < 0.05/Stenus: 0.64, P < 0.001), although the Shannon–Wiener Evenness for “Dianous” revealed no significant results for habitat (Spearman correlation coefficients for habitat: “Dianous”: P = 0.14/Stenus: 0.46, P < 0.05). When these correlation analyses for habitat were performed within each elevational range, significant correlations were seen between species frequency and Shannon niche width for Stenus (not for “Dianous”) at 400 -900 m (Spearman correlation coefficient: 0.64, P < 0.05, n = 14) and at 900 - 1400 m (Spearman correlation coefficient: 0.54, P < 0.05, n = 14).

Fig. 5
figure 5

Spearman Rank correlation between species frequency (i.e., the number of locations per species divided by the total number of collection sites) and elevational niche breadth. The Shannon–Wiener Evenness was taken from Online Resource 3. The species labels refer to the labels given in Fig. 4a

Full size image

Discussion

Mountain slopes resemble natural laboratories with gradients of ecological factors in a condensed format, making it possible to study the reaction norms of single clades and entire communities along continuously changing environmental parameters. The division of the main drivers nevertheless represents a great challenge, depending on their inter-connectedness. Indeed, both temperature and elevational area have previously been identified as major factors influencing numbers of species and community structure (Bishop et al. 2014; Wardhaugh et al. 2018). However, historical factors, i.e., climate change, niche conservatism, orogenesis, and immigration and evolutionary mechanisms (speciation and radiation), need to be considered (Plant 2009; Plant et al. 2011, 2012; Graham et al. 2014; Merckx et al. 2015) in addition to ecological explanations. All these factors are usually interwoven, and hence causal analyses require solid knowledge of the biology and evolutionary history of the organisms under study.

Although only two mountain regions were considered during our three-month field study, we collected 82 species (including seven new species descriptions), representing 56% of the known Thai Steninae fauna and 61% of that for Northern Thailand, thus demonstrating the high value of this region as a biodiversity hotspot. The detected large numbers of species might also be positively correlated with the high degree of endemism found in this area (Plant 2014). Although the current state of knowledge of the zoogeographic distribution of Steninae in Northern Thailand (and bordering countries) does not allow us to identify endemites with certainty, six of the collected “Dianous” and seven Stenus species are only known from the sampled mountain region.

The sampling events at PHK were 2.6 times lower than at ITH and comprised six times fewer individuals. Nevertheless, the total numbers of species collected in the two regions were similar. Since PHK lies 150 km further to the northwest of PHK, this is indicative of an increasing Steninae diversity toward the foothills of the Himalayas. The overall Sørensen similarity between the two regions was 63% (18 species were collected at PHK only) suggesting a high local diversification in this region. Accordingly, only six “Dianous” and seven Stenus species occurred in species frequencies of 20% or higher, whereas the vast majority occurred in only a few sampling events. Sørensen dissimilarity between these two mountain regions clearly revealed higher turnover values (compared with nestedness) further supporting the view that the overall dissimilarity between these two mountain regions is mainly caused by species turnover.

Our study has enabled us to assign each collected individual to a valid species, so that niche widths could be calculated, compared and integrated in the broader context of elevational gradients in species richness and community structure across a tropical mountain range. Whereas the integration of ecology and evolution has yet to be achieved in the study of elevational diversity gradients, evaluating niche metrics has the potential to contribute to our general understanding of the influence of both evolutionary and ecological determinants on local diversity patterns. We addressed four major aims connected to certain hypotheses that involved aspects such as niche differentiation and niche conservatism and their potential relationship to species frequency, environmental filtering, and colonization of the considered mountain range.

Aim 1: Species richness and species heterogeneity in relation to elevation

As we expected, and as established in other elevational gradient studies (Gathorne-Hardy et al. 2001; Hillebrand 2004; Malsch et al. 2008; Werenkraut and Ruggiero 2011), numbers of species clearly declined at higher elevations, particularly from 1400 m upwards (cf. Figs. 2, 3). One general physical reason for declining numbers of species at higher elevations lies in the reduced land area available for biota because of the steeper terrain toward the summits (Rahbek 1995; Körner 2000; Lomolino 2001). In comparison with absolute numbers of species, the much lower trend for decreasing species density (black bars in Fig. 2; Table 1) with elevation suggests that the observed elevational gradient in the number of species is superimposed by such an area effect. Notwithstanding, additional factors are probably involved because, for both of the highest elevational zones, our data show a slight decrease, even of species density (black bars in Fig. 2). Temperature is a major factor that decreases with elevation, influencing insect performance in various direct (e.g., metabolic activity, growth) and indirect (e.g., productivity, species interactions) ways (Hodkinson 2005); decreasing temperatures are usually assumed to act as an environmental filter that only allows physiological specialists to settle at higher elevations [cf. climatic hypotheses for elevational gradients in Grytnes and McCain (2013)]. Indeed, in additional laboratory observations, Stenus luteomaculatus sieved from the thick humus layer of the upper mountain regions only (cf. Fig. 4) had a clumsy appearance combined with a low consumption rate indicative of an “energy minimizing” strategy that possibly evolved in adaptation to a diminished food supply, as previously discussed for some Central European Stenus species (Betz and Fuhrmann 2001).

Since numbers of species are high, even at the lowest elevational range (400–900 m), no clear mid-elevational peak of species richness, as established for other groups of insects (Rahbek 1995, 2005; Fisher 1996, 2002; Lomolino 2001; Colwell et al. 2004; Plant 2009; Werenkraut and Ruggiero 2011; Plant et al. 2012; Yu et al. 2013), is visible when all microhabitats are taken together (Fig. 3a). However, if the four microhabitats are considered separately, the litter-inhabiting species (mostly Stenus, for which this habitat type is predominant among the collected species) and the species associated with rocks and logs in streams and waterfalls (mostly “Dianous”, for which this habitat type is predominant among the collected species) show a significant peak at mid-elevations (cf. Fig. 3c, d). Regarding the rock- and dead-log-associated species at waterfalls (Fig. 3d), this is clearly the result of the absence of richly structured waterfall environments at elevations above 1000 m. Consequently, in our DA, we have not established any elevation-specific species assemblages for them. In the litter-dwelling species, the elevational peak is higher (~ 1500 m) (Fig. 3c), corresponding to the evergreen lower montane rain forest being associated with a stable condensation (cloud) zone (Chayamarit and Puff 2007). Even during the dry season, these conditions provide more permanently moist soils (i.e., stable conditions) with favourable temperatures at these elevations. In addition, as suggested by the similar community structures between adjacent elevational ranges (cf. Table 3), such mid-elevational peaks might be simply accentuated by feedback among zonal communities (mid-domain effect), i.e., an elevational zone sandwiched between two species-rich zones might gain additional species richness by the turnover of species from both these zones (Lomolino 2001; Colwell et al. 2004; Dunn et al. 2007). In the vegetation-dwelling species (predominated by Stenus), we found a linear decline with elevation (Fig. 3b) (no mid-elevation peak). Here, several of the collection sites at the lowest elevations showed high numbers of species, so that this trend might be driven by the falling temperatures with increasing elevations and a possible shift of ground vegetation toward shrub-like (e.g. Elatostema) forms that might be colonizable by some specialists only (such as Stenus thanonensis beetles). Moreover, this trend might mirror the gradient of the declining area of the vegetation zone with elevation.

Unlike the numbers of species and densities, Simpson’s diversity index showed no significant trend with elevation. This result is probably partly attributable to some communities at lower and middle elevations being dominated by few species occurring in high numbers, thus biasing these indices (e.g., locations 1, 5, 19 in Online Resource 2).

In accordance with other such studies on insects (Kumar et al. 2009), community similarity decreased with increasing distance between the pre-defined four elevational ranges (cf. Table 3). On the one hand, this result is correlated with the decreasing numbers of species at higher elevations (notably, “Dianous” species tend to be totally absent from the highest elevational zone: Fig. 4). On the other hand, it also reflects the restricted elevational distribution of several species (cf. Fig. 4) suggesting (particularly for Stenus, which occurs up to the summit) a high species turnover from the species-rich lower elevations toward a species-poor but distinct montane fauna (see also Kumar et al. 2009 for paper wasps). This view is supported by (1) both components of the Sørensen index, which indicate that, in Stenus, the beta diversity toward the higher mountain range (> 1400 m) is determined in almost equal parts by species turnover and nestedness (Table 3) and (2) the results of our DA revealing Stenus species (assemblages) restricted to the two highest elevational ranges (cf. Fig. 1). These species might be especially adapted to the lower temperature regimes of these elevations and can thus be used as predictors for them. Interestingly, across the lowest (400–900 m) and the second lowest (900–1400 m) elevational zones, beta diversity is mainly determined by species turnover, whereas nestedness is negligible (Table 3). Taken together, these patterns indicate that, across the elevations, a replacement of some typical “lowland species” by “mountain species” occurs (cf. Fig. 4), whereby this “highland fauna” thins out with increasing elevation.

Aim 2: Elevation and habitat niche partitioning

Our quantitative niche breadth calculations regarding both elevation and habitat (cf. Fig. 4b, Online Resource 3) revealed different patterns of niche occupation amongst the investigated Steninae, i.e., some species showed (1) narrow niches for both habitat and elevation (e.g., S. luteomaculatus; present in the thick humus layer at upper elevations only), (2) a broad habitat niche combined with a medium elevational niche (e.g., S. srisukai; present in all of the pre-defined habitats, except for rocks and logs in running water), (3) a narrow habitat niche combined with a broad elevational niche (e.g., S. virgula; mainly associated with the litter layer but distributed across the entire elevational range), and (4) a broad habitat niche combined with a narrow elevational niche (e.g., S. pustulatus; present in all pre-defined habitats but not above 1400 m). Other species showed intermediate patterns (e.g., S. malickyanus). This variety of niche breadths, which is indicative of the co-occurrence of specialists and generalists, largely contributes to the ecological diversity across these mountains. If, for the collected Stenus species, the elevational distribution is a proxy for their thermal niche, the restricted (low or high) elevational preferences observed in several species are an indication that thermal niche differentiation has played a role in their evolution, lowering their dispersal (and thus gene flow) across the elevations. This agrees with hypotheses and studies predicting that rapid species diversification is associated with climatic-niche evolution, an aspect that is especially true in regions with limited seasonality such as the tropics (Janzen 1967; Ghalambor et al. 2006; Kozak and Wiens 2010).

In contrast to “Dianous” beetles, which are largely specialized on waterfall habitats, the litter layer and debris on the ground form the predominant habitat for the collected Stenus species (cf. Figure 4a, Online Resource 3). This is supported by our DA, which, according to species composition, separated the litter habitat from the other three habitats and established four well-distinguished Stenus assemblages, each consisting of species often co-occurring (cf. Online Resource 6). The (thick) humus layer probably provides the only habitat that can be colonized by Stenus beetles at higher elevations (1400–2400 m), with their lower temperatures.

Aim 3: Climatic-niche conservatism

All the collected “Dianous” and > 50% of the collected Stenus species belong to the Palearctic realm, and they or their ancestors have thus probably colonized the northern Thailand mountain area from (cooler) Himalayan regions further north. We found a significant indication that, among Stenus, Palearctic species are slightly more represented in the upper half of the mountain range (> 1400 m) (and less in the lower half) than expected, and vice versa for the Oriental species. This pattern might be a remnant of climatic-niche conservatism that, in some of these species, still seems to determine their elevational distribution.

Whereas the elevational niche breadth appears to be mainly determined by the preferred temperature regime in many Stenus species, the availability of suitable habitat structures associated with waterfall environments seems to be the crucial factor in “Dianous” (see also previous section). Most “Dianous” species show relatively low habitat niche widths (< 0.3) corresponding to their close association with sprinkled rocks around cascades (cf. Fig. 4, Online Resource 3). This is also supported by our DA in which the sprinkled rocks in streams and waterfalls form the only habitat type distinguishable by the occurring species assemblages. Such narrow habitat selection might be the main reason for our failure to collect significant “Dianous” quantities above 1400 m, where this habitat type is absent. Hence, the lack of “Dianous” at higher elevations in both ITH and PHK is most probably not an effect of the lower temperatures, but of the absence of suitable habitats. Indeed, in higher mountain ranges of the Himalayas (e.g., Laos), “Dianous” occurs at elevations up to 3800 m. Similarly, the lack of aquatic Hemerodromiini (Diptera: Empididae) at elevational extremes of the ITH is considered to have resulted from the paucity of suitable streams (Plant et al. 2012). Hence, although all collected “Dianous” species belong to the Palearctic realm, they are mostly restricted to the lower elevations (< 1400 m) (cf. Fig. 4). If they colonized this region from northern (cooler) Himalayan regions (like Palearctic Stenus), then, in this case and unlike some Stenus species, no remnant of climatic-niche conservatism seems to have been preserved. “Dianous” have adapted to the extreme conditions of waterfall environments that are not well exploited by Stenus, possibly leading to a situation of competitive release. Under such an evolutionary regime, (competition-driven) stabilizing selection for narrow thermal niches might have become rapidly reduced, so that this taxon has easily dispersed across elevational ranges down the slope and settled successfully in this competition-free zone. The process might have been (1) encouraged by a source–sink (mass effect) scenario (cf. Grytnes and McCain 2013) in which adult beetles and their propagules are permanently transported by a water stream from upper (source) to lower (sink) elevations and (2) accompanied by, at first, suboptimal temperature regimes, although the waterfall and stream environments might have maintained an (azonal) cool-temperate microclimate. Because of the competition-free situation in the area, sink populations can finally successfully establish themselves at lower elevations.

Aim 4: Niche breadth as a predictor for range size

Niche breath-range size relationships have been previously reported across various taxonomical groups and provide a potential general mechanism to explain the rarity of some species, despite others being common (Brown 1984; Slatyer et al. 2013). As expected, we found a positive correlation between the elevational niche width and species frequency, i.e., generalists tend to show more successful colonization of habitats across the investigated mountain areas than specialists (cf. Fig. 5). Regarding the habitat niche, this is only true for Stenus, because the habitat niche of “Dianous” beetles is largely restricted to waterfall environments. Even within a certain elevational range, namely the lowest elevational ranges (400–900 m and 900–1400 m), Stenus species frequency was positively correlated with habitat niche breadth. Thus habitat generalists have the potential to become more successfully established within certain elevational zones than specialists.

Considering elevational niche width as a proxy for temperature tolerance, the observed positive relationships between species frequency and elevational niche width suggest, from an evolutionary perspective, that colonization and dispersal processes have been accompanied by various patterns of climate-niche evolution: (1) increased climatic (thermal) tolerances enable certain species to disperse across elevations (both up- and downwards), larger geographic areas, and environments (cf. S. virgula) and may thus facilitate speciation; (2) rare species with only narrow climate (thermal) niche widths (possibly associated with patterns of strong climatic-niche conservatism) (e.g., S. luteomaculatus) will largely remain restricted to a narrow range of elevations, leading to increased isolation, which also promotes (endemic) speciation. Such patterns presumably have superimposed the overall radiation of Steninae into defined ecomorphs associated with the four major habitat types (i.e., bare ground, litter, vegetation, waterfalls) (cf. Betz et al. 2018); they open an approach toward more detailed considerations of the local and regional phylogeographic patterns involved in Steninae evolution.