Fluctuating Asymmetry as a Method of Assessing Environmental Stress in Two Predatory Carabid Species within Mediterranean Agroecosystems

: Fluctuating asymmetry (FA) is used in assessing the e ﬀ ect of environmental stress on the development stability of individuals by measuring small random deviations from perfect bilateral symmetry. Here, we checked for FA on two predatory carabid beetles, Pterostichus melas and Poecilus koyi , in order to evaluate species response to agricultural practices within Mediterranean agroecosystems, as well as FA as a method. The samples were collected in vineyards and olive groves, both under integrated pest management (IPM) and ecological pest management (EPM), and in pristine habitats in the Mediterranean region of Croatia. Geometric morphometrics (GMMs) were used to analyze the pronotum and abdomen shape variations and left–right asymmetries of each population. In respect to the FA measurements, analyzed species responded di ﬀ erently, with P. koyi displaying a lower intensity of FA than P. melas . On the other hand, P. melas beetles from vineyards showed a higher intensity of FA compared with populations from pristine habitats and olive groves. Accordingly, FA pointed out olive groves as potentially less adverse habitats to predatory carabids, keeping in mind the di ﬀ erent levels of asymmetry between the two species. Our study singled out P. melas as a more suitable species for further research, in the e ﬀ ect that di ﬀ erent agricultural practices can have their impact on non-target invertebrates analyzed by measuring the FA.


