Distribution and environmental determinants of darkling beetles assemblages (Coleoptera: Tenebrionidae) in Península Valdés (Argentinean Patagonia)

Abstract

Arid lands provide several products and services to humankind, but human activities affect this environment, increasing the extinction risk of the native species. Thus, to successfully conserve the biodiversity of these ecosystems, it is necessary to identify which environmental factors influence the spatial distribution of the organisms that offer these benefits. Darkling beetles play a relevant role on the functioning of deserts. Although these insects are among the most abundant and diverse in these environments, there is no agreement on the relative importance that different environmental factors have as determinants of their spatial distribution. In this work, we assessed the role of climate, vegetation, and soil variables as factors that determine distribution patterns of darkling beetles within the Natural Protected Area Península Valdés (Northeastern Patagonia). Five groups of environmental units were identified, each one with an exclusive tenebrionid species assemblage and different species diversity. The most influential environmental variables were temperature, precipitations, and soil texture. Results suggest that the magnitude of several ecosystem processes may vary among the groups of environmental units. We recommend prioritizing the conservation of the five groups of environmental units and incorporating the darkling beetles-environment relationship in future conservation strategies for arid Patagonia.

Key words
habitat and climate variables; arid Patagonia; spatial patterns; tenebrionids

INTRODUCTION

Arid and semi-arid lands are spatially heterogeneous habitats that support high biodiversity with multiple biological interactions (Whitford 2002, Ayal 2007, Baldi et al. 2017). It is recognized that this biodiversity provides several products and services to humankind, from economic gains (e.g. sheep wool) to atmospheric and climatic regulation (Whitford 2002). However, some human activities such as overgrazing and habitat fragmentation are degrading those habitats leading to biodiversity declines (Whitford 2002, Laity 2008, Ward 2009, Baldi et al. 2017). Preventing the loss of biological diversity is critical for an effective management and conservation strategy of those ecosystems (Northrup et al. 2013). Consequently, it is strongly important to increase the understanding of the actual factors that promote and maintain biodiversity and ecosystem processes in arid environments (Mazía et al. 2006, Martínez et al. 2018).

It is known that environmental conditions act as ecological filters of regional species pools, determining species composition at a given time and place (Lichti & Murphy 2010). In arid ecosystems the distribution of most species is affected by a combination of ambient temperature and rainfalls because these abiotic factors influence animal ecophysiology (Cloudsley-Thompson 2001, Ward 2009, Schowalter 2016). However, the high spatial variation in geological substrates and soil types in deserts may also play a key role as spatial filters, especially for organisms whose life cycle is partly developed in the soil (Ward 2009). Although related to previous factors, vegetation may also be essential to desert animals because, besides offering shelter and buffering microclimatic variation, it harbors greater concentrations of water, soil nutrients, and food (litter or potential prey) (Mazía et al. 2006, Ward 2009, Schowalter 2016). In this context, to successfully maintain the products and services provided by the biodiversity of arid environments, it is necessary to identify which environmental factors determine the presence of the living organisms that perform these benefits and develop sustainable management strategies including relations between biodiversity and environmental factors (Whitford 2002, Mazía et al. 2006, Northrup et al. 2013, Martínez et al. 2018).

Insects represent the majority of species in desert communities and are an integral part of their structure and dynamics (Schowalter 2016). Darkling beetles (Coleoptera: Tenebrionidae) are among the most abundant and diverse biomass contributing invertebrates in deserts (Cloudsley-Thompson 2001, Cepeda-Pizarro et al. 2005, Cheli et al. 2010, Matthews et al. 2010, Baldi et al. 2017). They comprise about 20,000 species around the world (Matthews et al. 2010) and have a relevant role in the functioning of desert ecosystems (Ayal 2007, Pizarro-Araya 2010, Bartholomew & El Moghrabi 2018). A number of studies in deserts around the world have stressed the importance of darkling beetles as detritivores in xeric ecosystems, suggesting that they play an important role in the cycle of important soil nutrients (Matthews et al. 2010). These insects can behave either as generalist herbivores (Flores & Debandi 2004, Cheli et al. 2009, Bisigato et al. 2015) or as consumers of the senescent vegetation (Ayal 2007, Pizarro-Araya 2010, Cheli 2009) or dead animals (Aballay et al. 2016, Cheli 2009). On the other hand, darkling beetles are important prey items for numerous vertebrates, transferring energy from low to higher trophic levels (Ayal & Merkl 1994, Flores, 1998, Formoso et al. 2012). Despite the extensive researches that exist on the ecology of darkling beetles, there is no agreement on the relative importance of the different environmental factors determining their spatial distribution (Krasnov & Shenbrot 1996). Some studies have documented that the distribution of darkling beetles mainly respond to variations in soil texture (Sheldon & Rogers 1984, Crawford 1988, Ayal & Merkl 1994, Matthews et al. 2010). Other works claim that vegetation cover and complexity are the main determinants of the darkling beetles diversity (Cepeda-Pizarro 1989, Parmenter et al. 1989a, Mazía et al. 2006, Liu et al. 2012). In contrast, other authors consider that climatic factors like temperature and precipitation/humidity may have the greatest influence on the spatial distribution of darkling beetles (Parmenter et al. 1989b, Flores 1998, Cloudsley-Thompson 2001, Carrara et al. 2011a, Rosas et al. 2019).

The Natural Protected Area Península Valdés (PV from now onwards) is situated at the northeastern portion of Patagonia (42º05’–42º53’S; 63º35’–65º04’W) in the Atlantic coast of Chubut, Argentina. This area, a UNESCO Natural World Heritage Site and Biosphere Reserve and one of the biggest arid protected areas of Argentina, has environmental singularities in climate (Coronato et al. 2017), geomorphology (Bouza et al. 2017a, b), and soils (Rostagno 1981, Bouza et al. 2017a, b) that determine several environmental units within the peninsula (Bertiller et al. 2017) which host a rich biodiversity (Baldi et al. 2017). Knowledge of terrestrial arthropods in PV has been considerably increased during the last decade (e.g. Cheli 2009, Cheli et al. 2010, 2013, Carrara et al. 2011b, Flores et al. 2011, Cheli & Martínez 2017, Baldi et al. 2017, Martínez et al. 2018). Previous studies of darkling beetles diversity in PV have focused on community composition and darkling beetles distribution patterns (Cheli et al. 2013, Flores et al. 2011, Carrara et al. 2011b), but relations among darkling beetles with climate, soil, and vegetation remain almost unknown. The purpose of this study was to analyze the spatial variations of tenebrionid assemblages in PV and to identify the main environmental factors that determine it. For this reason, the following questions will be addressed: Do species composition of tenebrionid assemblages and diversity vary among environmental units in PV? If so, which are the environmental variables that mainly determine this variation? Which is the relative contribution of soil, vegetation, and climate variables to the variation of tenebrionid species assemblages within PV?

METHODS

Study area

The study was conducted in the Natural Protected Area Península Valdés (PV), a wide plateau of low altitude ranging between 35 m below sea level because of a series of central salt flats and 80 m a.s.l. that covers approximately 4000 km² (Figure 1). Its climate is arid, characterized by hot, dry summers and relatively cold winters, with a mean annual temperature of 12ºC, and winds prevailing from the west (Coronato et al. 2017). Geomorphology, climate, and soil features show great variability inside the peninsula (Rostagno 1981, Coronato et al. 2017, Bouza et al. 2017a, b); as the result of this variability it was possible to recognize 16 environmental units within (Bertiller et al. 2017) (Table I).

Figure 1
Location of the study area: the Natural Protected Area Peninsula Valdés, Patagonia, Chubut, Argentina.

Table I
Environmental units of Peninsula Valdes where tenebrionid assemblages were studied, detailing the number of sampling sites in each. The dominant plant formation in each one is described on the right column. * The numbering of the environmental units in this table is the original number assigned to these units by Bertiller et al. (2017). The authors preferred not to vary this numbering so that interested readers can more easily find these units in the reference. Environmental units located outside the peninsula or on its isthmus (see Bertiller et al. 2017) were not considered in this study.

