Biogeographical Patterns of Earwigs in Italy

Simple Summary Italy plays a central role in the research on Europe’s biogeography because of its position in the middle of the Mediterranean, which is a global hotspot of diversity. This study investigated how present patterns of earwig species richness and composition in Italy are affected by climatic, geographical, and historical factors. Earwig richness does not decrease from the base to the tip of the Italian peninsula, which contrasts with the so-called ‘peninsula effect’. However, richness does not increase southward either, suggesting that southern regions did not play a crucial role as a refuge during Pleistocene glaciations. Inter-regional similarities in species composition between regions is more influenced by their geographical proximity than climatic similarity, although richness is positively correlated with precipitation, in accordance with earwig preferences for humid conditions. The similarity in species composition with central European fauna decreases southward, indicating possible exchanges between central Europe and Italy. The majority of the earwigs of the Italian fauna are either widespread across Europe and the Palearctic, or confined to the main Italian mountain ranges, that is, the Alps and the Apennines. The isolation of ancient earwig populations on mountains resulted in the development of a high proportion of endemics, making the Italian earwig fauna one of the richest in Europe. Abstract Placed in the center of the Mediterranean biodiversity hotspot, Italy plays a central role for the study of Europe’s biogeography. In this paper, the influence of climatic, spatial, and historical factors on current patterns of variation in earwig species richness and composition is investigated. The Italian earwig fauna is mainly composed of species which are either widely distributed in Europe and the Palearctic region or that are endemic to the Alps and the Apennines. Variation in species richness does not follow any obvious geographical patterns, but a positive influence of precipitation on richness is consistent with earwig preferences for humid climates. European mainland territories did not contribute substantially to the current biodiversity of Italian earwigs, which explains the lack of a distinct peninsula effect, although a southward decrease in similarity with the central European fauna was observed. However, southern areas did not exert a pivotal role during Pleistocene glaciations in determining current patterns of species richness. Variation in species composition among Italian regions can be mostly explained by geographical proximity, while climatic differences and historical (paleogeographical and paleoecological) events seem to have played a minor role. However, the isolation of ancient earwig stocks on Italian mountains led to the origin of a relatively large number of endemics, which makes the Italian earwig fauna one of the richest in Europe.


Introduction
Earwigs (Dermaptera) are a small order of hemimetabolous insects comprising about 2000 known species worldwide [1]. Most earwig diversity is concentrated in the tropical regions, whereas relatively few species inhabit the temperate zone of the Northern Hemisphere [2]. This has been interpreted as a result of the Gondwanan origin of earwigs, of ancient processes of allopatric speciation due to the effects of paleogeographical and paleoecological changes on the distribution of ancient populations. Prediction 5. Biogeographical similarities between Italian regions and adjacent countries reflect post-glacial colonization trajectories. Previous research revealed multiple colonization trajectories, both in tenebrionids [22] and odonates [24]. For tenebrionids, during the Pleistocene glacials, southern Italian regions (which have high levels of endemicity) acted as an important refugial center for thermophilous species of North African or east Mediterranean origin, while cold-adapted species coming from northern areas (Central Europe) dispersed into central and southern Italy through the main mountain ranges [22]. For the odonates, the Italian area may have played an important role as a refugial center, but with deglaciation, these insects moved towards central European areas, from which they re-colonized the Italian peninsula, leading to a pattern of southward decreasing biogeographical similarities with European faunas [24].
Thus, we can expect: (1) high similarities of Northwestern Italian regions with the French fauna through the Provencal area; (2) southward decreasing similarity with the Central European fauna; (3) high similarities of Northeastern Italian regions with Balkan territories through the Karst Plateau; and, finally, (4) high similarities of Southern regions with the Northern African fauna because of the possible persistence of a pre-Pleistocene (Tertiary) fauna and post-Pleistocene immigration of thermophilous species.

