Environmental Factors Drive Chalcid Body Size Increases with Altitudinal Gradients for Two Hyper-Diverse Taxa

Simple Summary The size of an organism is closely correlated with its physiological and ecological characteristics and strongly influences its fitness. Bergmann’s rule, originally widely applied in homeotherms, states that individuals living in a colder environment are larger than those living in a warmer region. Similar geographical patterns are also found in various groups of poikilotherms, but are still controversial, especially in Hymenoptera. In this investigation, we found a significant upward trend in body size with increasing elevation in two tiny groups of Chalcids (Pteromalidae and Eulophidae). The temperature and precipitation play a crucial role in their size variation. This result casts light on the environmental adaptation of parasitoids. Abstract Body size is the most essential feature that significantly correlates with insects’ longevity, fecundity, metabolic rate, and sex ratio. Numerous biogeographical rules have been proposed to illustrate the correlation between the body sizes of different taxa and corresponding geographical or environmental factors. Whether the minute and multifarious chalcids exhibit a similar geographical pattern is still little known. In this research, we analyzed morphological data from 2953 specimens worldwide, including the two most abundant and diverse taxa (Pteromalidae and Eulophidae), which are both composed of field-collected and BOLD system specimens. We examined forewing length as a surrogate of body size and analyzed the average size separately for males and females using two methods (species and assemblage-based method). To verify Bergmann’s rule, we included temperature, precipitation, wind speed and solar radiation as explanatory variables in a generalized linear model to analyze the causes of the size variation. We found that there was an increasing trend in the body size of Pteromalidae and Eulophidae with altitude. The optimal Akaike information criterion (AIC) models showed that larger sizes are significantly negatively correlated with temperature and positively correlated with precipitation, and the possible reasons for this variation are discussed and analyzed.


Introduction
Body size is one of the most remarkable characteristics of all organisms, since it is significantly correlated with the biological, physiological, and ecological traits involving longevity, fecundity, metabolic and reproduction rate, and population size [1][2][3]. Bergmann's rule [4], a widely recognized ecogeographic rule, states that animals living in colder climates are larger than those living in a warmer regions. The geographical pattern has been broadly verified in mammals [5] and birds [6], which have broad ecological niches ranging from tropical to arctic, and their body size is positively correlated with latitude.
By summarizing the relevant literature (Table 1), we found that most of the investigations into Hymenoptera were focused on ants, bees, and other larger wasps. The geographical trends vary among species and different morphological parts, and some ant species even showed a nonlinear trend with elevation [19]. Relevant research on miniature taxa only involved Encyrtidae [20] and Mymaridae [21], and only conducted a simple comparison of body-size differences, so the knowledge of small-sized taxa is still limited. The forewing length and hind tibia length are two generally used surrogates for body length in laboratory experiments, but this method lacks validation in multiple wild species. Meanwhile, they are closely correlated with longevity, fecundity, reproductive duration, and oviposition index [22,23], which provide insights into the relationship between size and environmental factors.
The various geographical trends are generally caused by the combined effect of multiple environmental factors [24]. Based on the results of several models' analyses, the explanatory mechanisms of these trends are usually correlated with temperature and precipitation in Diptera [25,26], Lepidoptera [13], Odonata [27,28], and Orthoptera [29]. Numerous hypotheses related to body-size control in parasitoids have been widely verified in laboratory experiments [30][31][32]. In consequence, it is critical to study these relationships under realistic natural conditions. As we all know, temperature significantly decreases with increasing elevation and is linearly correlated [33]; the most frequently verified hypotheses is the temperature-size rule [34,35], which states that a higher temperature may accelerate the development rate and result in adults with a smaller body size than those living in lower temperatures; therefore, insects living in lower temperatures tend to be relatively larger. Resources availability is also crucial to insect survival and development; insects tend to reduce their reproductive investment and increase energy storage to maintain survival under the stress of starvation [36]; thus, larger-sized insects have the potential to store more resources to resist that stress [37] and can also remain warm in cooler environments due to their lower surface-area-to-volume ratio. The wing-size polymorphism and variations are also closely related to altitude among flying insects. Wing loading and the wing aspect ratio are two common parameters for flight capability. The lower flight activity and larger wing size in a highland environment compensate for the thinner air encountered during flight [24]. However, it is still unknown whether these rules also apply to minute parasitic insects.
To address these gaps, we primarily focused on the geographical patterns of the forewing length of Pteromalidae and Eulophidae, based on a large-scale dataset in this analysis, which consists of field-collected samples from Northwestern China and data from the BOLD system. The main objectives of this research were to analyze the relationship between morphological characteristics and elevation; verify whether the body size increases with increasing elevation, as Bergmann's rule states; investigate the dominant environmental factors contributing to this trend, such as temperature and precipitation; and offer alternative explanations for the relationship, including the temperature-size hypothesis and the resource's availability hypothesis. We hope to provide insights that could improve the understanding of phenotypic plasticity in the future.