Environmental Unit*	Main Vegetation Formation	Number of sampling points
1	Perennial grass steppe of Sporobolus rigens and Nassella tenuis	5
2	Perennial grass steppe of Piptochaetium napostaense, Nassella tenuis and Plantago patagonica	7
3	Perennial grass steppe of Nassella tenuis, Nassella longiglumis with shrubs of Chuquiraga avellanedae	3
4	Shrub-perennial grass steppe of Chuquiraga erinacea subsp. erinacea and Nassella tenuis	7
5	Shrub-perennial grass steppe of Chuquiraga avellanedae and Nassella tenuis	18
7	Shrub steppe of Condalia microphylla and Lycium spp	3
8	Shrub steppe of Chuquiraga avellanedae and Chuquiraga erinacea subsp. erinacea	4
9	Shrub steppe of Chuquiraga avellanedae and Condalia microphylla	19
10	Shrub steppe of Schinus johnstonii and Lycium chilense	3
11	Shrub steppe of Chuquiraga avellanedae and Mulinum spinosum	12
12	Shrub steppe of Senecio filaginoides and Mulinum spinosum	10
13	Shrub steppe of Chuquiraga erinacea subsp. hystrix and Chuquiraga avellanedae	10
15	Shrub steppe of Cyclolepis genistoides, Chuquiraga avellanedae and Atriplex lampa	6
16	Shrub steppe of Chuquiraga erinacea subsp. Hystrix, Cyclolepis genistoides and Chuquiraga avellanedae whit perennial grasses	6
20	Mosaic: Perennial grass steppe of Sporobolus rigens and Nassella tenuis (1) and Dwarf shrub steppe de Hyalis argentea (6)	9
21	Mosaic: Perennial grass steppe of Piptochaetium napostaense, Nassella tenuis and Plantago patagonica (2) Shrub-perennial grass steppe of Chuquiraga avellanedae and Nassella tenuis (5)	12

Data

Darkling beetle data