Statistical Analyses
Before analysis, spatial autocorrelation between richness values and region centroids was tested by using the Moran I index, which showed no significant spatial autocorrelation (I = −0.077, p = 0.839).
Predictions 1a and 1b were tested by correlating (Pearson product-moment correlation coefficient) species richness with the centroid latitude of mainland regions. Prediction 2 was tested by evaluating the importance of the following climatic variables in determining species richness: average total annual precipitation (Pmean), average annual temperature (Tmean), mean minimum temperature (Tmin), mean maximum temperature (Tmax), and yearly temperature difference (∆T = Tmax − Tmin) [22,24,39,40]. Geographical parameters (i.e., area, latitude, and longitude of regions) were also considered. Geographical and climatic data were taken from Fattorini [22]. The importance of these variables in explaining species richness was tested using a multimodel inference approach by running models representing every possible combination of explanatory variables and then ranking them according to the values of the corrected Akaike Information Criterion (AICc). Models with a difference in the AICc values lower than two were averaged using both full and conditional averages. In the full average, regression coefficients for variables that are not included in a given model are set to zero, whereas conditional average only averages over the models where the parameter appears.
Predictions 3, 4a, and 4b were tested by correlating inter-regional biogeographical distances with geographical and climatic distances. For this purpose, biogeographical distances between regions were expressed using the overall ß-diversity (ßsor; that is 1-Sørensen index of similarity), the pure turnover component (ßsim; that is 1-Simpson index of similarity, which expresses compositional differences independently from the influence of nestedness), and the nestedness component (ßnest; that is ßsor-ßsim, which quantifies the part of the compositional change caused by ordered species loss) [59,60]. Geographical distances were calculated as distances between centroids. Climatic distances were calculated as Euclidean distances for the climatic variables reported above after standardization. Matrices were correlated using Mantel tests (to evaluate the effect of either geographical position or climatic conditions on biogeographical distances) and partial Mantel tests (to evaluate the effect of geography controlling for the effect of climate, and to evaluate the effect of climate controlling for the effect of geography) [22]. Additionally, biogeographical relationships between regions expressed by Sørensen and Simpson indexes were represented by using Non-Metric Multidimensional Scaling (NMDS), an ordination technique that allows the representation of dissimilarities among areas in a two-dimensional space, and which is therefore particularly effective at disclosing multiple relationships in biogeographical analyses [24,61,62]. The aim of NMDS is to find a spatial configuration of points in which the distances between pairs of points in the configuration matches as well as possible the original dissimilarities between the pairs. To this end, NMDS starts with random initial conditions and iteratively seeks new solutions which are evaluated using a measure of goodness of fit called stress (lower values indicating better fit). Procrustes distances were used to compare solutions until a minimum stress value was obtained. In the two-dimensional representation, the axis with the highest variance was standardized between 0 and 1, and the other axis was rescaled according to the first one. Next, the colors blue, green, yellow, and red were assigned to the four corners, and each area received an RGB (red, green, blue) color according to its position in the two-dimensional space.
Prediction 5 was tested by considering the distribution of each species of the Italian fauna in the following adjacent areas: Western Europe (fauna of France), Central Europe (faunas of Austria, Switzerland, and Germany), Eastern Europe (faunas of Slovenia, Albania, Bosnia, Croatia, Herzegovina, Montenegro, and mainland Greece), and Northern Africa (Tunisia and close areas). Next, the similarity (Sørensen and Simpson indexes) between these areas and the Italian regions was calculated and mapped. The Pearson correlation coefficient was used to test for correlation between these biogeographical similarities and latitude (use of Moran I tests showed non-significant spatial autocorrelations).
Finally, to describe the biogeographical composition of the Italian earwig fauna, species were assigned to chorotypes (i.e., groups of species with similar distributions [63,64]) following the nomenclature proposed by Vigna Taglianti et al. [65].