Introduction
An increase in the application of plant protection products, due to the need for constant production growth and following the effects of climate change, can be viewed as one of the major factors that poses a threat to overall insect biodiversity [1,2]. Intensive use of pesticides has been proven to be harmful to non-target insects [3]. In several cases, neonicotinoids were known to cause sublethal effects to bees and other insects [4][5][6]. This effect of pesticides is not only harmful to pollinators, but also to predator species such as ground beetles, who are important for controlling changes in prey populations and thus in biocontrol in agroecosystems [7,8]. Furthermore, ground beetles are important bio indicators of different human-induced environmental stresses, including agronomic production as well as urbanization, forestry, landscape fragmentation, and other anthropogenic impacts [9][10][11][12][13].
Intensive management in agricultural habitats provides constant stressors to organisms from the developmental to adult stages, especially by adding a broad range of pesticides to the systems. Fluctuating asymmetry (FA) is a common technique that measures random, small deviations from bilateral symmetry in organisms and is used in the assessment of different levels of stress, environmental or genomic, that affect organisms during their development, known as developmental instability (DI) or developmental noise (DN) [3]. Sukhodolskaya et al. [14] showed that, in the case of ground beetles, there is a species-specific sensitivity to FA, even among closely related species. Furthermore, it was shown that some species cannot be used as indicators of environmentally induced stress for measuring FA. Asymmetry of ground beetle populations was used in various studies to measure the effects of different agricultural practices [3,15,16] and to estimate the effects of urbanization [17,18] and intensive landscape changes on forestry [19][20][21][22]. There are differences between various authors' approaches; Elek et al. [18] demonstrated that ground beetles, common along a Danish urbanization gradient, do not seem to indicate differences in habitat quality by their level of FA, and Labrie et al. [15] found no differences in FA between sites with different agricultural practices. On the other hand, Weller and Ganzhorn [17] showed that proximity to urban areas increases the FA in species which are susceptible to urbanization. Furthermore, Weller and Ganzhorn [17] noticed that eurytopic species are less prone to environmental changes and thus have smaller levels of FA. This led us to the conclusion that disturbances in the environment, spanning from urbanization changes to pesticide application in agriculture, sometimes do not cause changes to FA in a way that different levels of disturbances, such as higher levels of pesticides, can be distinguished. It is also important to highlight the difference in FA between sexes, as it is a factor that has been omitted in many studies, but could have an impact FA scores [18,[23][24][25].
Olive groves and vineyards are a vital part of agricultural production in the Mediterranean region. Copper, as a plant protection product, has been used in this area for centuries and is still widely used nowadays when numerous farmers are turning to organic production and ecologically based pest management (EPM), in which the use of copper is permitted. On the other hand, integrated pest management (IPM) and conventional management (CM) use a spectrum of chemically synthetic pesticides that include active compounds, targeting different specific groups of organisms [15,18,[23][24][25]. The level of impact on insects is different depending on the group of pesticides, but their joined effect must not be neglected. Some research shows that organic farms have more species at three trophic levels compared to conventional farming, but this did not improve the mortality of pest species or, thus, the ecosystem service [26]. In addition, the meta-analysis research of Bengtsson et al. [27] showed that there was a higher species richness in organic farming than IPM farming. Additionally, not only the species richness, but the total fitness of predators, such as ground beetles, were revealed as important when compared between pest management types [28][29][30][31]. On the other hand, some authors reported no significant increase in carabid beetle diversity among differently managed fields [32,33]. It is important to keep in mind that in the agroecosystem, there is a broad list of factors influencing the predator dynamics and abundance, and that these results are sometimes due to an overall effect [29].
Even though olive groves and vineyards have been an integral part of the Mediterranean landscape and its diversity for thousands of years, there is still a lack of studies regarding the impact of management of stress on predatory invertebrate fauna. Therefore, we focused our research on two predatory ground beetles, Poecilus koyi (Germar, 1824) and Pterostichus melas (Creutzer, 1799), that are common in the studied sites, located in the Mediterranean area of Croatia. The selected predatory carabids play an important role within agroecosystems in controlling their prey abundance, including potential pest species, and thus contribute to biocontrol of pest populations [34][35][36][37]. Using FA, we wanted to measure stress induced by different agricultural practices on these two species, assuming that a higher intensity of FA could point out the type of practice that was more adverse toward predatory carabid beetles. Besides checking for the presence of FA within populations across the sites under different agricultural practices, we questioned if the selected species, both co-occurring in arable land across southwest Europe [38][39][40], showed the same patterns of FA, or if one of them was more sensitive to environmental changes. In addition, the study aimed to analyze the variation in the body shape and size between sexes and populations from different agroecosystems.

Sampling Sites
Sampling took place within four agricultural sites in Zadar County in the Mediterranean part of Croatia: vineyards, one site with ecological pest management (EV) and one with integrated pest management (IV), and olive groves with ecological pest management (EO) and integrated pest management (IO) (Figure 1). They were in pristine habitats (C) with no histories of pesticide treatment, represented by typical Mediterranean maquis and garrigue, with one site in Zadar County in 2019 and one in nearby Šibenik-Knin County. The study sites have a long history in agricultural use, with EO and EV being under EPM for a decade or more. Detailed information about agricultural practices during 2018 in the selected sites is listed in Table 1. To estimate the intensity of the FA within the sites, two predatory carabid species, P. koyi and P. melas, were sampled, and we proceeded with further analyses ( Table 2). At the control sites, only P. melas individuals were sampled in a substantial number after pooling them from two control sites for the analyses ( Table 2). Information on the amount of pesticides applied on the studied sites in the year of sample collection is given in Figure 2.   Table 1. Number of treatments, including the number of pesticides added and soil processing during 2018 for four managed sites. Added pesticides have been grouped according to the main active compounds, those being synthetic, biological, and copper ones. Details are given below. Site abbreviations are as follows: one vineyard with EPM (EV) and one with IPM (IV), and one olive grove with EPM (EO) and one with IPM (IO).

Treatments
Sites

Collection of Specimens
Specimens were collected from mid-April to mid-July and from mid-September to mid-November in 2018, using pitfall traps and actively searching and collecting by hand. Plastic containers with a volume of 300 mL were placed in the ground and used as pitfall traps. Traps were covered with elevated, flat stone slabs. The distance between traps was approximately 2 m in all sites, as they were dug in under olive trees or under the stumps of vines, depending on the site. For the preservation of trapped arthropods, an aqueous salt solution was used. Traps were emptied every two to three weeks, and the material was transferred in 80% ethanol and analyzed in the laboratory.
From the total of the ground beetles sampled, two common and numerous species, P. koyi and P. melas, were present at almost all sites, and for further analyses, a substantial number of individuals were selected for FA measurements ( Table 2). In the case of the P. melas from EO, there were not enough specimens, even after additional sampling efforts, so this population was not included in further analyses. This applied to P. koyi from the control site (C) as well. The control sample of P. melas was a bulk sample, assembled from samples collected at two control sites. Samples were pooled and analyzed as one. Specimens of P. koyi were photographed using an Epson Perfection V600 photo scanner, while specimens of P. melas were photographed using a digital Nikon D60 camera with a Sigma macro objective. In both species, males and females were separated, based on sex combs on the first pair of legs or, in a case that feature was insufficient or absent, external genital parts were used, and they were analyzed separately. For all specimens of P. koyi and P. melas, the ventral part of the body was photographed, while the dorsal part of the body was imprinted in clay to keep the specimen horizontally levelled during the process.

Shape Analysis
A total of 16 landmarks on the ventral part of the body were used for P. koyi and P. melas (Figure 3). Suggested landmarks were digitized using tpsDig 2.31 software [41]. To be able to evaluate the significance of fluctuating asymmetry, landmarks were digitized twice in the case of P. koyi and three times in the case of P. melas. Statistical analyses were performed in MorphoJ [42] and STATISTICA 13 (Statistica Inc. TIBCO Software). After Procrustes superimposition was done in MorphoJ, the coordinates of each landmark were used for further analyses [43]. To analyze the measurement error, Procrustes analysis of variance (Procrustes ANOVA) was performed on combined datasets of repeated landmark digitization for both P. koyi and P. melas. To estimate the differences in size between specimens from different sites and sexes, ANOVA and post hoc unequal honest significant difference (HSD) was performed in STATISTICA 13. Procrustes ANOVA was performed for the centroid size and shape, with sites and sexes as extra effects. Canonical variate analysis (CVA) was used to determine the relationships between groups of variables in our data set and for better visualization of separated (discriminated) groups. CVA was performed using symmetric components of the shape variation. In order to assess statistically significant differences across the sites, a pairwise comparison after 10,000 permutations was performed using Mahalanobis and Procrustes distances. The presence of FA in the species was tested using Procrustes ANOVA by checking individual*side interactions and mean squares of individual*side (MS ind*side ) interactions, corrected for error variance, were used as a measure representing FA. The T-test for independent samples was used to statistically confirm the differences in Procrustes FA scores between the species. In addition, to check if the differences in Procrustes FA scores (p < 0.05) observed between the populations were significant, ANOVA and Tukey HSD post hoc tests were performed after transformation, as using the natural logarithm had normalized the data (Kolmogorov-Smirnov p > 0.05). Levene's test confirmed the null hypothesis of equal variances for tested groups and populations (p > 0.05).

Results
The Procrustes ANOVA, applied on repeated measures of individuals in order to assess the measurement error (ME), showed that the ME was negligible (MSindiv > ME, MSind*side > ME, with MS standing for mean squares) between sets of measurements for both P. melas and P. koyi (Table 3). Subsequently, the same analysis indicated significantly high variations for both shape and size between populations and sexes (p < 0.0001) of P. melas, but just for shape in the case of P. koyi. Females were larger than males at every site for P. melas (Figure 4). The post hoc test showed a significant difference in size between females from IV and the control (p = 0.00216) and females from EV and the control (p = 0.00427), in the case of P. melas. The P. koyi species showed no differences between males and females from different sites, or between the sexes within the same site.   Table 2.
A canonical variate analysis (CVA), applied to the symmetrical components of shape variation in P. melas, showed that the Mahalanobis distance significantly differed (p < 0.0001) between populations from studied sites, except for males from the control sites, which did not differ from other groups. Furthermore, there was a significant Mahalanobis distance between sexes within the same population (p < 0.0001) (Table 4, Figure 5). The results of the CVA with Mahalanobis and Procrustes distances (Tables 4 and 5, Figure 6) for P. koyi showed high morphological variability, especially in males, across the sites (p < 0.0001).  FA values (mean squares value of ind*side interaction, corrected for error variance) differed between the studied species, with P. melas displaying higher FA values (Figure 7, with the T-test on normalized Procrustes FA scores yielding t-value = 5.97, p < 0.05). Regarding P. melas, populations from all three agroecosystems had higher Procrustes FA scores than the control population, with the FA of the population from EV being the highest and differing significantly from other populations (ANOVA and Tukey HSD Post hoc, F = 11.01, p < 0.05). Individuals from the olive grove had lower FA values than those from the vineyards, but a bit higher than individuals from the control sites (Figure 7). There was no significant difference in FA values between the males and females of P. melas within or between the sites (Factorial ANOVA, p < 0.05). In P. koyi, on the other hand, males were more asymmetrical at every site except EV (Table 6, Figure 7), but there were no significant differences in FA values between its populations, neither for the site nor for sex and management (IPM vs. EPM) effects (Factorial ANOVA, p > 0.05). Regression of the Mahalanobis and Procrustes FAs with individual shapes was performed in order to see the magnitude in FA values, showing olive groves as indeed less influential on FA than vineyards (Figure 8).

Discussion
The FA, as a subtle deviation from bilateral symmetry, has long been of interest to researchers studying the effects that changes in environmental quality have on organisms [44]. In this study, we employed FA measurements on two predatory carabids, P. melas and P. koyi, which were common and relatively abundant in the study area [40]. The FA was detected in every population, regardless of species, site, or sex, but FA values differed between the species and between populations (Table 6, Figure 7). Besides differences in FA, Procrustes ANOVA and CVA revealed a significant variation in shape and size of the tested populations of both species between study sites, as well between sexes for P. melas.
Canonical variate analysis showed significant differences in shape between populations for both P. melas and P. koyi, with the exception of P. koyi females from the olive groves, which did not statistically differ from each other (Tables 4 and 5). Significant shape variation between sexes and among sites was noted for other species of beetles [18,[45][46][47][48][49]. In this study, P. melas females from both vineyards were larger than females from the control. The same observation was made for males, but it was not statistically significant. It is generally considered that environmental conditions cause variations in size between populations, whereas variations in shape reflect variations in the genetic composition [50,51]. In addition, it is possible that P. melas, a species with weak dispersal power, developed intraspecific variability in shape and size as a result of different habitat conditions at the study sites. The fact that the control group for P. melas was pooled from two different locations likely increased the shape variation. Alibert et al. [46] also found females to be larger than males, but there was no size difference between populations from different sites. Fitness in females is often size-related, and is profoundly influenced by conditions during larval development [52]. Such female-biased sexual dimorphism is commonly observed in invertebrate taxa, and it is hypothesized to be a result of a positive correlation between fecundity and female abdomen size [45].
Higher FA scores within P. melas populations from vineyards than those in olive orchards may indicate practices run within olive groves as potentially less adverse to predatory carabids. Contrary to our assumptions that IPM would produce higher levels of asymmetry within populations due to a broader spectrum of pesticides being applied [23,53,54], higher FA values were detected in both species from sites with ecological pest control practices, although for P. koyi, those differences were not as prominent as in P. melas, which appeared to be the less robust species. However, these results corroborate with findings from some previous research where FA levels did not differ between populations of ground beetle species among continuous forest and fragmentation habitats in a plantation area [22], or with different types of management in orchards [15].
A higher FA score for the P. melas from EV can be due to the frequency of mechanical tillage of the soil being the highest at this study site (Table 1). Carabids, as a part of the soil fauna, can be affected by mechanical disturbances of the soil, especially during their development. Benítez et al. [3] reported higher FA scores in populations within annual arable lands as more unstable than in a perennial agroecosystem. Additionally, these results are supported by the fact that this population has the highest density at the stated study site, as higher densities in populations enhance intraspecific competition, which reduces individual food availability and increases the FA [22,55]. Furthermore, Nattero et al. [6] showed that pesticide treatment levels of asymmetry decreased due to the higher mortality rates of less adapted individuals, often those with higher asymmetry levels. This potentially can explain higher FA values at sites with EPM, compared with those using IPM, meaning that anthropogenic factors (e.g., pesticide use and tillage) act as selective pressures that favor more symmetric individuals. On the other hand, individuals within control populations from pristine habitats develop without anthropogenic factors disrupting their developmental stability and, because of that, are the least asymmetric. Floate and Fox [56] offered a differential mortality hypothesis as an explanation for this phenomenon. After laboratory experiments conducted on dipterans, they observed that only individuals that survived stress could be captured later, which means our sample consisted of less affected individuals. Several previous studies showed a correlation between survival probability and the level of FA in insects [57][58][59]. Keeping in mind that field conditions are even more variable than laboratory conditions, more research is needed to establish FA as an indicator of the harmful effects of pesticides on predatory organisms. Research on two different populations of P. melas in Croatian agroecosystems did not reveal any significant difference between FA values when compared with the control [3], and [60] did not detect any relationship between FA and heavy metal pollution in populations of the wolf spider Pirata piraticus. Furthermore, more common species can be less susceptible to stress in their environment [17], which could also explain why there was little difference in the FA values between P. koyi populations, as it was one of the dominant species at our research sites. The same observation was made for rodents from two different farms, with the conclusion that the effects of the farming practice would affect generalist species less than it would affect specialists [5]. P. melas, as a predator from natural forest ecosystems, was able to adapt to conditions in Mediterranean areas managed for agricultural purposes, such as olive groves [13]. It is considered that this adaptation was boosted by the turnover of specimen from surrounding natural environments to agricultural lands. This could also be the reason why P. melas was more susceptible to asymmetry changes, due to agricultural practices in areas where this turnover was made more difficult. Furthermore, the reason for this difference in asymmetry levels between the two species could be due to the fact that P. koyi is more adapted to different types of agricultural lands, where often this species is among the most abundant ones. Different responses between species in terms of FA, and species sensitivity and adaptation to agricultural practices, may be strongly connected to its historical presence in the natural environment surrounding arable land [17,61]. Mazeed [62] demonstrated that a honeybee species not native to the area had a higher FA than the native species, meaning that the geographical location and adaptation to different environmental conditions, to which each organism is subjected, may affect the degree of asymmetry.
Differences in FA between sexes were noted in previous studies [18,23,24], where males of studied species had been regularly less asymmetric compared with females. In our study, P. koyi males were more asymmetrical at every site except EV, while no differences between sexes were noted in P. melas with the exception of the control group, where the males were a bit more asymmetrical. Since the FA in P. koyi was much lower than in P. melas, and the lowest FA score was recorded for the P. melas control group, we believe that the differences in FA between sexes could be ignored in this case. This small deviation in FA between sexes within the control group could be due to pulling P. melas from two different pristine locations.

Conclusions
As expected, the lowest FA intensity was found in the control population of P. melas, which goes in favor of FA as a method to be used in assessing certain environmental stresses caused by agricultural practices. However, the results indicate that different species do not have the same response in respect to FA. P. koyi showed no patterns in FA between populations from different agricultural sites, while the populations of P. melas from olive groves displayed a lower FA than populations from vineyards, and the FA values were closer to the control group from the pristine habitat. Thus, we recommend P. melas as a test species for future studies on environmental stress using FA. Furthermore, FA singled out olive groves as a potentially less adverse habitat to predatory carabids when compared with vineyards, regardless of the pest management (IPM or EPM). In future studies, the FA of P. melas can be used as a method to indicate agricultural practices that are closer to natural.