Darkling beetle records inside PV were compiled from 86 sampling points of previous works in the study area (Cheli et al. 2010, 2013, Flores et al. 2011, Carrara et al. 2011b, Martínez et al. 2018), plus a trip collecting to PV carried out in January 2010 that comprised 48 sites in the 16 terrestrial environmental units within PV, i.e. three sites per environmental unit. Thus, a total of 134 sampling points (86 + 48 points) have been included. In each site, a 20-minutes visual inspection was conducted in an approximately 0.5 ha by three observers who scrutinized the ground, shrubs, rocks, or any other object on the ground that could offer shelter for darkling beetles. All individuals (live or dead) observed during were collected (more details about the collecting process can be found in Cheli et al. 2016a). To increase sampling representativeness, surveys in each environmental unit were performed on sites located at least 10 km from each other. In order to avoid possible bias, analyses were performed only with those darkling beetle species with clearly epigeal habits. Individuals were identified to species level based on reviews and keys by Kulzer (1955, 1963), Flores (1999, 2004), Flores et al. (2011) and comparisons with material deposited in the IADIZA collection. We followed the classification proposed by Matthews et al. (2010) and Kamiński et al. (2021) to assign species and genera into subfamilies and tribes. Following the recent taxonomic changes suggested by Silvestro & Flores (2016), Nyctelia nodosa (Germar, 1823) mentioned in previous works (Cheli et al. 2010, 2013, Carrara et al. 2011b) was treated here as Nyctelia picipes (Billberg, 1815). After a recent study of some Kulzer & Kaszab`s type specimens at the Hungarian Natural History Museum in Budapest, Hungary, we present an improved list of darkling beetles species from Peninsula Valdés (Table II). Here we refer to Hylithus kovaksi Kaszab and Psectrascelis hirtus Kulzer instead of Hylithus tentyroides (Lacordaire) and Psectrascelis sulcicollis (Waterhouse) listed in our previous work (Carrara et al. 2011b). The tenebrionid specimens collected were deposited in the Entomological Collections of IPEEC and IADIZA. Each tenebrionid record was assigned to an environmental unit (Bertiller et al. 2017) from georeferenced maps and satellite images of PV using QGIS software (QGIS Development Team 2018). Because tenebrionid records in PV did not came from a long-term standardized sampling, we decided to use incidence data (species presence/absence). With the complied records of darkling beetles species, from both previous works (Cheli et al. 2010, 2013, Flores et al. 2011, Carrara et al. 2011b, Martínez et al. 2018) and the collecting campaign, a matrix of species incidence by environmental unit was constructed (Table II) to perform all statistical analyses.

Table II
Tenebrionid species (and their supraspecific taxonomic classification) recorded for the environmental units within Peninsula Valdes. The numbering of the environmental units is the same as that in Table I. * Species acronym used in Figure 2. Environmental units: 0 = species absent; X = present.

Subfamily	Tribe	Species	Acronym*	Envionmental unit
				1	2	3	4	5	7	8	9	10	11	12	13	15	16	20	21
Pimeliinae	Edrotini	Hylithus kovaksi	Hyli_kova	X	X	X	X	X	0	X	X	0	X	0	X	0	0	X	X
Pimeliinae	Nycteliini	Epipedonota cristallisata	Epip_cris	X	X	X	X	X	X	X	X	0	X	X	X	X	X	X	X
Pimeliinae	Nycteliini	Mitragenius araneiformis	Mitra_aran	0	X	X	0	X	0	X	X	X	X	0	X	0	X	X	X
Pimeliinae	Nycteliini	Nyctelia circumundata	Nycte_circ	X	X	0	X	0	0	0	0	0	X	X	0	X	0	X	X
Pimeliinae	Nycteliini	Nyctelia darwini	Nycte_darw	0	0	0	0	X	0	0	0	0	0	0	0	0	0	0	0
Pimeliinae	Nycteliini	Nyctelia dorsata	Nycte_dors	X	X	0	X	0	0	X	0	0	X	X	X	X	X	X	X
Pimeliinae	Nycteliini	Nyctelia picipes	Nycte_pici	0	X	X	X	X	X	0	X	X	X	0	X	0	0	X	X
Pimeliinae	Nycteliini	Patagonogenius collaris	Pata_colla	0	0	0	X	0	0	0	0	0	X	X	X	X	X	0	0
Pimeliinae	Nycteliini	Patagonogenius quadricollis	Pata_quad	0	0	0	0	0	0	0	0	0	X	0	0	X	0	0	0
Pimeliinae	Nycteliini	Psectrascelis hirtus	Psec_hirt	0	0	0	0	0	0	0	0	0	X	X	X	X	X	0	0
Pimeliinae	Physogasterini	Pimelosomus sphaericus	Pime_spha	0	0	0	0	0	0	0	0	0	0	0	0	0	X	0	0
Pimeliinae	Praocini	Calymmophorus patagonicus	Caly_pata	0	0	0	0	0	0	0	0	0	0	0	0	X	X	0	0
Pimeliinae	Praocini	Calymmophorus peninsularis	Caly_peni	0	0	0	0	0	0	0	0	0	0	0	0	X	0	0	0
Pimeliinae	Praocini	Plathestes kuscheli	Plat_kusch	X	X	0	0	X	0	0	0	X	X	X	X	X	X	0	0
Pimeliinae	Praocini	Praocis (Hemipraocis) sellata granulipennis	Prao_sgra	X	0	0	X	0	X	0	0	0	0	0	X	X	0	0	X
Pimeliinae	Praocini	Praocis (Hemipraocis) sellata peninsularis	Prao_spe	X	0	0	0	X	0	0	0	0	0	X	X	X	0	0	0
Pimeliinae	Praocini	Praocis (Hemipraocis) sp.1	Prao_sp	0	0	0	0	0	0	0	0	0	0	0	X	0	0	0	0
Pimeliinae	Praocini	Praocis (Orthogonoderes) argentina	Prao_arge	0	0	0	0	0	0	0	0	0	X	0	0	0	0	0	0
Pimeliinae	Stenosini	Ecnomoderes bruchi	Ecno_bruc	X	0	0	0	0	0	0	X	0	0	0	0	0	0	X	X
Pimeliinae	Trilobocarini	Salax lacordairei	Sala_laco	X	0	0	0	0	0	0	0	0	0	X	0	0	0	0	0
Blaptinae	Opatrini	Blapstinus punctulatus	Blap_pun	X	X	X	X	X	0	X	X	0	X	0	0	0	0	X	X
Tenebrioninae	Scotobiini	Emmallodera crenaticostata crenaticostata	Emma_cren	0	0	0	0	X	0	0	0	0	X	X	0	X	0	0	0
Tenebrioninae	Scotobiini	Emmallodera hirtipes	Emma_hirt	X	X	X	X	X	X	X	X	0	X	X	X	X	X	X	X
Tenebrioninae	Scotobiini	Leptynoderes nordenskioldi	Lepty_nord	X	0	0	X	0	0	0	X	0	0	0	X	0	0	X	X
Tenebrioninae	Scotobiini	Leptynoderes tuberculata	Lepty_tube	0	0	0	0	0	0	0	0	0	0	0	X	X	0	X	X

Environmental variables

Habitat variables were compiled from three sources: a) specific bibliography: Rostagno 1981, Bertiller et al. 2017, Bouza et al. 2017a, b; b) unpublished environmental information kindly offered by Drs. M. Rostagno and P. Bouza (IPEEC, CCT CENPAT CONICET, pers. com.); and c) satellite images of Moderate Resolution Imaging Spectroradiometer (MODIS). Compiled habitat information was classified into three classes: 1- Climatic variables: Precipitation; Diurnal mean Temperature; Nocturnal mean temperature. 2- Soil variables: Dominant type soil; Geomorphological system (Uplands and plains, Great endorheic basins and Coastal zone); Geomorphological Subsystem (Terrace Levels, Aeolian fields, Piedemont Pediments and Bajadas, Coastal Piedemont Pediments, Pleistocene beach ridges, Holocene beach ridges); Number of soil horizons; Depth of “A” soil horizon (topsoil); Dominant topographic Slope; pH; Texture of the most superficial 20 cm of the soil profile (% gravel, % sand, % clay and % silt); Soil organic carbon content. 3- Vegetation variables: Normalized Difference Vegetation Index (NDVI); Enhanced Vegetation Index (EVI); Actual Evapotranspiration and Potential Evapotranspiration; Functional type of steppe (perennial grass steppe, dwarf shrub steppe, shrub steppe, perennial grass-shrub steppe); Number of vegetation strata (grassy, Dwarf-shrubs, low shrubby (<50 cm), tall shrubby (>50 cm) strata); Total vegetation cover (%); Vegetation maximums height (cm); Plant Species Richness (Table III).

Table III
Habitat variables recorded for the environmental units within Peninsula Valdes. * The numbering of the environmental units is the same as that in table I. References: Dominant Slope: 1: flat slope (2-2%), with slow runoff // 2: gently slope (2-4%), with moderate runoff // 3: inclined (4-10%), with fast runoff. Soil Type (Soil Taxonomy): 1: Natrargids Xerollic (BT clay horizon) // 2:Torripsamment Xeric (O) + Calciorthids Xerollic (Sin BT) // 3: Torripsamment Xeric + Natrargids Xerollic // 4: Torripsamment Tipic (saline horizon) // 5: Torripsamment Tipic + Salorthids Aquollic (saline soils) // 6: Loose sand blanket + Torripsamment Xeric // 7: Rodados Patagonicos (gravels) without fine material + Natrargids Xerollic. Steppe Functional type: 1: Perennial grass steppe // 2 Dwarf-shrub steppe // 3 Perennial grass-shrub steppes // 4 Shrub-perennial grass steppes // 5 Shrub steppe. Geomorphological System: 1: System A: uplands and plains // 2: System B: great endhoreic basins // 3: System C: coastal zone. Geomorphological Subsystem: 1: Terrace Levels // 2: Aeolian fields // 3: Piedemont Pediments and Bajadas // 4: Coastal Piedemont Pediments // 5: Pleistocene beach ridges // 6: Holocene beach ridges.

Environmental unit *	NDVI	EVI	Precipitation (mm)	Actual Evapotranspiration	Pot. Evapotranspiration	Temp. Day (K°)	Temp .Night (K°)	Steppe functional type	Veg. strata	Soil type	Geomorph. system	Geomoph. subsystem	Dominant Slope	Veg. cover (%)	Veg. height (cm)	Veg. Species Richness	Depth soil horizon A (cm)	Num. Soil horizons	Gravels (%)	Sand (%)	Silt (%)	Clay (%)	Soil Organic Carbon	PH
1	16762400	10454000	11.461	8.9	256.192	307.672	284.769	1	1	3	1	2	1	75	30	6	20	1	5	92.4	6.6	1	0.32	7.72
2	22585428.57	13932571.43	15.799	14.306	249.806	307.619	284.129	1	2	2	1	1	1	65	20	13	10.67	2	3.37	73.48	13.18	13.3	0.723	8.057
3	18661333.33	9653333.333	14.209	12.878	265.556	308.411	284.376	1	3	7	3	5	1	85	30	11	10	2.7	9.67	56.78	20.56	22.66	0.831	8.154
4	19902000	12112857.14	9.959	12.622	265.933	305.397	284.823	4	3	4	3	4	2	55	80	11	13.33	1.7	14.17	72.4	15.8	11.8	0.702	7.845
5	20935777.78	12088777.78	15.212	8.748	257.644	308.637	284.8	4	2	1	1	1	1	50	55	10	8.5	2	2.53	71.51	12.41	16.04	0.74	8.083
7	20429000	12210000	14.209	12.167	257.2	306.093	284.57	5	2	4	3	4	2	60	90	12	13.33	1.7	14.17	72.4	15.8	11.8	0.702	7.845
8	21356000	12610000	10.849	11.408	261.242	309.767	284.582	5	2	4	1	1	2	70	75	19	8.5	2	2.53	71.51	12.41	16.04	0.74	8.083
9	21174800	12255900	11.775	10.075	261.813	309.142	285.07	5	2	1	1	1	1	55	80	12	8.5	2	2.53	71.51	12.41	16.04	0.74	8.083
10	14370666.67	8984666.667	19.632	12.711	263.411	307.267	283.911	5	2	7	3	6	2	40	75	8	1	1	73	26.7	28.2	45.1	0.22	7.98
11	16399833.33	9557500	13.803	9.093	272.12	305.226	284.175	5	3	2	3	4	1	50	95	12	13.33	1.7	14.17	72.4	15.8	11.8	0.702	7.845
12	19835600	12316200	28.476	7.753	267.77	307.367	285.48	5	3	6	3	4	3	50	90	11	13.33	1.7	14.17	72.4	15.8	11.8	0.702	7.845
13	17857200	10988400	12.873	8.17	257.61	306.399	284.833	5	2	5	2	3	3	70	115	20	10.33	2.3	6.17	70.2	22.1	7.7	0.65	7.643
15	24952305.56	12430555.56	25.48	8.061	273.806	303.78	285.963	5	3	4	3	4	3	60	80	9	13.33	1.7	14.17	72.4	15.8	11.8	0.702	7.845
16	18835333.33	11119083.33	19.118	9.461	263.922	307.39	285.152	5	3	4	2	3	3	50	80	11	17	2	10	69.8	18.2	12	0.71	7.43
20	19876666.67	10956222.22	9.912	8.404	260.681	308.071	284.633	2	2	3	1	2	1	80	50	9	20	1	5	92.4	6.6	1	0.32	7.72
21	20474000	11470000	9.245	9.37	260.978	307.721	284.435	3	2	2	1	1	1	65	20	16	10.67	2	3.37	73.48	13.18	13.3	0.723	8.057

Statistical data analysis

We applied non-metric multidimensional scaling (NMDS) to visualize the variation in composition of the darkling beetle community among environmental units within PV. The purpose of this technique is to perform an ordination of the samples in function of their species similarity (Clarke & Warwick 2001, Legendre & Legendre 2012). We performed the analysis based on a matrix of biological similarity, using the Sørensen index as a measure of distance on the taxa incidences (Legendre & Legendre 2012). Differences among the groups evidenced by the NMDS ordination were tested using a one-way analysis of similarity (ANOSIM) (Clarke & Warwick 2001, Legendre & Legendre 2012). This test allows comparing assemblages as function of factors based on distance measurements. ANOSIM significance was determined using 999 random permutations of group membership.

To identify tenebrionid species that characterize variation in beetles’ assemblages among environmental units, a Principal Components Analysis (PCA) was performed (Borcard et al. 2011, Legendre & Legendre 2012). It is known that PCA can be affected by double-zero cases (Ruokolainen & Blanchet 2014); to avoid this inconvenience, data were transformed using Hellinger distance (Borcard et al. 2011, Legendre & Legendre 2012, Ruokolainen & Blanchet 2014). We evaluated how many PCA axes were significant using the broken-stick distribution. PCA axes with larger percentages of explained variance than the broken-stick model were considered significant (Borcard et al. 2011, Legendre & Legendre 2012).

In order to relate habitat variables with the structure of the darkling beetle community, a Redundancy Analysis (RDA) was performed. Environmental variables used in RDA analyses were previously standardized. To avoid multicollinearity, those variables with correlations higher than r = 0.7 were previously excluded from analyses. To assess the relative importance of soil, vegetation, and climatic variables on changes in the structure of the darkling beetle community, a variation partitioning analysis (VARPART) was performed via RDA (Borcard et al. 2011, Legendre & Legendre 2012). Significance of global RDA analysis, individual axes, and for each testable fraction of the variation partitioning were evaluated using unrestricted Monte Carlo permutation tests with 499 permutations (Borcard et al. 2011, Legendre & Legendre 2012).

Following Chao et al. (2014) and Hsieh et al. (2016), darkling beetle diversity was studied using integrated sample-size and coverage-based rarefaction/extrapolation sampling curves with Hill numbers based on incidence data (qΔ). Hill numbers for incidence data are based on the relative species incidence in the sites that constitute the studied assemblages (Chao et al. 2014). To incorporate the full effect of relative species incidences on diversity estimation, curves were plotted for the Hill numbers q = 0, q =1 and q = 2. The parameter q determines the sensitivity of qΔ to the relative species incidences. When q = 0, the incidence of individual species is not considered, so that 0Δ indicates simply species richness. If q = 1, species are weighted in proportion to their incidence (more weight on “typical” species) and it is homologous to the exponential of the Shannon index (Shannon diversity). When q = 2, the estimated number of effective species is similar to the inverse of the Simpson index (Simpson diversity) giving more weight to species widely present in the community (dominant species) (Chao et al. 2014). A diversity profile (a plot of qΔ vs. q from q = 0 to =3) was also performed (Chao et al. 2014, Hsieh et al. 2016). The slope of the diversity profile curve reflects the unevenness of species relative incidences. The more uneven the distribution of relative species incidences, the more steeply the curve declines (Chao et al. 2014, Hsieh et al. 2016). The relationship between sample coverage and sample size was studied using a sample completeness curve for both smaller rarefied samples and larger extrapolated samples (Chao et al. 2014, Hsieh et al. 2016). All curves were plotted with their 95% confidence intervals. If their intervals do not overlap, the curves are statistically different (Chao et al. 2014, Hsieh et al. 2016).

We performed NMDS, ANOSIM, PCA, RDA and VARPART analyses using metaMDS, anosim, rda and varpart functions of the vegan package (Oksanen 2019) for R (R Core Team 2018). We evaluated how many PCA axes were significant with the evplot function for R (Borcard et al. 2011). Diversity analyses were performed using the iNEXT package for R (Hsieh et al. 2016).

RESULTS

Taxonomic Composition

The 134 sampling points compiled into the 16 environmental units within PV, contained 378 records of tenebrionid beetles. As a result, the presence of 25 species of tenebrionid beetles belonging to three subfamilies and eight tribes was found inside PV (Table II). Forty-four percent of the species (11 species) were registered in up to 25% of the environmental units. Five species (20% of the species) were recorded exclusively in one environmental unit: Calymmophorus peninsularis Flores & Cheli (in shrub steppes of Cyclolepis genistoides, Chuquiraga avellanedae and Atriplex lampa), Nyctelia darwini Waterhouse (in shrub-herbaceous steppes of Chuquiraga avellanedae and Nassella tenuis), Pimelosomus sphaericus Burmeister (in shrub steppes of Chuquiraga avellanedae, Cyclolepis genistoides and Chuquiraga erinacea subsp. hystrix), Praocis (Orthogonoderes) argentina Kulzer (in shrub steppes of Chuquiraga avellanedae and Mulinum spinosum) and Praocis (Hemipraocis) sp.1 (in shrub steppes of Chuquiraga erinacea subsp. hystrix and Chuquiraga avellanedae). Twenty-four percent of species were present in 26-50% of the environmental units, and an equal percentage in 51-75%. Only two species (8%) were widely distributed among environmental units (between 76-100%): Emmallodera hirtipes Kulzer and Epipedonota cristallisata (Lacordaire), both of which were present in all environmental units, except for the shrub steppe of Schinus johnstonii and Lycium chilense (Table II).

Variation in the tenebrionid species assemblages

The NMDS evidenced that the specific composition of the darkling beetle community varies among different environmental units within PV, consequently five tenebrionid species assemblages can be identified (Figure 2). The first was present in a group of environmental units conformed by units 1, 4, 20, 21 (G1); the second in units 11, 12, 13, 15, and 16 (G2); units 2, 3, 5, 8 and 9 conform the third (G3); species assemblage G4 was only found in environmental unit 7 (shrub steppes of C. microphilla and Lycium sp.), and G5 in unit 10 (shrub steppes of S. johnstonii and Lycium chilense). Figure 2 shows a distant location of environmental units 7 and 10 with respect to the rest of the members of the other groups (G1, G2 and G3), evidencing that both units have extremely dissimilar tenebrionid compositions. Analysis of similarities confirmed that the composition of darkling beetles species varied significantly among the five tenebrionid species assemblages (ANOSIM: R = 0.76; p = 0.001; 999 permutations, Table IV).

Figure 2
NMDS diagram showing the ordination of the environmental units within PV according to their tenebrionid species composition. Polygons represent each of the three main tenebrionid species assemblages identified by ANOSIM (Group 1 = triangles; Group 2 = squares; Group 3 = circles) while assemblages present in only one environmental unit are represented by a hexagon (Group 4) and a diamond (Group 5).

Table IV
Differences in similarity among tenebrionid species assemblages present in the five groups of environmental units analyzed with ANOSIM. Note: * Groups 4 and 5 were excluded from pairwise test because performing multiple comparisons ANOSIM requires the groups to be composed by more than one element in size (Legendre & Legendre 1998, Clarke & Warwick 2001).

Groups of environmental units	Global Test	Pairwise Tests*
		Groups	R Statistic	P
G1 (11, 12, 13, 15,16)		1 vs 3	0.812	0.008
G2 (1, 4, 20, 21)	Global R = 0.756	1 vs 2	0.734	0.016
G3 (2, 3, 5, 8, 9)	( p = 0.001)
G4 (7)		3 vs 2	0.541	0.008
G5 (10)

The first two PCA axes explained together almost 50% of the variation in the tenebrionid species composition, evidencing the same ordination of environmental units as NMDS did (Figure 3a). The broken-stick distribution model showed that only these two axes explain a significant proportion of the variation in the darkling beetle species composition among environmental units within PV (Figure 3b). The first axis, explaining 29.3% of the variability in darkling beetle species composition, essentially distinguished group 2 (clumped on the right of this axis), from the other environmental units (located on the left of the ordination diagram). This group of environmental units (G2) had a tenebrionid assemblage mainly characterized by the presence of Emmallodera crenaticostata crenaticostata Blanchard, Psectrascelis hirtus, Patagonogenius collaris (Kulzer), Praocis (Hemipraocis) sp. 1, and to a lesser extent Salax lacordairei Guérin-Méneville, Calymmophorus peninsularis, Calymmophorus patagonicus Bruch, Patagonogenius quadricollis (Fairmaire), Pimelosomus sphaericus, and Plathestes kuscheli Kulzer (Figure 3a). The second PCA axis, accounting for 19.9% of the darkling beetle variability, mainly distinguished G3 and G5 (in their negative portion) from G1 (located in the positive sector of this axis). G3 was mainly characterized by the presence of Nyctelia picipes, Mitragenius araneiformis Curtis and to a lesser extent Plathestes kuscheli. G1 typically Leptynoderes nordenskioldi Kulzer, Leptynoderes tuberculata Curtis, Nyctelia circumundata Lesne, Nyctelia dorsata Fairmaire, Ecnomoderes bruchi Gebien occurred, and to a lesser extent Hylithus kovaksi and Blapstinus punctulatus Solier (Figure 3a). The composition of tenebrionid assemblage in the environmental unit 10 (G5) was similar to that of G3 but with singular absence of Emmallodera hirtipes and Epipedonota cristallisata (both species were present in all environmental units, except for this environmental unit). In addition, the tenebrionid species assemblage in environmental unit 7 (G4) was similar to that in G2 but whit the notable absence of B. punctulatus.

Figure 3
A. Ordination of the environmental units as a function to their tenebrionid species composition in PV on the plain defined by the first two axes of PCA. The tenebrionid species assemblages are represented by different symbols (Group 1 = triangles; Group 2 = squares; Group 3 = circles; Group 4 = hexagon; Group 5 = diamond). Species acronyms as in Table II. B. Significance of PCA axes applying the broken-stick distribution. PCA axes with larger percentages of variance than the broken-stick variances were considered significant.

Habitat characteristics responsible for the variation in the composition of tenebrionid species assemblage

RDA analyses evidenced that the combination of environmental variables that explained the greatest variability in tenebrionid species composition within PV was constituted by: pH, Diurnal mean Temperature, Number of soil horizons, Actual Evapotranspiration, % clay, Precipitation, Nocturnal mean temperature, Soil organic carbon content, and NDVI. These variables explained 77.2% of the variation in the tenebrionid species composition and showed a significant relationship for the overall test on all constrained axes (Pseudo F = 2.2626; p = 0.001; with 999 permutations). The ordination space defined by the first two RDA axes (Figure 4) explained 45.7% of the total variation in darkling beetle species composition and evidenced the same ordination of environment units and darkling beetle assemblages than the first two axes of PCA (Figure 3a).

Figure 4
Ordination of the tenebrionid species (continuous lines), habitat variables (dotted lines) and environmental units (triangles, squares, and circles) on the plain defined by the first two axes of RDA. The tenebrionid species assemblages are represented by different symbols (Group 1 = triangles; Group 2 = squares; Group 3 = circles; Group 4 = hexagon; Group 5 = diamond). Numbers designate species: 1. Praocis (Hemipraocis) sellata granulipennis // 2. Nyctelia dorsata // 3. Nyctelia circumundata // 4. Leptynoderes tuberculata // 5. Salax lacordairei // 6. Praocis (Hemipraocis) sellata peninsularis // 7. Patagonogenius quadricollis // 8. Patagonogenius collaris // 9. Psectrascelis hirtus // 10. Calymmophorus patagonicus // 11. Emmallodera crenaticostata crenaticostata // 12. Plathestes kuscheli // 13. Mitragenius araneiformis // 14. Nyctelia picipes // 15. Blapstinus punctulatus // 16. Hylithus kovaksi // 17. Ecnomoderes bruchi // 18. Leptynoderes nordenskioldi // 19. Epipedonota cristallisata // 20. Emmallodera hirtipes. Variables acroyms: SOC: Soil organic carbon content // NT: Nocturnal mean temperature // DT: Diurnal mean Temperature // Precip: Precipitation // R.Evap: Real Evapotranspiration // NumHor: Number of soil horizons // NDVI: Normalized Difference Vegetation Index // Clay: % of clay in soil texture // PH: soil pH.

The environmental relationship assessed by the first RDA axis explained 28% of the variation in the tenebrionid species composition (Pseudo F = 7.3806; p = 0.001; 999 permutations) showing that the tenebrionid assemblage present in group G2 was positively correlated with nocturnal mean temperature and precipitation, and negatively with soil pH, diurnal mean temperature, number of soil horizons, and actual evapotranspiration (Figure 4). At the same time, the second RDA axis explained 17.7% of the variation in the composition of tenebrionid assemblages (Pseudo F = 4.6689; p = 0.027; 999 permutations) and evidenced that tenebrionid species assemblages in G3 and G5 were positively correlated to clay percentage in soil and real evapotranspiration; while those in G1 were related with greater NDVI and soil organic carbon content (Figure 4). Finally, tenebrionid species assemblage in G4 responded to intermediate environmental characteristics between G1 and G3.

The Variation Partitioning Analysis showed that climate and soil features were the greatest determinants of the tenebrionid variability within PV, explaining 15% and 14% respectively; while vegetation explained only 1% (Figure 5). The covariation (interaction) between climate and soil variables added 9% to the explained variability in the darkling beetle species composition, while the triple-interaction climate-soil-vegetation added only 6 % more (Figure 5).

Figure 5
Variation partitioning of RDA for the relationship among habitat variables and composition of tenebrionid species assemblages in PV. Three sets of explanatory variables were included: climate, soil and vegetation variables. The percentages corresponded to the variability in tenebrionid assemblage composition explained by each set of explanatory variables and their interactions.

Tenebrionid species diversity variation among environmental units of PV

The sample coverage for the recorded darkling beetle fauna of PV was 0.97 (CV: 0.707), indicating that sampling was complete, and therefore, the tenebrionid inventory is reliable. Figure 6 shows that 95% confidence intervals of the empirical and theoretical diversity curves overlapped entirely in the profile. Consequently, a reliable darkling beetle species diversity is estimated for the sampled community. At the same time, the slope of the estimated diversity profile reflected a moderate unevenness of darkling beetle species incidences in PV (Figure 6).

Figure 6
Tenebrionid diversity profiles curves plotting Hill numbers ^qΔ (∞) as a function of order q, 0 ≤ q ≤ 3. The slope of the diversity profile curves reflects the unevenness of species relative incidences. When the more uneven the distribution of relative species incidences is, the more steeply the curve declines (Chao et al. 2014, Hsieh et al. 2016). The continuous line represents the theoretical diversity profile, while the dotted line is the empirical diversity profile. Both curves were plotted with their 95% confidence intervals.

The sample coverage for assemblages in G1, G2 and G3 were respectively 93.33%, 92.35%, and 86.49 (Figure 7). Additionally, this curve showed that 95% confidence bands for the three groups widely overlapped among them, indicating that sampling completeness is nearly similar for the three darkling beetle assemblages (Figure 7). For G1 the effective number of species in the assemblage for q = 0, q = 1 and q = 2 were respectively 16, 14.3 and 13.4; for G2: 23, 19.3 and 17; while for G3: 14, 11.3 and 10 (solid triangles, circles and squares in Figure 8).

Figure 7
Sample completeness curve. The three main tenebrionid species assemblages identified by ANOSIM are represented by different symbols (Group 1 = triangles; Group 2 = squares; Group 3 = circles). It was not possible to include sample completeness curves for G4 and G5 groups because statistical analyses require groups of environmental units to be composed by more than one element in size (Chao et al. 2014, Hsieh et al. 2016). All curves were plotted with their 95% confidence intervals. Continuous lines represent the interpolated and dotted lines the extrapolated portions of the curve.

Figure 8
Sample-size and coverage-based rarefaction/extrapolation sampling curves with Hill numbers based on incidence data (q = 0, 1, 2) among the three main tenebrionid species assemblages identified by ANOSIM (Group 1 = triangles; Group 2 = squares; Group 3 = circles). It was not possible to include sample-size and coverage-based rarefaction/extrapolation sampling curves for G4 and G5 groups because statistical analyses require groups of environmental units to be composed by more than one element in size (Chao et al. 2014, Hsieh et al. 2016). All curves were plotted with their 95% confidence intervals. Continuous lines represent the interpolated and dotted lines the extrapolated portions of the curve. Upper: Sample-size-based rarefaction/extrapolation sampling curves. Bottom: Coverage-based rarefaction/extrapolation sampling curves.

The sample-size-based curves of the G1, G2 and G3 showed that 95% confidence bands of darkling beetle richness (q = 0) widely overlap among the three assemblages, indicating that species richness is nearly similar for all (Figure 8a). The estimate of the effective number of species in the assemblages considering their incidences (q = 1), showed that almost all 95% confidence intervals of G2 do no overlap (except for very small sizes) with those of G1 and G3 (Figure 8b). Consequently, diversity of this assemblage is significantly greater than the other two. At the same time, when q = 2, in G3 it is lower than G1 and G2 (Figure 8c). From the comparison among coverage-based curves, the number of effective species estimated for q = 0, 1, and 2 show similar ordering of diversities as in the sample-size-based curve (Figure 8d, 8e, 8f). However, these coverage-based curves clearly show that G2 presented greater darkling beetle species diversity than the other two assemblages, since their 95% confidence bands never overlapped when q = 1 and almost never when q = 2 (Figure 8e, 8f). Briefly, from both sample-size-based and coverage-based sampling curves, darkling beetle species richness (q = 0) did not differ among groups of environmental units, while the estimated number of effective species for q = 1, and 2 decreased from G2 > G1 > G3 (Figure 8). It was not possible to statistically compare the tenebrionid species diversity of assemblages present in G4 and G5 with the remaining groups of environmental units because analyses require groups of environmental units to be composed by more than one element in size (Chao et al. 2014, Hsieh et al. 2016). However, table II showed that the number of registered tenebrionid species in both environmental units was lower than in the other units (4 and 3 respectively).

DISCUSSION

This work evidenced that darkling beetles make a discretionary use of habitat types based on the environmental factors present in PV. Variations in assemblage composition patterns among environmental units were mainly determined by the variation in climatic and edaphic features, and to a lesser extent by vegetation characteristics. These findings confirmed those of Carrara et al. (2011b), who suggested that tenebrionids are distributed in discrete units within the peninsula. At the same time, our findings about environmental factors that determine distribution patterns of darkling beetles in PV are in agreement with those evidenced by other epigeal beetle communities both in northwestern (Mazía et al. 2006, Werenkraut & Ruggiero 2012) and southern Patagonia (Rosas et al. 2019). Sheldon & Rogers (1984), in a study at a similar scale carried out in Arizona, USA, found that the majority of darkling beetle species occurred in all or most of the environmental units. However, most of the tenebrionid species inhabiting the peninsula showed narrow distributions. This fact evidences their relatively constricted environmental preferences, which is especially evident in the five species that were found only in one environmental unit. In the present study few tenebrionid species were widely distributed among environmental units, which suggest that darkling beetles from Península Valdés may have narrow tolerance limits to environmental variations than those in the American northern hemisphere.

The present work allows clustering environmental units within the peninsula into five groups, each one with an exclusive darkling beetles species assemblage. The first group of environmental units is represented by established sandy aeolian fields, where perennial grass steppes are the dominant vegetation and photosynthetic activity and soil organic carbon content are the greatest of the peninsula (Rostagno 1981, Bouza et al. 2017a, b). The second group of environmental units is composed by piedmont pediments of the coastal zones and piedmont pediments and “bajadas” of great endorheic basins (salt flats), where more abundant rains and more elevated night temperatures than the remaining environment units of the peninsula are present (Rostagno 1981, Bouza et al. 2017a, b). The third group is constituted by the alluvial plains (terraces) of the peninsula, where daytime temperatures are higher, rains are lower, soils profiles are more complex (due to the presence of hardy horizons that hinder root development), and consequently Chuquiraga avellanedae is the dominant shrub (Rostagno et al. 1981). The fourth and fifth tenebrionid species assemblages are only present in environmental units G4 in coastal piedemont pediments, and G5 in Holocene beach ridges, respectively. It should be noted that in both assemblages the smaller number of tenebrionid species were recorded. It is also remarkable the low incidence of species in the tenebrionid assemblage present on the north-eastern tip of the peninsula (G5) where small mixed patches of Schinus polygamus and Lycium chilense are the only vegetation. Considering the previous discussion on the variables that most affect the distribution of these beetles in the peninsula, it may be caused by the convergence of some particular features that determine a very unfavorable environment for the establishment of the majority of tenebrionid species inhabiting PV. This environment, the youngest within PV, is constituted entirely by beach ridges recently originated during the Holocene from successive coastal marine deposits. Consequently, their poorly developed soils are still large gravelly areas of bare soil without any vegetation nor litter cover, and extremes temperatures (Bouza et al. 2017a, Bertiller et al. 2017).

Considering the ecological relevance that darkling beetles have in arid ecosystems and taking into account the restricted distribution observed for seven species inside PV (C. peninsularis, N. darwini, P. sphaericus, P. (Orthogonoderes) argentina, P. (Hemipraocis) sp.1, S. lacordairei and C. patagonicus), to maximize the chances of an effective preservation of the biodiversity and ecological processes in the region, we strongly recommend to prioritize the conservation efforts on at least some of the environmental units inside each of the groups identified in this study. Since many of the environmental units within PV are similar to other areas in the northeastern Patagonia (Beeskow et al. 1997), the same environmental factors determining the structure and diversity of darkling beetles species in PV may be expected for the rest of the region; however, more studies would be necessary to corroborate this thought.

In agreement with other authors, daytime temperatures and rainfall were among the most important environmental variables that explained the distribution of darkling beetles in PV (de los Santos et al. 2002, Werenkraut & Ruggiero 2012, Liu et al. 2012, Bartholomew & El Moghrabi 2018, Rosas et al. 2019). Since both temperatures and rainfall also determine much of the spatial variation in the structure of beetle communities in the northwest (Werenkraut & Ruggiero 2012) and the southernmost extreme of Patagonia (Rosas et al. 2019), it is expected that these variables will be very important for these insects in all arid Patagonia. According to Parmenter et al. (1989b), we found that temperature was more influential on darkling beetles distribution than precipitation. This is also supported by our field observations on PV about several tenebrionid species performing behaviors related to thermoregulation, e.g. burrowing in the substrate, “shuttling” between open areas exposed to sunshine and shade places, and “stilting”-elevating the body as high as possible above the hot substrate- (Cloudsley-Thompson 2001, Cheli & Martínez 2017). However, in accordance with the results obtained by other authors (e.g. Marino 1986, Crawford 1988, de los Santos et al. 2002), precipitation was also very important to explain variation in darkling beetles distribution patterns inside PV. It is known that conserving corporal water is crucial for the survival of these insects in arid ecosystems (Flores 1998, Cloudsley-Thompson 2001, de los Santos et al. 2002), and the fact that more than 90% of the darkling beetles species recorded in PV have characteristics compatible with high or very high adaptations to retain water (Carrara et al. 2011b) is in agreement with those ideas.

At the same time, soil pH, the number of soil horizons, and soil texture were also very important for darkling beetle community patterns within PV. These findings agree with other authors who pointed out that edaphic variables are important determinants of the structure of darkling beetle communities (e.g. Sheldon & Roger 1984, Ayal & Merkl 1994, Krasnov & Shenbrot 1996, de los Santos et al. 2002, Perner et al. 2005, Werenkraut & Ruggiero 2012, Rosas et al. 2019). Generally, it is assumed that soil variables have mainly indirect effects on insects through their influence on the composition of the plant community (Perner et al. 2005, Werenkraut & Ruggiero 2012); however, the variation partitioning procedures performed in this study showed that there is no shared effect between soil and vegetation variables on the explanation of spatial distribution patterns of darkling beetles in PV. In addition, most of the tenebrionid species present in PV belong to the subfamily Pimeliinae (Flores 1998, Matthews et al. 2010) and to the Tenebrioninae tribe Scotobiini (Kulzer 1955, Matthews et al. 2010) (Table II), which have hypogeal larvae (Flores 1998, Matthews et al. 2010, Silvestro & Michat 2016, Cheli G.H. pers.Obs.). Consequently, we argue that soil variables mostly have a direct influence on the spatial distribution patterns of darkling beetles within PV. This is additionally supported by the fact that in our study the higher incidence and diversity of tenebrionid species were present in the group of environmental units with the most soft-packed soils. These soils have low adhesion between their particles due to their medium size (Rostagno 1981, Bouza et al. 2017b). Consequently, may be easier for darkling beetles to dig in softer soil in their search for shelter, oviposition suitable places, and food supply (Crawford 1988, Krasnov & Shenbrot 1996, de los Santos et al. 2002, Carrara et al. 2011b, Flores et al. 2011). In contrast, soils with high percentage of clay and several horizons (in general related to the presence of argillic horizons) are very hard-packed (see Rostagno 1981 and Bouza et al. 2017b) and thus difficult to dig for several tenebrionid species (Doyen & Tschinkel 1974, Crawford 1988, Krasnov & Shenbrot 1996). This is in agreement with the lower incidence and diversity of tenebrionid species observed in those environmental units in PV with hard-packed soils.

Several authors argue that vegetation, especially shrub cover, is the main determinant of habitat selection for darkling beetles at local scales (see Parmenter et al. 1989a, Mazía et al. 2006, Liu et al. 2012). However, the present study, carried out on a regional scale, did not show that percentage of plant cover, plant species richness nor vegetation structure mainly determine distribution patterns of tenebrionid beetles within PV. Recently Martínez et al. (2018), working locally in steppes of southern PV, found that variables related to vegetation influenced more than soil on the abundance and diversity of epigeal arthropods. At the same time, Mazía et al. (2006), also at local scale, showed a strong influence of shrub patches on darkling beetles habitat use in a steppe of northwest Patagonia. Inconsistencies among our findings with those of the above mentioned studies in Patagonia may be due to the fact that both species distribution patterns and ecological processes are scale-dependent. Thus, the processes to which tenebrionid species distribution patterns respond may be different according to the scale of analyses (Colombini et al. 2005). More research comparing arthropod distribution patterns at different scales in the region will help to clarify these considerations.

The present study identified species with high habitat specificity and fidelity (i.e., characteristic species, like C. peninsularis, L. nordenskioldi, P. hirtus, N. picipes). Although changes in the abundance of characteristic species are useful for monitoring within the habitat to which they are specific, they provide no information on the direction of monitoring ecological changes, because they are restricted to a single ecological state (Mc Geoch et al. 2002). Therefore, in order to maximize the information on habitat quality extracted from bioindicator assemblages and to improve the efficiency of bioindication systems, those authors suggest to study those species having different degrees of preference for different ecological states, i.e. with intermediate habitat specificity (called detector species) jointly with characteristic species. These species can indicate the direction of environmental change by their prevalence as the habitat changed from their less towards their more preferred habitat state. Following this framework, to estimate the direction of possible environmental changes in PV, we strongly recommend to incorporate in futures monitoring of tenebrionid assemblages those species distributed in several environmental units (e.g. P. kuscheli and B. punctulatus) jointly with characteristic species in futures monitoring of tenebrionid assemblages.

Finally, the introduction of domestic livestock in the arid Patagonia produced changes in vegetation and soil that altered fundamental ecosystem processes, increasing desertification and biodiversity loss (Beeskow et al. 1997, Bisigato et al. 2005, 2009, Chartier et al. 2011, Rossi & Ares 2012). Moreover, despite being a nature reserve, PV is not the exception (Cheli et al. 2016b, Baldi et al. 2017, Cheli & Martínez 2017). Consequently, it is to be expected that anthropogenic disturbances, as grazing, would interact with environmental factors determining the composition of tenebrionid assemblages in the area. In this context, taking into account the ecological relevance of tenebrionid beetles in PV (Cheli et al. 2009, Carrara et al. 2011b, Flores et al. 2011, Baldi et al. 2017, Cheli & Martínez 2017), we recommend incorporating these darkling beetle assemblages into future conservation strategies for arid Patagonia, including characteristic and detector species, their relationship with environmental factors and anthropogenic disturbances.

ACKNOWLEDGMENTS

We warmest thanks to Drs. Mario Rostagno, Pablo Bouza, and MSc Ana Maria Beeskow who kindly shared with us many unpublished Península Valdés environmental data, and whose helpful comments have greatly improved early drafts of the present manuscript. We thank to Alejandra Bueno for her language assistance and two anonymous reviewers for their valuable comments that greatly improved this manuscript. We deeply thank IPEEC (CCT CENPAT-CONICET) and Fundación Vida Silvestre Argentina for providing facilities and logistic support; Association of owners of Península Valdés (PROPENVAL); Administrator of the protected area, and Fundación Vida Silvestre Argentina who allowed access to the study area. We are grateful to Secretaría de Turismo y Áreas Naturales Protegidas and Secretaría de Fauna of the Chubut province for providing collecting permissions. Field work was supported by Agencia Nacional de Promoción Científica y Tecnológica (PICT 2012-2660), CONICET (PIP 112-200801-00162; PIP 112-201101-00987; PUE-IPEEC-CONICET Nº 22920160100044) and the International Barcode of Life Project (iBOL). Trip to Budapest Museum of GEF in 2018 to study Kulzer and Kaszab’s type specimens was supported by CONICET, Argentina and by a grant awarded by The Ernst Mayr Grants Committee (Harvard University).

REFERENCES

• ABALLAY FH, FLORES GE, SILVESTRO VA, ZANETTI NI & CENTENO ND. 2016. An illustrated key to and diagnoses of the species of Tenebrionidae (Coleoptera) associated with decaying carcasses in Argentina. Ann Zool 66: 703-726. doi: 10.3161/00034541ANZ2016.66.4.021.
• AYAL Y. 2007. Trophic structure and the role of predation in shaping hot desert communities. J Arid Environ 68: 171-187.
• AYAL Y & MERKL O. 1994. Spatial and temporal distribution of tenebrionid species (Coleoptera) in the Negev Highlands, Israel. J Arid Environ 27: 347-361.
• BALDI R, CHELI G, UDRIZAR SAUTHIER DE, GATTO A, PAZOS GE & AVILA LJ. 2017. Animal Diversity, Distribution and Conservation. In: Bouza P & Bilmes A (Eds), Late Cenozoic of Península Valdés, Patagonia, Argentina. Springer Earth System Sciences. Springer, Cham., p. 263-303.
• BARTHOLOMEW A & EL MOGHRABI J. 2018. Seasonal preference of darkling beetles (Tenebrionidae) for shrub vegetation due to high temperatures, not predation or food availability. J Arid Environ 156: 34-40. DOI: 10.1016/j.jaridenv.2018.04.008.
• BEESKOW AM, DEL VALLE HF & ROSTAGNO CM. 1997. Los sistemas fisiográficos de la región árida y semiárida de la Provincia de Chubut. Delegación Regional SECYT, Puerto Madryn, Chubut, 144 p.
• BERTILLER MB, BEESKOW AM, BLANCO PD, IDASZKIN YL, PAZOS GE & HARDTKE L. 2017. Vegetation of Península Valdés: Priority Sites for Conservation. In: (Eds), Late Cenozoic of Península Valdés, Patagonia, Argentina. Springer Earth System Sciences. Springer, Cham., p. 131-159.
• BISIGATO AJ, BERTILLER MB, ARES JO & PAZOS GE. 2005. Effect of grazing on plant patterns in arid ecosystems of Patagonian Monte. Ecography 28: 561-572.
• BISIGATO AJ, SAÍN CL, CAMPANELLA MV & CHELI GH. 2015. Leaf traits, water stress, and insect herbivory: Is food selection a hierarchical process? Arthropod-Plant Inte 9: 477-485.
• BISIGATO AJ, VILLAGRA PE, ARES JO & ROSSI BE. 2009. Vegetation heterogeneity in Monte Desert ecosystems: A multi-scale approach linking patterns and processes. J Arid Environ 73: 182-191.
• BORCARD D, GILLET F & LEGENDRE P. 2011. Numerical ecology with R. Springer: New York, New York, NY, p. 319.
• BOUZA P, BILMES A, DEL VALLE H & ROSTAGNO CM. 2017a. Late Cenozoic Landforms and Landscape Evolution of Península Valdés. In: Bouza P & Bilmes A (Eds), Late Cenozoic of Península Valdés, Patagonia, Argentina. Springer Earth System Sciences. Springer, Cham., p. 105-129.
• BOUZA P, RÍOS I, ROSTAGNO CM & SAÍN C. 2017b. Soil–Geomorphology Relationships and Pedogenic Processes in Península Valdés. In: (Eds), Late Cenozoic of Península Valdés, Patagonia, Argentina. Springer Earth System Sciences. Springer, Cham., p. 161-190.
• CARRARA R, CHELI GH & FLORES GE. 2011b. Patrones biogeográficos de los tenebriónidos epígeos (Coleoptera: Tenebrionidae) del Área Natural Protegida Península Valdés, Argentina: implicaciones para su conservación. Rev Mex Biodivers 82: 1297-1310.
• CARRARA R, VAZQUEZ D & FLORES GE. 2011a. Habitat specificity can blur the predictions of species−energy theory: a case study of tenebrionid beetles adapted to aridity. J Arid Environ 75(8): 703-710.
• CEPEDA-PIZARRO J, PIZARRO-ARAYA J & VÁSQUEZ H. 2005.Variación en la abundancia de Arthropoda en un transecto latitudinal del desierto costero transicional de Chile, con énfasis en los tenebriónidos epígeos. Rev Chil Hist Nat 78: 651-663.
• CEPEDA-PIZARRO JG. 1989. Actividad temporal de tenebriónidos epígeos (Coleoptera) y su relación con la vegetación arbustiva en un ecosistema árido de Chile. Rev Chil Hist Nat 62: 115-125.
• CHAO A, GOTELLI NJ, HSIEH TC, SANDER EL, MA KH, COLWELL RK & ELLISON AM. 2014. Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecol Monogr 84: 45-67. <https://doi.org/10.1890/13-0133.1>.
• CHARTIER M, ROSTAGNO M & PAZOS GE. 2011. Effects of soil degradation on infiltration rates in grazed semiarid rangelands of northeastern Patagonia, Argentina. J Arid Environ 75: 656-661.
• CHELI GH. 2009. Efectos del disturbio por pastoreo ovino sobre la comunidad de artrópodos epigeos en Península Valdés (Chubut, Argentina). PhD thesis. Universidad Nacional del Comahue, Centro Regional Universitario Bariloche, Río Negro, Argentina.
• CHELI GH, CORLEY J, BRUZZONE O, DEL BRIO M, MARTÍNEZ F, MARTÍNEZ ROMAN N & RÍOS I. 2010. The ground-dwelling arthropod community of Península Valdés (Patagonia, Argentina). J Insect Sci 10:50. <www.insectsicence.org/10.50>
• CHELI GH, CORLEY J, CASTILLO LD & MARTÍNEZ F. 2009. Una aproximación experimental a la preferencia alimentaria de Nyctelia circumundata (Coleoptera: Tenebrionidae) en el Noreste de la Patagonia. Interciencia 34: 771-776.
• CHELI GH, FLORES GE, MARTÍNEZ ROMÁN N, PODESTÁ D, MAZZANTI R & MIYASHIRO L. 2013. Tenebrionid beetle’s dataset (Coleoptera, Tenebrionidae) from Península Valdés (Chubut, Argentina). ZooKeys 364: 93-108.
• CHELI GH & MARTÍNEZ FJ. 2017. Artrópodos terrestres, su rol como indicadores ambientales. In: Udrizar sauthier DE, Pazos GE & Arias AM (Eds), Reserva de Vida Silvestre San Pablo de Valdés: 10 años protegiendo el patrimonio natural y cultural de la Península Valdés. Fundación Vida Silvestre Argentina & CONICET, Buenos Aires, Argentina, p. 98-117.
• CHELI GH, MARTÍNEZ FJ, PAZOS G, UDRIZAR SAUTHIER D, CORONATO F & FLORES GE. 2016a. Epigean Tenebrionids (Coleoptera: Tenebrionidae) from Leones and Tova-Tovita Islands (Chubut, Argentina) and its comparison with the mainland assemblage. Ann Zool 66: 631-642. <https://doi.org/10.3161/00034541ANZ2016.66.4.015>.
• CHELI GH, PAZOS GE, FLORES GE & CORLEY JC. 2016b. Efecto de los gradientes de pastoreo ovino sobre la vegetación y el suelo en Península Valdés, Patagonia Argentina. Ecología Austral 26: 200-211.
• CLARKE KR & WARWICK RM. 2001. Change in marine communities: an approach to statistical analysis and interpretation, 2nd edition. PRIMER-E, .
• CLOUDSLEY-THOMPSON JL. 2001. Thermal and water relations of desert beetles. Naturwissenschaften 88: 447-460.
• COLOMBINI I, FALLACI M & CHELAZZI L. 2005. Micro-scale distribution of some arthropods inhabiting a Mediterranean sandy beach in relation to environmental parameters. Acta Oecol 28: 249-265.
• CORONATO F, PESSACG N & ALVAREZ MP. 2017. The climate of Península Valdés within a regional frame. In: Bouza P & Bilmes A (Eds), Late Cenozoic of Península Valdés, Patagonia, Argentina. Springer Earth System Sciences. Springer, Cham.
• CRAWFORD CS. 1988. Surface-active arthropods in a desert landscape: influences of microclimate vegetation and soil texture on assemblage structure. Pedobiologia 32: 373-385.
• DE LOS SANTOS A, DE NICOLÁS JP & FERRER F. 2002. Habitat selection and assemblage structure of darkling beetles (Col. Tenebrionidae) along environmental gradients on the island of Tenerife (Canary Islands). J Arid Environ 52: 63-85.
• DOYEN JT & TSCHINKEL WR. 1974. Population size, microgeographic distribution and habitat separation in some tenebrionid beetles (Coleoptera). Ann Entomol Soc Amer 67: 617-626.
• FLORES GE. 1998. Tenebrionidae. In: Morrone JJ & Coscarón S (Eds), Biodiversidad de Artrópodos Argentinos volumen 1: 232-240. Ediciones Sur, La Plata, Argentina.
• FLORES GE. 1999. Systematic revision and cladistic analysis of the Neotropical genera Mitragenius Solier, Auladera Solier and Patagonogenius gen. n. (Coleoptera: Tenebrionidae). Entomol Scand 30: 361-396.
• FLORES GE. 2004. Systematic revision and cladistic analysis of the Patagonian genus Platesthes (Coleoptera: Tenebrionidae) Eur J Entomol 101: 591-608.
• FLORES GE, CARRARA R & CHELI GH. 2011. Three new Praociini (Coleoptera: Tenebrionidae) from Peninsula Valdés (Argentina), with zoogeographical and ecological remarks. Zootaxa 2965: 39-50.
• FLORES GE & DEBANDI GO. 2004. Tenebrionidae. In: Catálogo de insectos fitófagos de la Argentina y sus plantas asociadas. Ed. Sociedad Entomológica Argentina, Buenos Aires, p. 197-201.
• FORMOSO AE, TETA P & CHELI GH. 2012. Food habits of the Magellanic horned owl (Bubo virginianus magellanicus) at southernmost Patagonia, Argentina. J Raptor Res 46: 401-406.
• HSIEH TC, MA KH & CHAO A. 2016. iNEXT: an R package for rarefaction and extrapolation of species diversity (Hill numbers). Methods Ecol Evol 7: 1451-1456. doi: 10.1111/2041-210X.12613.
• KAMIŃSKI MJ ET AL. 2021. Reevaluation of Blapimorpha and Opatrinae: addressing a major phylogeny-classification gap in darkling beetles (Coleoptera: Tenebrionidae: Blaptinae). Syst Entomol 46: 140-156. <https://doi.org/10.1111/syen.12453>.
• KRASNOV B & SHENBROT G. 1996. Spatial Structure of Community of Darkling Beetles (Coleoptera: Tenebrionidae) in the Negev Highlands, Israel. Ecography 19: 139-152.
• KULZER H. 1955. Monographie der Scotobiini (Zehnter Beitrag zur Kenntnis der Tenebrioniden). Ent Arb Mus G Frey 6: 383-478.
• KULZER H. 1963. Revision der südamerikanischen Gattung Nyctelia Latr. (Col. Teneb.) (24 Beitrag zur Kenntnis der Tenebrioniden). Entomologische Arbeiten aus dem Museum Georg Frey 14: 1-71.
• LAITY J. 2008. Deserts and desert environments. Wiley-Blackwell Publishing, 342 p. ISBN 978-1-57718-033-3.
• LEGENDRE P & LEGENDRE L. 2012. Numerical ecology, third English edition. Elsevier, Amsterdam, p. 1006.
• LICHTI NI & MURPHY MT. 2010. Selection ratios on community aggregated traits estimate ecological filters imposed on species by sites. Ecology 91: 347-354.
• LIU J, LI F, LIU C & LIU Q. 2012. Influences of shrub vegetation on distribution and diversity of a ground beetle community in a Gobi Desert ecosystem. Biodivers Conserv 21: 2601-2619.
• MARINO PC. 1986. Activity Patterns and Microhabitat Selection in a Desert Tenebrionid Beetle (Coleoptera: Tenebrionidae). Ann Entomol Soc Am 79: 468-471.
• MARTÍNEZ F, CHELI GH & PAZOS G. 2018. Structure of ground-dwelling arthropod assemblages in vegetation units of Área Natural Protegida Península Valdés, Patagonia, Argentina. J Insect Conserv 22: 287-301.
• MATTHEWS EG, LAWRENCE JF, BOUCHARD P, STEINER WE & ŚLIPIŃSKI A. 2010. Tenebrionidae Latreille, 1802. In: Kristensen NP & Beutel RG (Eds), Handbook of Zoology. Arthropoda: Insecta. Volume 2: Coleoptera, Beetles. Morphology and Systematics (Elateroidea, Bosthichiformia, Cucujiformia partim). DeGruyter, Berlin.
• MAZÍA NC, CHANETON E & KITZBERGER T. 2006. Small-scale habitat use and assemblage structure of ground-dwelling beetles in a Patagonian shrub steppe. J Arid Environ 67: 177-194.
• MC GEOCH MA, VAN RENSBURG BJ & BOTES A. 2002. The verification and application of bioindicators: a case study of dung beetles in a savanna ecosystem. J Appl Ecol 39: 661-672.
• NORTHRUP JM, HOOTEN MB, ANDERSON CHR, ANDAND JR & WITTEMYER G. 2013. Practical guidance on characterizing availability in resource selection functions under a use–availability design. Ecology 94: 1456-1463.
• OKSANEN J. 2019. Multivariate analysis of ecological communities in R: Vegan tutorial. <https://cran.r-project.org/web/packages/vegan/vegan.pdf>. Accessed 17 Nov 2019.
  » https://cran.r-project.org/web/packages/vegan/vegan.pdf
• PARMENTER RR, PARMENTER CHA & CHENEY CD. 1989a. Factors Influencing Microhabitat Partitioning among Coexisting Species of Arid-Land Darkling Beetles (Tenebrionidae): Behavioral Responses to Vegetation Architecture. Southwest Nat 34: 319-329.
• PARMENTER RR, PARMENTER CHA & CHENEY CD. 1989b. Factors influencing partitioning in arid-land darkling beetles (Tenebrionidae): temperature and water conservation. J Arid Environ 17: 57-67.
• PERNER J, WYTRYKUSH C, KAHMEN A, BUCHMANN N, EGERER I, CREUTZBURG S, ODAT N, AUDORFF V & WEISSER WW. 2005. Effects of plant diversity, plant productivity and habitat parameters on arthropod abundance in montane European grasslands. Ecography 28: 429-442.
• PIZARRO-ARAYA J. 2010. Hábitos alimenticios del género Gyriosomus Guérin-Méneville, 1834 (Coleoptera: Tenebrionidae): qué comen las vaquitas del desierto costero? IDESIA 28: 115-119.
• QGIS DEVELOPMENT TEAM. 2018. QGIS Geographic Information System. Open Source Geospatial Foundation Project. <http://qgis.osgeo.org>.
  » http://qgis.osgeo.org
• R CORE TEAM. 2018. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. <https://www.R-project.org>.
  » https://www.R-project.org
• ROSAS YM, PERI PL, CARRARA R, FLORES GE, PEDRANA J & MARTÍNEZ PASTUR G. 2019. Potential biodiversity map of darkling beetles (Tenebrionidae): Environmental characterization, land uses and analyses of protection areas in Southern Patagonia J Insect Conserv. 23: 885-897. <https://doi.org/10.1007/s10841-019-00170-w>.
• ROSSI MJ & ARES JO. 2012. Close Range Stereophotogrammetry and Video Imagery Analyses in Soil Ecohydrology Modelling. The Photogrammetric Record 27: 111-126.
• ROSTAGNO CM. 1981. Reconocimiento de los suelos de la Península Valdés. Publicación Nº 44. Centro Nacional Patagónico, Chubut Argentina (inédito), 24 p. Puerto Madryn.
• RUOKOLAINEN L & BLANCHET G. 2014. Introduction to Ecological Multivariate Analysis. University of Helsinki, 79 p. <https://blogs.helsinki.fi/laruokol/files/2016/10/IEMA.pdf>.
  » https://blogs.helsinki.fi/laruokol/files/2016/10/IEMA.pdf
• SCHOWALTER TD. 2016. Insect Ecology. An ecosystem approach. Elsevier, 762 p.
• SHELDON JK & ROGERS LE. 1984. Seasonal and habitat distribution of tenebrionid beetles in shrubsteppe communities of the Hanford site in eastern Washington. Environ Entomol 13: 214-220.
• SILVESTRO VA & FLORES GE. 2016. Notes on types of some species described by Billberg and Germar belonging to the South American genera Scotobius Germar, 1823 and Nyctelia Berthold, 1827 (Coleoptera: Tenebrionidae), with new synonymies and taxonomic implications. Ann Zool 66: 653-663.
• SILVESTRO VA & MICHAT MC. 2016. Description of the egg and first-instar larva of Scotobius pilularius Germar (Coleoptera: Tenebrionidae: Scotobiini), and analysis of larval primary chaetotaxy. Ann Zool 66: 681-691.
• WARD D. 2009. The Biology of Deserts. Oxford University Press, Oxford, UK, 352 p.
• WERENKRAUT V & RUGGIERO A. 2012. Altitudinal variation in the taxonomic composition of ground-dwelling beetle assemblages in NW Patagonia, Argentina: environmental correlates at regional and local scales. Insect Conserv Diver 6: 82-92. doi: 10.1111/j.1752-4598.2012.00190.x.
• WHITFORD WG. 2002. Ecology of desert systems. Academic Press, an Elsevier Science Imprint, San Diego, California, 343 p. ISBN 0-12-747261-4.

Publication Dates

• Publication in this collection
  10 May 2021
• Date of issue
  2021

History

• Received
  10 Aug 2020
• Accepted
  14 Feb 2021