Results
Earwig species richness did not show a latitudinal pattern (correlation between number of species and latitude: r = 0.026, p = 0.929, Figure 1b). However, very low values were observed in Apulia (region 14), Sicily (region 15), and Sardinia (region 16), which are particularly arid regions, and average total annual precipitation (Pmean) was the most important factor influencing earwig species richness in Italy (Table 1). Thus, in contrast with Prediction 1b, species richness did not decrease southward in response to decreasing rainfall, although precipitation was an important predictor of richness, as expected according to Prediction 2. A possible positive influence of temperature appeared when the conditional average was used.
The NMDS with Sørensen index (Figure 2a) showed a strong separation of region 2 (Trentino-Alto Adige) from all other regions. Alpine regions 1, 3, and 4 formed a distinctly recognizable group. The regions south of the Po River formed a complex group in which the Tyrrhenian regions (regions 8, 9, 10, 11, and 12) appeared to be closely associated. Apulia (region 14) occupied a very distinct position, while the three islands appeared close to the adjacent Italian mainland areas: Sicily (region 15) was close to Calabria (region 13) and Campania (region 12); Sardinia (region 16) was close to regions facing the northwestern sector of the Tyrrhenian Sea (regions 5 and 6); and Corsica (region 17) was close to Liguria (region 6). The results of NMDS with Simpson index (Figure 2b) were very similar to those obtained with the Sørensen index in showing: (1) a strong separation of region 2 (Trentino-Alto Adige) from all other regions, (2) the presence of a group of Alpine regions (regions 1, 3, and 4), and (3) complex relationships between the regions south of the Po river and the three islands. Sicily (region 15) resulted to be very close to Calabria (region 13), whereas Sardinia (region 16) was very close to Corsica (region 17) and Liguria (region 6). These results indicate a major role exerted by current geographical settings (the position of the main mountain ranges, i.e., the Alps at north and the Apennines at south, and geographical proximity) more than by paleogeographical and paleoecological history, therefore being contrary to Prediction 4b.
Partial Mantel tests (Table 2) showed that correlations between biogeographical dissimilarities (calculated as either Sørensen or Simpson index) and geographical position were significant even after controlling for climate, whereas correlations between biogeographical dissimilarities and climate were not significant after controlling for geographical position. This indicates an importance of the geographical position independently from the effect of climate (thus supporting Prediction 4a), whereas (in contrast with Prediction 3) there was no relevant influence of climate on biogeographical dissimilarities independently from the geographical position. The similarity between Italian regions and adjacent areas expressed by the Sørensen and Simpson indices pointed to some recognizable geographical patterns (Figures 3 and 4).
The  (Figure 4b), whereas those for the Eastern fauna were observed in region 1 (which includes the Eastern Alps) and region 14 (Apulia, the easternmost Italian region) (Figure 4c). Values for the North African fauna were rather uniformly high because most of the few species included in this category are very widespread (Figure 4d). For the Central European fauna, Sørensen similarities decreased southward, whereas no significant correlation was found between latitude and similarity for the Western European and the Eastern European fauna (Table 3). A slightly negative correlation was found between the latitude and Sørensen similarities for the African fauna (Table 3). No significant correlation was found for the Simpson index. These results are in partial agreement with Prediction 5.