Study Area and Specimen Measurement
This study was based on a compiled dataset of 2953 specimens, which included 1406 Pteromalids (52 genera and 81 species) and 1547 Eulophids (54 genera and 92 species), including endemically common and cosmopolitan species. More detailed information and the composition of the dataset is shown in Table S1.
The field-collected specimens were primarily from Northwestern China (mainly covering the endemic habitats of Xinjiang, Gansu, Ningxia, and western Inner Mongolia). The area provides full-scale altitudinal gradient transect ranges from 220 m in the basin to 4336.9 m in the Altun mountain national nature reserve (Table S2). The specimens were collected using net-sweeping and yellow pan traps (malaise traps were used at two alpine stations) during July and August from 2017 to 2020, then labeled and identified definitely to the level of genus or species following a combination of classical taxonomic guides [52][53][54][55][56], and preserved in the Insect Collection of the College of Life Science and Technology, Xinjiang University, Urumqi, Xinjiang, China (ICXU). To increase the sample quantity and geographic coverage, we selected specimens from the BOLD system (http://www.boldsystems.org/. accessed on 1 February 2021), with a complete morphological graph (with scale bar), taxonomy, and collection site, to add to the dataset. The final dataset covers 217 collection sites worldwide (sites less than 10 km apart were merged into one site). We also extracted the GPS coordinates of each point for subsequent analysis. The measured specimens finally used for analysis were distributed across all the continents except the poles (Figure 1a).
We photographed each collected specimen using a Nikon SMZ25 stereomicroscope (Nikon, Japan) in conjunction with an NIS-Elements D 4.30.00 image-collecting system using 10× magnification to capture the entire body in a single image. To avoid deviation, for accurate measurements, and to increase the reliability of the data, any severely shriveled or damaged specimens were ignored and the average value of multiple measurements for each characteristic was calculated. We measured the body length (BL) from the top of the head to the end of the abdomen, excluding the ovipositor protrusion, the mesosoma length (MEL) from the anterior edge of the pronotum to the posterior margin of the propodeum, the forewing length (FL) from the tegula proximally to the outermost point of its distal end, and the hind tibia length (HTL) (Figure 1b). The forewing was spread on the microscope slides for accurate measurement. The wing area (WA) was measured by drawing the contour of the forewing, starting and ending at the tegula ( Figure 1b). After summarizing the relevant literature, we found that the wing loading (WL) and wing aspect ratio (WAR) also have a certain relationship with altitude [13,57]. Since it is difficult to obtain accurate body mass data on parasitoids, the WL was replaced by MEL 3 /WA 2 [58], and the WAR was evaluated by FL 2 /WA [59]. All the measurements were conducted using ImageJ/Fiji 1. 46. stations) during July and August from 2017 to 2020, then labeled and identified definitely to the level of genus or species following a combination of classical taxonomic guides [52][53][54][55][56], and preserved in the Insect Collection of the College of Life Science and Technology, Xinjiang University, Urumqi, Xinjiang, China (ICXU). To increase the sample quantity and geographic coverage, we selected specimens from the BOLD system (http://www.boldsystems.org/. accessed on 1 February 2021), with a complete morphological graph (with scale bar), taxonomy, and collection site, to add to the dataset. The final dataset covers 217 collection sites worldwide (sites less than 10 km apart were merged into one site). We also extracted the GPS coordinates of each point for subsequent analysis. The measured specimens finally used for analysis were distributed across all the continents except the poles ( Figure 1a). We photographed each collected specimen using a Nikon SMZ25 stereomicroscope (Nikon, Japan) in conjunction with an NIS-Elements D 4.30.00 image-collecting system using 10× magnification to capture the entire body in a single image. To avoid deviation, for accurate measurements, and to increase the reliability of the data, any severely shriveled or damaged specimens were ignored and the average value of multiple measurements for each characteristic was calculated. We measured the body length (BL) from the top of the head to the end of the abdomen, excluding the ovipositor protrusion, the mesosoma length (MEL) from the anterior edge of the pronotum to the posterior margin of the propodeum, the forewing length (FL) from the tegula proximally to the outermost point of its distal end, and the hind tibia length (HTL) (Figure 1b). The forewing was spread on the microscope slides for accurate measurement. The wing area (WA) was measured by

Environmental Data
Based on the previous predictions and hypotheses of the size-related Bergmann's rule, wind speed has a certain effect on parasitoids' flight ability [60], and the lower water-loss hypothesis for larger insects in drought climates [61,62], we only considered annual mean precipitation, temperature, wind speed and solar radiation to analyze their effects on the morphological traits. The environmental variables used in our study were obtained from WorldClim-Global Climate and Weather data (version 2.1 available at https: //www.worldclim.org/data/worldclim21.html. accessed on 10 April 2022.), developed by Fick and Hijmans [63]. The four environmental factors for each site were extracted using the 'surface tool' in ArcGIS Spatial Analyst [64], based on the geographical coordinates of the specimens used in our study; the latter two variables were obtained by calculating the average value of 12 months of data; all the digital information was acquired at the spatial level resolution of 30 arcseconds.

Statistical Analysis
We identified the best morphological characteristics that can be used as a proxy for body length by using an OLS regression of seven morphological parameters from males and females of Pteromalidae and Eulophidae. After establishing that the FL is more closely related to the BL than other characteristics (Pteromalidae: Cor = 0.931; Eulophidae: Cor = 0.953), we analyzed the geographical patterns of FL with altitude for both sexes, using two methods at the family level: (1). In the species-based method, the arithmetic means of FL were calculated separately for each species per location, and species with fewer than three records were excluded; a total of 2134 (1233 Eulophidae and 901 Pteromalidae) specimens were used in this analysis. (2). For the assemblage-based method, the arithmetic means of FL were calculated for all specimens from the same location, ignoring the species; some sites with fewer than three records were excluded; 2706 (1420 Eulophidae and 1286 Pteromalidae) specimens were analyzed by this method. This method aims to reflect the complete situation of populations at a specific location and is more ecologically meaningful. The generalized linear model (GLM) is more inclusive, allowing for dependent variables to be non-normally distributed. We used the FL of the different sexes as dependent variables in the GLM for the species-based and assemblage-based methods, and the corresponding environmental factors of each site were used as the explanatory variables. A gamma distribution of errors and a log link function were used in the GLM, while the stepAIC () function was applied to avoid collinearity and simplify the models according to the Akaike information criterion (AIC) values.
We analyzed the relationship between the FL and altitude for six dominant subfamilies with the largest sample sizes. In the genus-level analysis, we accessed the relationship between the altitude (dependent variable) and the FL of 16 widely spread genera (sample size: at least 30 specimens) from Pteromalidae and Eulophidae (independent variables), taking into account the interaction between sex and altitude. We indicated the correlation structure using the corAR () function to account for the spatial autocorrelation in the data and also used the AIC value to determine the best-fit model.
The models and statistical analyzes were performed in R (version 4.0.5, R Core Team, 2021) and Excel; the data visualization was conducted by using ggplot2, ggpubr, and ggpmisc; and the model results were checked by using the ANOVA function in the "car" package [65]. The R model code refers to the code in the related articles [25,26].

Results
The regression results of the FL and HTL on the BL for females and males of 51 Pteromalidae and 48 Eulophidae species confirmed the validity of the FL as the best surrogate for body size compared to the HTL (Figure 2, Figure S1). Meanwhile, the WL was considerably negatively correlated with the FL (Pteromalidae: cor = −0.639; Eulophidae: cor = −0.604), while there was no significant correlation between the WAR and the FL (Pteromalidae: cor = 0.024; Eulophidae: cor = −0.059). The FL increased significantly with altitude and decreased with temperature in both taxa for both sexes, and in both species and assemblage-based methods ( Figure 3, Table 1).
The results of the species-based method revealed a more pronounced difference between sexes for both taxa; the mean FL of females was obviously larger than that of males, indicating that larger females dominated at all altitudinal ranges. Pteromalids are generally larger than Eulophids ( Figure S2, all sexes combined; mean FL of Pteromalids = 1.51 mm, mean FL of Eulophids = 1.05 mm, t. test: p < 0.001); the FL size frequency of both taxa demonstrated significant right skewness ( Figure S3), and the frequency peaks in both sexes transferred from left to right, indicating that the proportion of individuals with larger FL tends to be greater at a higher altitude, and this was remarkably evident in the males of both taxa. The slope values showed that the trend of FL variation with altitude was more pronounced in the assemblage-based method (Figure 3).
A stepwise best subset regression, including four environmental factors and sex, was conducted. The results of the GLM model on the FL from the two methods showed that it was negatively related to sex in both taxa (males have smaller FL, p < 0.01), and positively correlated with the annual mean precipitation (p < 0.01) and annual mean solar radiation, while negatively related to the annual mean temperature (p < 0.001) ( Table 2). The annual mean wind speed was excluded and not included in any of the models with the lowest AIC. The temperature consistently provided the greatest independent explanation for the FL variation in both methods. omalidae and 48 Eulophidae species confirmed the validity of the FL as the best surrogate for body size compared to the HTL (Figure 2, Figure S1). Meanwhile, the WL was considerably negatively correlated with the FL (Pteromalidae: cor = −0.639; Eulophidae: cor = −0.604), while there was no significant correlation between the WAR and the FL (Pteromalidae: cor = 0.024; Eulophidae: cor = −0.059). The FL increased significantly with altitude and decreased with temperature in both taxa for both sexes, and in both species and assemblage-based methods ( Figure 3, Table 1).  The results of the species-based method revealed a more pronounced difference between sexes for both taxa; the mean FL of females was obviously larger than that of males, indicating that larger females dominated at all altitudinal ranges. Pteromalids are generally larger than Eulophids ( Figure S2, all sexes combined; mean FL of Pteromalids = 1.51 mm, mean FL of Eulophids = 1.05 mm, t. test: p < 0.001); the FL size frequency of both taxa demonstrated significant right skewness ( Figure S3), and the frequency peaks in both sexes transferred from left to right, indicating that the proportion of individuals with  Table 2. The generalized linear model (GLM) regression of relationships between FL of Pteromalidae and Eulophidae and environmental variables. AMTEM: annual mean temperature ( • C), AMPRE: annual mean precipitation (mm), AMSR: annual mean solar radiation (kJ/m 2 ), AMWS: annual mean wind speed (m/s). (** </ =0.01; *** </ =0.001). Five common subfamilies (Pteromalidae: Asaphinae, Miscogasterinae, Pteromalinae; Eulophidae: Entedoninae, and Tetrastichinae) showed an obvious trend (p < 0.05) of the FL increasing with altitude, except for Eulophinae (p > 0.05) ( Figure S4), mainly due to the relatively large-sized Euplectrus, Cirrospilus, and Pnigalio generally being distributed at low altitudes ( Figure 4). However, after removing the data of these genera, the FL of Eulophinae also demonstrated an increasing trend with altitude (p < 0.05). The altitudinal trend of Miscogasterinae and Entedoninae is particularly significant, and the variation in Halticoptera (p < 0.001) and Entedon (p < 0.001) made the greatest respective contributions ( Figures S5 and S6). Based on the analysis of the genera with a sample size greater than 20 records, we found that most of the genera had a wide altitudinal range, between 200 and 4000 m; some alpine-dwelling genera, including Selderma, Callicarolynia, Stenomalina, Thinodytes, and Entedon, have a larger population at higher altitudes (≥3000 m), while some genera were only distributed at lower altitudes (≤1500 m), such as Norbanus, Philotrypesis, Dzhanokmenia, and Omphale. We used the GLS model to validate Bergmann's rule for 16 genera with widespread distribution and large populations (distributed in the altitudinal range of 200-4000 m and ≥ 30 records). The results indicated that the FL of 11 genera significantly increased with altitude among all 16 genera, while some genera demonstrated contrasting patterns (Table  3, Figures S5 and S6). The trend of Pteromalus was not significant (female: p = 0.244; male: p = 0.728), especially for its males, mostly due to their comparatively larger size in Pteromalidae, and the fact that there were also several larger species found at low altitudes. There was no significant correlation between the FL of Homoporus (female: p = 0.662; male: p = 0.707), Neochrysocharis (female: p = 0.411; male: p = 0.005), and Pediobius (female: p = 0.422; male: p = 0.281) and altitude in general. The FL of Chrysocharis (female: slope = −1.71; male: slope = −1.54) showed a reverse trend with increasing elevation. All the interactions between sex and altitude were not significant. We used the GLS model to validate Bergmann's rule for 16 genera with widespread distribution and large populations (distributed in the altitudinal range of 200-4000 m and ≥ 30 records). The results indicated that the FL of 11 genera significantly increased with altitude among all 16 genera, while some genera demonstrated contrasting patterns (Table 3, Figures S5 and S6). The trend of Pteromalus was not significant (female: p = 0.244; male: p = 0.728), especially for its males, mostly due to their comparatively larger size in Pteromalidae, and the fact that there were also several larger species found at low altitudes. There was no significant correlation between the FL of Homoporus (female: p = 0.662; male: p = 0.707), Neochrysocharis (female: p = 0.411; male: p = 0.005), and Pediobius (female: p = 0.422; male: p = 0.281) and altitude in general. The FL of Chrysocharis (female: slope = −1.71; male: slope = −1.54) showed a reverse trend with increasing elevation. All the interactions between sex and altitude were not significant. Table 3. Coefficient and ANOVA of the best fit GLS models on FL of 16 genera from Pteromalidae and Eulophidae (CE = coefficients, n = number of specimens). * denotes level of significance (* </ =0.05; ** </ =0.01; *** </ =0.001).

Discussion
Much of the previous research has documented the applicability of Bergmann's rule in multiple insect groups at the species or family level, while no general pattern has been detected in parasitoids, even for the same group, and there has been little in-depth related research on minute parasitic wasps. In this study, we examined the geographical patterns in two hyper-diverse groups of Chalcidoidea and analyzed the dominant environmental factors responsible for this variation. We found that the FL of Pteromalidae and Eulophidae increased with elevation in both sexes at the family level (a pattern that follows Bergmann's rule), with temperature and precipitation as the main environmental drivers.
Our results, based on the two different methods, might be explained by some of the hypotheses mentioned above. The average size of chalcids from these two taxa is larger at higher altitudes than at low altitudes. Temperature was obviously the most statistically explanatory variable in both taxa, explaining most of the variance in body size among other environmental factors in our study, and may play a more direct role. This phenomenon is also in line with numerous laboratory experiments, which show that Pnigalio soemius (Eulophidae: Eulophinae) raised at lower temperatures have larger FL [31], and are similar to Eretmocerus warrae (Hymenoptera: Aphelinidae). Larger insects usually have a smaller surface-area-to-volume ratio, which may help them maintain their thermal energy and reduce the water-loss rate [61], and also facilitate the storage of the energy needed to maintain normal activities [66]. Larger individuals are generally more fertile, which facilitates the reproduction of more offspring to maintain the population in the shorter reproductive periods in alpine regions. As stated in the starvation-resistance hypothesis, larger parasitoids are able to store more energy, such as lipids and glycogen, which are involved in egg production and other physiological mechanisms. The species' diversity and coverage of flowering plants in the alpine area were significantly lower than those at lower altitudes, meaning that parasitoids need to restore more energy to cope with the extreme environment; thus, larger species are more suitable for survival in the alpine region. Precipitation also showed a high correlation with the FL variation; it is generally inversely proportional to elevation and is supposed to be negatively correlated with FL, as with temperature, while our results are contrary to the above speculation. Numerous investigations' results about the effect of precipitation on body size are also controversial. Precipitation not only indirectly affects the growth of the host plants of parasitoids, and thus their population density, but also directly affects the water-loss rate of insects [67]. Therefore, precipitation may play an indirect role and interact with other environmental factors in complex ways. The geographical trend in size might be the result of interactions between various environmental factors. In addition to the factors considered above, the oxygen level (P O2 ) also declines significantly with altitude, and also has effects on the insect's body size, combined with temperature [68]. The host identity and host plants also influence the parasitoid's body size [69]. However, the various parasitoids analyzed in this study do not have specific host information, and the data on diet and DNA barcode also remain incomplete. This will also be a future research direction.
During the extensive field collection, we found brachypterous Pteromalidae and apterous Encyrtidae species in the alpine area (≥3500 m), while no similar species were found at low altitudes. We also found that the wing loading of Pteromalidae and Eulophidae was inversely proportional to the FL (Pteromalidae: cor = −0.639; Eulophidae: cor = −0.604) and decreased significantly (p < 0.001) with increasing altitude. A similar trend was also found in Drosophila [70], Bumbus [40], and moth [13], mainly because a lower wing loading is generally associated with a superior flight maneuverability, reduced energy being required for flight, and enhanced flight duration. Therefore, having this feature will help them live normally in the alpine environment. However, wing aspect ratio, also an important indicator of flight performance, did not significantly correlate with either body length or altitude in our research.
We found that the populations of Pahcyneuron, Halticoptera, and Selderma from the Pteromalidae, and Entedon, Neochrysocharis, and Diglyphus from the Eulophidae were highly dominant among other groups in the alpine region. Eleven out of sixteen genera followed Bergmann's rule, and only the genera Pteromalus, Homoporus, Neochrysocharis, and Pediobius showed no significant trend, while Charysocharis showed a transverse trend. Pteromalus is one of the large-sized genera of Pteromalidae; the average FL of the Pteromalus specimens in our data is 1.782 mm(female), 1.35 mm(male), and the average FL of all the Pteromalidae is 1.652 mm(female), 1.422 mm(male). The larger individuals in this genus are also abundant at a lower altitude, resulting in an insignificant FL trend with elevation. Two female specimens of Chrysocharis pubicornis from Norway, found at elevations of 65 and 63 m, are larger than the average FL of the genus (female: 1.777 > 1.03 mm, male: 1.128 > 0.997 mm); the FL of the genus showed no significant trend with altitude when we removed these pieces of data. Such contradictory trends may also be caused by the interaction between various factors or there being less diversity in alpine species, and need further information on these genera to resolve this problem. Through the statistical analysis of numerous species, we also discovered some species with a wide distribution across altitudes, including Asaphes vulgaris, Pachyneuron korlense, Pachyneuron aphidis, Diglyphus isaea, Diaulinopsis arenaria, and Neochrysocharis formosus, etc., which will provide objects for further study of the adaptation mechanism.
In our research, the geographical trends in the size of two abundant chalcids were analyzed along a wide range of altitudinal gradients. The most seminal result is that the alpine region tends to harbor larger chalcids (following Bergmann's rule) than the lower altitudes. Temperature and precipitation are the dominant environmental drivers of this variation, and the results indicate that temperature may be more directly related to body-size variation with altitudinal gradients than precipitation. Other hypotheses, such as resource availability, can also be used to explain the trend. The combination and interaction of these factors play an important role in the formation of body-size patterns. Additional information, such as genetic data (a test of whether body-size response to altitude shows a phylogenetic signal or there are positive selected genes) and more data on different species is needed to help further understand the mechanisms that cause body-size variations. Our results may provide a cornerstone for relevant future work.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/insects14010067/s1, Table S1: Overview of the numbers of analyzed specimens in Pteromalidae and Eulophidae. Figure S1: The correlation between seven morphological characteristics of both sexes in Pteromalidae(a) and Eulophidae(b); Figure S2: Frequency distribution histogram of the FL of Pteromalidae(a) and Eulophidae(b) from all elevation gradients; Figure S3: Histogram of the FL of Pteromalidae(a,b) and Eulophidae(c,d) from four elevation ranges; Figure  S4: The relationship between forewing length and elevation in six subfamilies from Pteromalidae (a-c) and Eulophidae (d-f); Figure S5: The relationship between forewing length and elevation in eight genera(n ≥ 30 and with wide altitudinal range) of Pteromalidae; Figure S6: The relationship between forewing length and elevation in eight genera (n ≥ 30 and with wide altitudinal range) of Eulophidae. File S1: Forewing length data of Pteromalidae and Eulophidae which used for GLM and GLS model in R analysis.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author, and also available in supplementary material here.