Discussion
Italian earwigs do not show any distinct latitudinal gradient in species richness, which contrasts with previous findings showing either an increasing diversity southward (as observed in dung beetles and tenebrionids [22,74]), nor an increasing diversity northward (as found in birds, small mammals, ondonates, carabid beetles, hydradephagan beetles, and ants [24,42,75]). Thus, earwigs follow neither Prediction 1a ('latitudinal gradient') nor Prediction 1b ('peninsular effect'). This unexpected result may be explained by the possible contrasting effects of the spatial distribution of temperatures (which increase southward) and humidity (which increases northward) [76]. In general, earwigs tend to prefer warm and humid climates. In Italy, however, the warmest areas (located in the southern parts of the country) are also the driest, whereas the most humid regions (located in the north) are also the coldest, as shown by the distribution of Köppen-Geiger climate types [76]. Although inter-regional biogeographical relationships were correlated with climatic dissimilarities, this correlation disappeared when controlling for geographical position, which contrasts with the hypothesis of a primary role exerted by climatic conditions in shaping variations in species composition, as postulated by Prediction 3, but supports a strong influence of geographical position independently from climate, as postulated by Prediction 4a. This indicates that earwig species distributions in Italy were more profoundly shaped by geographical barriers/connections than by climate. The results of the ordination analysis (NMDS) indicate: (1) a major division between the Alpine and the Apennine faunas; (2) a distinct separation of region 2 (Trentino-Alto Adige, which is the coldest region) from all other regions; (3) a high similarity between Sicily (region 15) and the adjacent mainland (Calabria, region 13); (4) a high similarity between Corsica (region 17) and the adjacent mainland (Liguria, region 6); (5) a strong similarity between Sardinia (region 16) and both Corsica and Liguria. These patterns point to the role of the main geographical characteristics of the study area in shaping biogeographical similarities, whereas the possible importance of historical factors (Prediction 4b) seems to be secondary. Even the similarity between Corsica and Sardinia (already observed in other insects, such as butterflies [48], tenebrionids [22], carabids [46], and chrysomelids [46], but not in odonates [24]), and their similarity with the adjacent mainland, seem to be more easily interpretable for earwigs as a result of geographical proximity than as evidence of their shared paleogeographical history as a microplate (Corsardinia: Corsica + Sardinia) detached from the Provencal area [22]. However, it is important to note that about 40% of the native earwig fauna of Italy (including Sicily, Corsica, and Sardinia) is endemic to this area (and this proportion would be even higher considering that the distribution of Chelidurella pseudovignai Kočárek & Kirstová, 2020 and C. mutica (Krauss, 1886) out of Italy is restricted to areas close to the borders of the study area), and relatively high proportions of the endemic species are found in most of the regions. The earwigs endemic to Italy can be classified into two main groups whose distributions correspond to the two main mountain ranges: the Alpine endemics (Chelidura aptera (Megerle in Charpentier, 1825), Chelidurella poggii Capra, 1982, and C. vignai Galvagni, 1995) and the Apennine endemics (Chelidurella caprai Vigna Taglianti, 1993, Pseudochelidura galvagnii Vigna Taglianti, 1999, P. orsinii (Gené, 1833), Forficula apennina Costa, 1881, and Forficula silana Costa, 1881) [4,8], which is consistent with the recognized importance of mountain ranges for cladogenetic processes in these animals in Europe [77].
In agreement with Prediction 5 and findings for other groups [22,24], Italian earwigs showed a southward decline in their similarity with the Central European fauna and an opposite pattern for their similarity with the North African fauna. However, in contrast with Prediction 5, there were no clear relationships with the Western and Eastern European faunas. Additionally, no correlation was found when the pure turnover component (Simpson index) was used. Looking at the distribution of the species in the Italian territory, it appears that the Italian fauna is composed of a mixture of species with very restricted distributions on one hand (with many endemics), and widely distributed species on the other hand (with species widely distributed in Europe, the Palearctic region, or even broadly). This suggests that the Italian territory was invaded by stocks of species with high dispersal capabilities without distinct colonization routes. However, some ancient events of colonization resulted in the isolation on mountain areas of the ancestors that evolved into the current endemic species.

Conclusions
The Italian earwig fauna is mainly composed of either species which are widely distributed in Europe and the Palearctic region, or species endemic to the two main mountain ranges, that is, the Alps and the Apennines. Variation in species richness does not follow any obvious geographical pattern, but a positive influence of precipitation on richness is consistent with the ecological preferences of these insects for humid climates. Contrary to what could be expected on the basis of the 'peninsula effect', earwig richness did not decrease southward along the Italian peninsula, thus indicating that European mainland territories did not contribute much to the current biodiversity of Italian earwigs, at least not through distinctly recognizable pathways (with exception of a certain southward decrease of similarity in species composition with the central European fauna). However, an opposite pattern of increasing diversity southward was not observed either, which indicates that the possible refugial role of southern areas during Pleistocene glaciations did not exert a pivotal role in determining current patterns of species richness in Italy. The variation in species composition among Italian regions can mostly be explained by geographical proximity, while climatic differences and historical (paleogeographical and paleoecological) events seem to be less important. However, the isolation of ancient earwig stocks on Italian mountain regions led to the origin of a relatively large number of endemics, which makes the Italian earwig faunas one of the richest in Europe.
Funding: This research received no external funding.
Data Availability Statement: Species distribution data are given in Table S1. Geographical and climatic data are available from Fattorini [25].