*2.3. Microsatellite Analysis*

For comparison of different honey bee colonies on the autosomal level, we have 14 microsatellite loci described in [61]. The choice of loci used in the analysis was made

according to the most frequent microsatellite loci used in a number of different studies dealing with the variability of these genetic markers. The selected loci, primer pairs used for amplification, and the corresponding annealing temperatures are presented in Table S1. The microsatellite loci were amplified in PCR reactions in which forward primers were labeled with a fluorescent dye (Table S1). PCR was performed in four reactions that differed in annealing temperature (Table S1) using the following program: one cycle of initial denaturation at 94 ◦C for 5 min, after which there were 30 cycles of 35 s at 94 ◦C, 35 s at annealing temperature (Table S1) and 35 s at 72 ◦C. The final elongation step was performed at 72 ◦C for one hour. Loci were amplified in a MiniAmp Plus Thermal Cycler (Applied Biosystems, ThermoFisher Scientific) in four multiplex reactions. The amplification was carried out in a volume of 20 µL with the following final concentrations of reaction components: 1 × Taq Buffer with (NH4)2SO4, 2.5 mM MgCl2, 0.8 mM dNTP mix, 1 U of reverse Taq polymerase (all components were produced by Thermo Fisher Scientific, EU), and 5 pmol of each forward and reverse primer. For the amplification of the microsatellite loci, 1.5–1.9 ng of DNA was used. To verify the reliability of the data, 10% of samples were reamplified for the second time.

To use data from the study [58] we performed calibration by reanalyzing 10 samples and 9 microsatellite loci from this data set. The DNA of the same worker honey bees used in [58] was processed in the same way as the samples from Serbia.

### *2.4. Fragment Analysis*

For fragment analysis, the first and second multiplex reactions were multipooled, and the third and fourth multiplex reactions were multipooled. All reactions were mixed in equal volumes and plated as one reaction in a volume of 1 µL. Each amplification mix contained seven different loci. GeneScan 600 LIZ size standard was used to score alleles (Applied Biosystems, Warrington, UK). Fragment analysis was performed on the 3130 Genetic Analyzer (Applied Biosystems, UK). Data were analyzed using Gene Mapper Software (Life Technologies, Foster City, CA, USA).

#### *2.5. Statistical Analyses*

The standard parameters of genetic diversity for microsatellite loci (number of alleles, allelic size range, average gene diversity over loci, number of alleles based on a minimal sample size (obtained by rarefaction), number of private alleles based on a minimal sample size (obtained by rarefaction), observed (HO) and expected (HE) heterozygosity, random match probability (RMP), and the mean number of pairwise differences (MPD)) were calculated using Arlequin ver. 3.5.2.2 software [62] and HP-Rare 1.1 [63]. The RMP parameter is used to express the probability that two randomly sampled individuals from a population have a matching genotype and is calculated as the sum of square frequencies [64]. MPD is a parameter that represents the measure of differences between all pairs of haplotypes in the sample. Arlequin ver. 3.5.2.2 software was also used to assess genetic differentiation among populations by analysis of molecular variances (AMOVA) and to estimate the pairwise population and overall *FST* and *FIS* values. The statistical significance of all performed tests was assessed with 10,000 permutations. The matrix of pairwise population *FST* values was visualized using a multidimensional scaling method (nonmetric MDS) implemented in the PAST 3.25 software [65], and the R functions connected with Arlequin ver. 3.5.2.2 software. The Hardy–Weinberg equilibrium was tested using Arlequin ver. 3.5.2.2 software with 1,000,000 steps in MC and 100,000 dememorization steps. To correct the probabilities when multiple tests were performed simultaneously, we performed a sequential Bonferroni test for the Hardy–Weinberg equilibrium. The linkage disequilibrium between the pairs of loci was estimated using the likelihood ratio test in Arlequin ver. 3.5.2.2 software with 10,000 steps in MC and 10,000 dememorization steps.

The number of genetic clusters represented in the sample was estimated with STRUC-TURE v 2.3.4 software [66–68]. For the analysis, the admixture model was used with a burn length of 10,000 and a Markov chain Monte Carlo (MCMC) of 100,000 randomizations. The range of the possible number of clusters (K) was from 1 to 10, with a series of 10 runs for each K. The results obtained by STRUCTURE were analyzed by the STRUCTURE harvester [69]. To detect the number of K groups that best fit the data set, this software used results generated by the STRUCTURE software to create a plot of the mean likelihood value per K value and calculated the highest value of the second-order rate of change (∆ K) using the of Evanno method [70]. The model choice criterion, LnP (D), implemented in the STRUCTURE which detects the true K as an estimate of the posterior probability of the data for a given K, was evaluated as well. The most likely scenario was chosen and used to graphically plot both the individuals and populations analyzed.

The observed distances among samples are presented using discriminant analysis of principal components (DAPC) [71]. This method consists of performing the linear discriminant analysis (LDA) on the principal components analysis' (PCA) transformed matrix. In the case of samples from Serbia, only LDA was performed on the first 32 PCs and in the case of all analyzed populations on the first 55 PCs which cumulatively conserve 98.9% of the total variance. The number of retained PCs was estimated using randomly repeated cross-validation (100 iterations), which consisted of performing DAPC on 90% randomly sampled training set observations (stratified sampling was used so the training set consisted of 90% of the observations from each population) after retaining 10–183 PCs and using the obtained model to predict the groups (populations) in the remaining 10% of samples (test set). Average prediction success per group was used as a metric. Additionally, the PCA transformed matrix (all 183 PCs) was used to find the optimal number of clusters using Ward's method [72]. We tested 2–50 clusters, and the optimal number of clusters was chosen using BIC statistics using the "diffNgroup" method. This method uses Ward's clustering method to split the differences between successive values of the BIC summary statistic into two groups to differentiate sharp decreases from mild decreases or increases. The retained K was the one before the first group switch. Thus, estimated clusters of samples were compared with the a priori defined populations.

#### **3. Results**

### *3.1. PCR-RFLP*

The size of the PCR-amplified COI segment for the RFLP analysis was 1029 bp. Digestion with both *Nco*I and *Sty*I did not show a restriction pattern characteristic for mtDNA lineage found in *A. m. macedonica*. Since no restriction sites were observed after RFLP analysis, we presume that all individuals in our sample belong to *A. m. carnica* [46,61].

#### *3.2. Genetic Diversity Analysis*

#### 3.2.1. Genetic Diversity Analysis for 14 Microsatellite Loci in the Serbian Sample

The standard diversity parameters for the sampled localities in Serbia for all 14 analyzed microsatellite loci are presented in Table 1 and Table S2. The average numbers of alleles' observed heterozygosity and average gene diversity over loci were the highest in L, the lowest values for these parameters were found in T, Vr, and FG, respectively. Considering the mean number of private alleles based on minimal sample size, the highest number was observed in L, but values in other analyzed localities were in the same range. The observed heterozygosity was generally lower than expected and *FIS* values varied between −0.04 in T and 0.19 in Vr (Table S3). It is interesting to note that the departure from the Hardy–Weinberg equilibrium coincided with significant heterozygote deficiencies, especially for locus A43 in all localities except FG. Furthermore, in all localities except L, observed heterozygosity before Bonferroni corrections for selectively adaptive locus Ap249 was significantly lower than expected (Table S2), and after correction it remained significant for T, FG, S, and DP. Linkage disequilibrium was also observed for some pairs of loci, mostly prominent for two selectively adaptive loci (Ap249 and B124) in three northern and one southern locality (Table S4).


**Table 1.** Standard diversity parameters for sampled localities in Serbia for 14 microsatellite loci.

N—number of genotyped individuals, Na—the average number of alleles, Ar—allelic size range, Agd—average gene diversity over loci, Ho—observed heterozygosity, He—expected heterozygosity, RMP—random match probability, MPD—mean number of pairwise differences, G–W—Garza–Williamson index, Ar8—number of alleles based on a minimal sample size of 8 diploid individuals, and Apr—number of private alleles based on a sample of 8 diploid individuals. Part of the results are published in Proceedings of the The 1st International Electronic Conference on Entomology session Apiculture and Pollinators, 1–15 July 2021, MDPI: Basel, Switzerland, doi: 10.3390/IECE-10720.

#### 3.2.2. Population Genetic Analysis for Nine Microsatellite Loci in All Sampled Localities

The standard diversity parameters for all sampled localities for nine analyzed microsatellite loci are presented in Tables 2 and S5. The average number of alleles was the highest in Hungary and the lowest in the Polish sampled site, Wroclaw. The average gene diversity over loci was the highest in Poland and Spain and the lowest in Serbian populations. Heterozygosity excess was observed in all populations except those in Serbia, both for all analyzed loci and individual loci per population as well as for *FIS* values (Tables S5 and S6). Linkage disequilibrium analysis was also performed for this set of data and the results are presented in Table S7.

**Table 2.** Standard diversity parameters for all analyzed localities for 9 microsatellite loci.


N—number of genotyped individuals, Na—the average number of alleles, Ar—allelic size range, Agd—average gene diversity over loci, Ho—observed heterozygosity, He—expected heterozygosity, MPD—mean number of pairwise differences, RMP—random match probability, G–W—Garza–Williamson index, Ar10—number of alleles based on a minimal sample size of 10 diploid individuals, and Apr—number of private alleles based on a sample of 10 diploid individuals.

1

### *3.3. Population Structure*

3.3.1. Population Structure Based on 14 Microsatellite Loci in the Serbian Sample

The average number of pairwise differences between and within Serbian localities together with *Nei's* distances is visualized in Figure 2a and pairwise *FST* is visualized in Figure 2b.

**Figure 2.** Matrixes of the average number of pairwise *Nei's* (**a**) and *FST* (**b**) distances based on the analysis of 14 microsatellite loci for localities in Serbia. (**a**) The average number of pairwise differences between populations is presented above diagonally, the average number of pairwise differences within the population is presented diagonally, and *Nei's* distances are presented below diagonally. (**b**) Statistically significant *FST* values are marked with an asterisk (\*). Part of the results are published in Proceedings of the The 1st International Electronic Conference on Entomology session Apiculture and Pollinators, 1–15 July 2021, MDPI: Basel, Switzerland, doi: 10.3390/IECE-10720.

Differences between some pairs of localities were consistent in all analyses with the Vr locality showing statistically significant pairwise differences with all analyzed localities. Statistically significant differences between localities were observed for pairs of the south (L-SP and L-T), all north and south/north (L-FG, L-DP, L-S, V-S, SP-S, and SP-DP) comparisons (Figure 1, Table S8). Overall, the north showed greater population differences than the southern regions, and there is no clear pattern of differentiation between the geographical regions.

The AMOVA performed across all 14 loci showed a low but significant value of genetic variance between localities (0.047), with 1.42% of the genetic variance being attributed to the variation among localities (Table 3). When localities were grouped according to their geographical region, the percentage of variation was higher within regions than among them. Additionally, when localities were grouped according to their region, the percentage of variation among localities within regions remained statistically significant, while differentiation between geographical regions could not be observed (Table S9).

The results of the analysis performed by the DAPC method are shown in Figure 3, Supplementary Figure S1, and the visualization by the MDS plot shows the positioning of populations in two dimensions (Figure 4). Although individuals from T, SP, L, and DP tend to cluster separately from others, DAPC analysis showed that individuals from geographically remote localities are grouped in cluster overlaps, indicating similarity between them. Moreover, assignment to the previous predesigned group was relatively

low, with *p* ranging 0.2–0.3, indicating admixture. The MDS plot placed localities separately from each other, which is in correlation with AMOVA, suggesting the presence of distinct genetic variability in all analyzed localities. However, there is no clear grouping of localities according to their geographical region, which is also in concordance with AMOVA.

**Table 3.** AMOVA results when all localities in Serbia were analyzed without grouping.


d.f.—degrees of freedom, SS—the sum of squares, and *p*-statistical significance (statistically significant values are in bold).

**Figure 3.** Discriminant analysis of principal components. The first and second linear discriminants are presented in the plot.

The STRUCTURE Harvester showed that *K* = 2 is the most likely scenario (Supplementary Figure S2). The same number of clusters was inferred by LnP (D) analysis (Supplementary Figure S3). Both clusters inferred by STRUCTURE are equally distributed in all sampled localities. Additionally, the number of clusters inferred with the DAPC method was 8, with mixed distribution across localities (Supplementary Figure S4).

3.3.2. Population Structure Based on Nine Microsatellite Loci in All Sampled Localities

The average number of pairwise differences between and within localities grouped according to their geographical region and subspecies together with *Nei's* distances is visualized in Figure 5a, and pairwise *FST* is visualized in Figure 5b. As expected, individuals from Spain were shown to be the most separated from others, but the separation between region and subspecies can also be observed since statistically significant *FST* values were obtained (Table S10). The same conclusion can be inferred when all sampled localities were compared separately (Supplementary Figures S5 and S6, and Table S11). It is interesting to note that some differences between Serbian localities disappear but that there is a clear distinction between Serbian localities and other analyzed localities and subspecies. Moreover, very low but statistically significant *FST* value was detected between southern and northern Serbia.

**Figure 4.** Non-metric multidimensional scaling plot of *FST* distances between localities in Serbia. The goodness of fit is expressed with the stress value which is 0.1519 for this data set. Population pairwise *FST* values are presented in Table S8.

The AMOVA performed across nine loci showed a high and significant value of genetic variance among localities (0.28) with 10.79% of the genetic variance being attributed to the variation among the localities (Table 4). A negative value of differences among individuals within localities indicates that individuals in any given sampled population are mostly uniform and closely related to each other. Additionally, when AMOVA was performed with a different grouping of localities and subspecies, differences among geographical regions and subspecies remained significant, indicating regional differentiation that reflects subspecies and geographical distribution (Table S12).

The results of AMOVA were further corroborated by DAPC analysis (Figure 6 and Supplementary Figure S7) and the positioning of the populations in two dimensions in the MDS plot (Figure 7). When localities were grouped according to geographical region and subspecies, clear differentiation could be observed. As expected, Spain's population is the most separated from the others. Buckfast individuals from Hungary are closer to Italian individuals than Hungarian ones, and the Hungarian population is relatively homogeneous as previously reported. Serbian localities were in a cluster overlap and separated from other analyzed populations. The alternative grouping of localities and subspecies does not change the relative relations among the analyzed localities; localities from the same geographical region and subspecies were always clustered together and separated from others (Supplementary Figures S5–S7).

**Figure 5.** Matrixes of the average number of pairwise *Nei's* (**a**) and *FST* (**b**) distances based on the analysis of 9 microsatellite loci when all localities were grouped according to geographical region and subspecies. (**a**) The average number of pairwise differences between populations is presented above diagonally, the average number of pairwise differences within the population is presented diagonally, and *Nei's* distances are presented below diagonally. (**b**) Statistically significant *FST* values are marked with an asterisk (\*). **Figure 5.** Matrixes of the average number of pairwise *Nei's* (**a**) and *FST* (**b**) distances based on the analysis of 9 microsatellite loci when all localities were grouped according to geographical region and subspecies. (**a**) The average number of pairwise differences between populations is presented above diagonally, the average number of pairwise differences within the population is presented diagonally, and *Nei's* distances are presented below diagonally. (**b**) Statistically significant *FST* values are marked with an asterisk (\*).



**Table 4.** AMOVA results when all localities were analyzed without grouping. d.f.—degrees of freedom, SS—the sum of squares, and *p*-statistical significance (statistically significant values are in bold).

**Source of Variation d.f. SS Variance Components Percentage of Variation** Among localities 16 286.781 0.2801 **10.79 (***p* **= 0.00)**  Among individuals within localities 557 1034.257 −0.4601 −17.72 Within individuals 574 1594.000 2.7770 106.93 Total 1147 2915.037 2.5970 d.f.—degrees of freedom, SS—the sum of squares, and *p*-statistical significance (statistically significant values are in bold). The STRUCTURE Harvester showed that *K* = 8 is the most likely scenario (Figure 8), since the LnP (D) showed that *K* = 8 best fits the data even though ∆ K suggested *K* = 2 has the highest probability (Figure S8). In our data, *K* = 8 gives the most plausible distribution of inferred genetic clusters, which were specifically distributed among the individuals in the populations originating from different geographical regions or subspecies. Furthermore, the number of clusters inferred with the DAPC method was four, with specific distribution of clusters across sampled geographical regions (Supplementary Figure S9).

The results of AMOVA were further corroborated by DAPC analysis (Figure 6 and Based on all analyses it can be concluded that strong geographical differentiation exists between analyzed geographical regions and their corresponding subspecies.

Supplementary Figure S7) and the positioning of the populations in two dimensions in the MDS plot (Figure 7). When localities were grouped according to geographical region and subspecies, clear differentiation could be observed. As expected, Spain's population is the most separated from the others. Buckfast individuals from Hungary are closer to Italian individuals than Hungarian ones, and the Hungarian population is relatively homogeneous as previously reported. Serbian localities were in a cluster overlap and separated from other analyzed populations. The alternative grouping of localities and subspecies does not change the relative relations among the analyzed localities; localities from

**Figure 6.** Discriminant analysis of principal components when all localities were grouped according to their geographical region and subspecies. The first and second linear discriminants are presented in the plot.

**Figure 7.** Non-metric multidimensional scaling plot of *FST* distances between analyzed geographical regions and subspecies. The goodness of fit is expressed with the stress value which is 0.0426 for this data set. Population pairwise *FST* values are presented in Table S10.

**Figure 8.** (**a**) L(K) mean for the assumed number of genetic clusters. (**b**) Proportions of inferred STRUCTURE clusters (*K* = 4 and *K* = 8). (**c**) Proportions of the inferred STRUCTURE clusters (*K* = 4 and *K* = 8) from the individuals. 1—Hungary, 2—Spain, 3–7—Poland (3—August forest, 4—Krakow, 5—Bialowieza, 6—Wrocław, and 7—Siedlice), \* 8—Buckfast lineage from Hungary, 9—Italy, 10–13— Southern Serbia (10—Leskovac, 11—Vlasina, 12—Tromeda, and 13—Stara Planina), and 14–17— ¯ Northern Serbia (14—Fruška Gora, 15—Subotica, 16—Deliblatska pešˇcara, and 17—Vršac).

#### **4. Discussion**

The modernization of beekeeping practices and rapidly growing numbers of beehives in Serbia have invoked the need to re-study previously described genetic variability in the Serbian honey bee population [23,24,26,44,45,53]. Most of the previous genetic studies on the Serbian honey bee, even the most recently published ones, are based on samples from the first decade of the 21st century. Since then, significant changes in beekeeping practices together with stricter implementation of Serbian legislation and an increased number of beehives have led to changes in genetic variability in Serbian honey bee populations, as suggested in [46,61]. Therefore, we examined 14 microsatellite loci in Serbian worker bees from eight different localities to further shed light on the current status of the genetic diversity of Serbian honey bees. Furthermore, we compared nine microsatellite loci in our sample with previously published samples from Hungary, Poland, Spain, and Italy [58] to infer broader genetic relations between different *A. mellifera* populations and subspecies. Our results suggest that the Serbian honey bee population is relatively homogenous with preserved local differences and separated from the other populations analyzed.

The parameters of genetic diversity in Serbian localities are relatively high but lower than in other analyzed localities. Moreover, reference localities showed significant heterozygosity excess, while in Serbian localities Ho was in line with He, and for some loci, heterozygote deficiency may have been observed. These results, together with the G–W

index, indicate that the Serbian honey bee population did not experience recent bottleneck events. For some loci, departures from Hardy–Weinberg and linkage disequilibrium were observed, which may be an indication of recent gene flow from other subspecies or populations [73].

Although significant *FST* values were obtained between some pairs of localities, there is no clear pattern that indicates a south/north geographical distribution of microsatellite loci in the Serbian honey bee population. The results of AMOVA analysis suggest that grouping according to region may indicate some geographical distribution since the percentage of variation among localities within groups slightly decreases, but the value is not statistically significant (0.32, *p* = 0.119). Together with the equal distribution of two clusters inferred by STRUCTURE analysis in all localities and cluster overlapping inferred by DAPC analysis, the presented results indicate population admixture and a relatively homogenous population. However, local differences are still preserved since significant *FST* values can be observed between pairs of localities, and some differentiation may be observed according to the position of population in DAPC and MDS landscapes.

Our results could not confirm the presence of *A. m. macedonica* and north/south differences between individuals from different parts of the country previously reported for the Serbian honey bee population [23,26,53] but are in concordance with recently published work on the genetic diversity of the mtDNA tRNAleu-cox2 for the same sample [46]. Both uniparental and biparental markers showed that although local specific genetic variants and weak regional differences can be observed, previously reported regional differences indicative of subspecies distribution could not be confirmed. The formation of an admixture population may be one of the reasons behind the presented results. Extensive hybridization between *A. m. carnica* and *A. m. macedonica* subspecies in the central part of Serbian territory was previously described [23,26], and it is possible that the hybridization zone expanded reflecting recent changes in beekeeping practices as was shown by mtDNA data. The absence of north/south regional differentiation may also be partially attributed to the intensification of migratory beekeeping, since apiaries from the south are transported to the north during the flowering season of agricultural plants. As there is no human control of mating between individuals from different apiaries, when migratory apiaries return to the region where the stationary apiaries sampled in this study are located, admixture may be propelled. The Serbian leading beekeeper organization strongly encourages strict implementation of Serbian legislation that only *A. m. carnica* subspecies can be present in apiaries which, with a growing number of *A. m. carnica* queen manufacturers, may also contribute to the observed loss of *A. m. macedonica* and the admixture of the Serbian honey bee population.

Population structure analysis showed that each geographical region and each corresponding subspecies are separated. *A. m. iberica* from Spain is represented by its own cluster in STRUCTURE analysis as well as clearly separated in MDS and DAPC plots. The Italian *A. m. ligustica* is represented by its own structure cluster which is present in the Hungarian Buckfast sample as expected. The same result was obtained from DAPC and MDS analysis. The Polish populations of *A. m. carnica* and *A. m. mellifera* were located close to each other but still separated. The Hungarian population is well separated from the other populations studied, and although relatively homogeneous, some *A. m. ligustica* introgression may be observed, as previously reported by [58]. Serbian populations are well separated from others with significant overlap between individuals from the south and north, although very low but significant *FST* value can be observed. Structure analysis showed weak but still detectable introgression of *A. m. ligustica* alleles, which is in concordance with our field data that some illegal importation of Italian bees occurred in the past, since this subspecies has been one of beekeepers' favorites. Two distinct clusters can be observed in the Serbian honey bee population, and they are almost equally distributed among localities, suggesting population admixture.

Our results suggest that, as already shown in many studies [74], geographical distance together with environmental factors maintain the specific genetic diversity of *A. mellifera* subspecies within any given geographical region. However, this genetic diversity is under constant anthropogenic influence due to the modernization of beekeeping practices, such as migratory beekeeping, importation of foreign queens, and even legal practices [4,33,35,37]. Serbia is the natural area of contact between warmer climates preferring *A. m. macedonica* and colder ones preferring *A. m. carnica* and, although relief and ecological differences exist between these two parts of the country, climate conditions are favorable for both subspecies, and the main reason between their distinct distribution may lie in isolation by distance. The distances between the southern and northern parts of Serbia may be too great for bees, but for beekeepers they are rather small and they readily travel 500 km in the flowering season for different plants. Together with a vast increase in the number of beehives and beekeepers in the past decade [54] and legislation that specifically allows breeding of a single subspecies, it is not surprising that the previous composition of the diversity of honey bees in Serbia has changed. However, specific local genetic variability may still be retained since differences between analyzed localities can be observed.

Unfortunately, socioeconomic interests and beekeeper preferences for more productive and gentler individuals have often taken precedence over the conservation of locally native subspecies [27,37,39], leading to the predominance of admixture populations in humandominated areas [4]. Our results suggest that this scenario happened in Serbian honey bee populations and that for the above-mentioned reasons Serbia now harbors a distinct hybrid honey bee population. Further analysis that will include honey bee populations from eastern and western parts of Serbia are needed in order to better understand the pattern of genetic variability of managed honey bees in Serbia, so that the best managing strategies, with the goal of preserving the existing genetic diversity, can be implemented.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/insects13020180/s1, Figure S1: Discriminant analysis of principal components. The first three linear discriminants are presented in the plot; Figure S2: (a) Delta *K* values for the assumed number of genetic clusters. (b) Proportions of inferred STRUCTURE clusters (*K* = 2). (c) Proportions of the inferred STRUCTURE clusters (*K* = 2) from the individual worker bees; Figure S3: Ln values of probability for the assumed number of genetic clusters; Figure S4: Distribution of clusters according to the DAPC method and inffered number of 8 clusters; Figure S5: Non-metric multidimensional scaling plot of *FST* distances between 8 localities in Serbia (*A. m. carnica*) and other analysed *A. mellifera* populations from Hungary (*A. m. carnica* and Buckfast), Poland (Krakow and Wroclaw–*A. m. carnica*; August forest, Bialowieza and Siedlice–*A. m. mellifera*) Spain (*A. m. iberica*) and Italy (*A. m. ligustica*); Figure S6: Non-metric multidimensional scaling plot of *FST* distances between 8 localities in Serbia (*A. m. carnica*) and other analysed *A. mellifera* populations from Hungary (*A. m. carnica* and Buckfast), Poland (Krakow and Wroclaw–*A. m. carnica*; August forest, Bialowieza and Siedlice–*A. m. mellifera*) Spain (*A. m. iberica*) and Italy (*A. m. ligustica*); Figure S7: Discriminant analysis of principal components. The first three linear discriminants are presented in the plot; Figure S8: Delta *K* values for the assumed number of genetic clusters; Figure S9: Distribution of clusters according to the DAPC method and inffered number of 4 clusters; Table S1: List of loci used in genetic analyses, primers used for their amplification, fluorescent dyes used for tagging the forward primers, annealing temperature (Tm) used in each reaction for the amplification of specific microsatellite loci with the combination of loci amplified together in multiplex reactions I, II, III and IV; Table S2: Parameters of genetic diversity calculated per locus per populations including Garza Williamson index, Hardy Weinberg equilibrium (1,001,000 steps done)– Serbian localities; Table S3: Population specific Fis indices for 10,100 permutations (Serbian localities); Table S4: Tables of significant linkage disequilibrium (Serbian localities); Table S5: Parameters of genetic diversity calculated per locus per populations including Garza Williamson index, Hardy Weinberg equilibrium (1,001,000 steps done) (all populations); Table S6: Population specific Fis indices for 10,100 permutations (all populations); Table S7: Tables of significant linkage disequilibrium (all populations); Table S8. Pairwise population Fst (below diagonal) and Fst *p* values (above diagonal) between the populations based on the variability of 14 microsatellite loci found in 8 different localities of *Apis mellifera carnica* in Serbia; Table S9: Outcomes of AMOVA analysis based on the variability of 14 microsatelite loci when population sample from Serbia was grouped according to the geographical region: North vs South; Table S10: Pairwise population Fst (below diagonal) and Fst *p* values (above

diagonal) between the populations based on the variability of 9 microsatellite loci when Serbian localities were grouped according to their geographical region and other *A. mellifera* populations from: Hungary (*A. m. carnica* and Buckfast), Poland (*A. m. carnica* and *A. m. mellifera*), Italy (*A. m. ligustica*) and Spain (*A. m. iberica*); Table S11: Pairwise population Fst (below diagonal) and Fst *p* values (above diagonal) between the populations based on the variability of 9 microsatellite loci found in 8 different localities of *Apis mellifera carnica* from Serbia and other *A. mellifera* populations from: Hungary (*A. m. carnica* and Buckfast), Poland (*A. m. carnica* and *A. m. mellifera*), Italy (*A. m. ligustica*) and Spain (*A. m. iberica*); Table S12: Outcomes of AMOVA analysis based on the variability of 9 microsatelite loci when population sample from Serbia was compared with other analysed populations.

**Author Contributions:** Conceptualization, M.T., L.S. and S.D.; methodology, M.T., A.P., K.E., P.E., M.M., V.T. and S.D.; software, M.T. and S.D.; validation, M.T., A.P., K.E., P.E., L.S., M.M., V.T. and S.D.; formal analysis, M.T., A.P., K.E., P.E., M.M., V.T. and S.D.; investigation, M.T., A.P., K.E., P.E. and S.D.; resources, M.T., A.P., K.E., P.E., S.K., A.O. and S.D.; data curation, M.T. and S.D.; writing—original draft preparation, M.T. and S.D.; writing—review and editing, M.T., A.P., K.E., P.E., L.S., S.K., A.O. and S.D.; visualization, M.T. and S.D.; supervision, S.D. and L.S.; project administration, S.D.; funding acquisition, S.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Science Fund of the Republic of Serbia, PROMIS, #GRANT No 6066205, SERBHIWE for M.T., A.P., P.E., K.E., L.S., S.D., and the Ministry of Education, Science and Technological Development of the Republic of Serbia, grant number 451-03-9/2021-14/ 200007 for M.T., P.E., A.P., K.E. and S.D.

**Data Availability Statement:** Data are available on request from the corresponding author.

**Acknowledgments:** We are thankful to three anonymous reviewers and the Academic Editor whose comments significantly improved our work. Part of the data from this paper were previously published in Proceedings of the 1st International Electronic Conference on Entomology, 1–15 July 2021, MDPI: Basel, Switzerland, doi:10.3390/IECE-10720, Tanaskovi´c, M.; Patenkovi´c, A.; Eri´c, K.; Eri´c, P.; Stanisavljevi´c, L.; Davidovi´c, S., 2021. Microsatellite analysis of Apis mellifera from northern and southern parts of Serbia. In *The 1st International Electronic Conference on Entomology session Apiculture and Pollinators; 1–15 July 2021*. Basel: MDPI.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**


**Paula Souto 1,2,\* ,†, Gonçalo Abraços-Duarte 1,2,† , Elsa Borges da Silva 1,3 and Elisabete Figueiredo 1,2,\***


**Simple Summary:** Omnivorous predators, such as some mirids, are important biological control agents in several vegetable crops since they are generalists and can survive in the crop in the absence of prey. *Nesidiocoris tenuis* is a mirid used worldwide and its phytophagy is well known, which is not the case for the Palearctic *Dicyphus cerastii.* To use the latter in biological control it is crucial to evaluate the damage it causes to plants. We compared these two mirid species, under laboratory and semi-field conditions, regarding the damage they caused to plants and fruit, and their location on the plant versus on the fruit. Both species produced plant damage (scar punctures on leaves and necrotic patches on petioles) and caused flower abortion, at a similar level, however, only *N. tenuis* produced necrotic rings. Overall, *N. tenuis* females produced more damage to tomato fruit than *D. cerastii*. There was an increased frequency of *D. cerastii* females found on the plants over time, which did not happen with *N. tenuis*. Our results suggested that although *D. cerastii* caused less damage to tomato fruit than *N. tenuis*, it did feed on the fruit and could cause floral abortion, which requires field evaluation and caution in its use.

**Abstract:** Despite their importance as biological control agents, zoophytophagous dicyphine mirids can produce economically important damage. We evaluated the phytophagy and potential impact on tomato plants of *Dicyphus cerastii* and *Nesidiocoris tenuis*. We developed a study in three parts: (i) a semi-field trial to characterize the type of plant damage produced by these species on caged tomato plants; (ii) a laboratory experiment to assess the effect of fruit ripeness, mirid age, and prey availability on feeding injuries on fruit; and (iii) a laboratory assay to compare the position of both species on either fruit or plants, over time. Both species produced plant damage, however, although both species produced scar punctures on leaves and necrotic patches on petioles, only *N. tenuis* produced necrotic rings. Both species caused flower abortion at a similar level. Overall, *N. tenuis* females produced more damage to tomato fruit than *D. cerastii*. There was an increased frequency of *D. cerastii* females found on the plants over time, which did not happen with *N. tenuis*. Our results suggested that, although *D. cerastii* caused less damage to fruit than *N. tenuis*, it still fed on them and could cause floral abortion, which requires field evaluation and caution in its use in biological control strategies.

**Keywords:** omnivorous predator; *Nesidiocoris tenuis*; *Dicyphus cerastii*; plant damage; zoophytophagy; fruit injury; protected crops; biological control; tomato

#### **1. Introduction**

Zoophytophagous mirid species (Hemiptera: Miridae) are important biological control agents in several crops. Dicyphine (Miridae: Bryocorinae: Dicyphini) species, such

**Citation:** Souto, P.; Abraços-Duarte, G.; da Silva, E.B.; Figueiredo, E. Half Friend, Half Enemy? Comparative Phytophagy between Two Dicyphini Species (Hemiptera: Miridae). *Insects* **2022**, *13*, 175. https://doi.org/ 10.3390/insects13020175

Academic Editor: Nickolas G. Kavallieratos

Received: 31 December 2021 Accepted: 4 February 2022 Published: 7 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

as *Nesidiocoris tenuis* (Reuter), and several species of the genera *Macrolophus* Fieber and *Dicyphus* Fieber, are used worldwide as generalist predators [1] on several vegetable crops, both in conservation and augmentative biological control strategies.

*Dicyphus cerastii* Wagner is a Palearctic mirid, reported in the Mediterranean Basin, which spontaneously colonizes Portuguese greenhouses [2]. It is currently being evaluated as a candidate biological control agent (BCA), since it can feed on several horticultural pests [3]. *Nesidiocoris tenuis* is currently commercialized and released to control whiteflies and *Tuta absoluta* (Meyrick) in Mediterranean greenhouses [1,4,5]

Dicyphine mirids may resort to phytophagy in periods of prey scarcity [6,7], and to obtain water [8] and nutrients [9] from plants. Despite being advantageous as a feeding strategy, phytophagy can have negative effects in an agronomical context. Plant feeding may lead to a decrease in predation activity [9,10], and, more importantly, phytophagy can cause damage of economic importance, such as necrotic rings in stems and leaf petioles, as well as flower or fruit abortion, and punctures in the fruit [6,11–14]. This is particularly evident with *N. tenuis*, which is often the target of pesticide sprays that are used to control its populations when there is a risk of plant damage occurring, a practice that negatively impacts other natural enemies present on crops.

The increasing demand for food products without pesticide residues, combined with the need to control pests, highlights the urgency for sustainable alternatives that reduce negative effects on both the consumer and the environment's health. The use of predatory mirids in biological control has been very successful in protected tomato crops (e.g., [5,15,16]), despite the damage produced by some species. Therefore, to enhance biological control in tomato crops, research should focus on less phytophagous yet efficient dicyphine predators.

Plant damage can greatly vary with host plant species [6], mirid species [6,17,18], and even among populations of the same species [19,20]. Therefore, understanding the risk that phytophagy represents to crops is a key point in the evaluation of a new candidate BCA [13], such as *D. cerastii*, that will help decision making while choosing the best suited mirid species. Although damage caused by *N. tenuis* has been studied (e.g., [6,11,21,22]), little is known about those induced by *D. cerastii.* The latter has been observed producing chlorotic punctures on leaves [4] and also necrotic damage on tomato stems and leaf petioles, as well as feeding punctures on fruit (our pers. obs.). However, the influence of its damage on plant development and possible economic impacts has never been evaluated or compared with other mirid species.

We hypothesized that *D. cerastii* and *N. tenuis* display different phytophagous behavior. Within this context, the aim of this study was to compare the phytophagy and potential impact on tomato production of these two Dicyphini species. For this, we characterized the type of plant damage produced by these species on tomato plants. Then, we assessed the effect of fruit ripeness, mirid age, and prey availability on feeding injuries on tomato fruit. Finally, we compared the location of both species on either tomato fruit or plants.

#### **2. Materials and Methods**

#### *2.1. Rearing of Mirid Predators*

Stock colonies of both species (*D. cerastii* and *N. tenuis*) are maintained at the Instituto Superior de Agronomia (ISA), Lisbon, Portugal, on tobacco plants (*Nicotiana tabacum* L.). The colony of *D. cerastii* was started with individuals from different geographical sites in Portugal (Fataca, Ferreira do Zêzere, Lisbon, and Póvoa de Varzim) and is frequently refreshed with individuals, mostly from the Oeste region (Mafra and Silveira). The colony of *N. tenuis* was started with individuals from the Oeste region (Silveira) and from Koppert Biological Systems (The Netherlands). For rearing details, see [3]. Young adult females (between 1 and 7 days old), for all three bioassays, were obtained from the regular collection of large nymphs from breeding cages that were transferred to separate cages, where they could reach adulthood. For nymph experiments (see Section 2.3, Fruit Damage), 4th/5th instar nymphs were collected from immature rearing cages.

#### *2.2. Phytophagy in Semi-Field Conditions*

Phytophagy was observed for *N. tenuis* and *D. cerastii* in semi-field conditions in a greenhouse at ISA's *campus*. For this, mesh cages (1.5 m high and 1.0 m wide) were used. In each cage, there were two tomato plants (cv. Vayana), each one in a 15 L pot. Plants were fertilized and watered as needed, using an organic fertilizer solution (Húmus Líquido Horta®, SIRO, Mira, Portugal). When plants had 4 to 5 developed leaves (ca. 30–40 cm high), six couples of *D. cerastii* or *N. tenuis* were released. Control cages were set without any insects. Each treatment (12 *D. cerastii*, 12 *N. tenuis* and a control) had five replicates. To simulate the natural presence of prey, a teaspoon (ca. 4 g) of a mix of *Artemia* spp. (Anostraca: Artemiidae) and *Ephestia kuehniella* Zeller (Lepidoptera: Pyralidae) eggs (Entofood®, Koppert Biological Systems) was sprinkled evenly on the plants at the time of insect release, and every two weeks after that. Six weeks after release, all insects in the cages were captured into vials containing 96% ethanol and counted. Plant damage, expressed as the total number of necrotic rings or necrotic patches, was recorded for every cage. Flower abortion was counted as the proportion of missing flowers/total number of flowers in the flower rachis, for every cage. When trusses already had small fruit, these were counted as flowers. The mirids were released on 28 April 2021, and the assay ended on 11 June 2021.

#### *2.3. Fruit Damage*

The following factors were considered to compare the puncture level on fruit: mirid species (*D. cerastii* and *N. tenuis*), mirid age (nymph and adult female), tomato ripeness (unripe and fully ripe), and availability of food and/or water. For each species, three adult females or three nymphs were placed in plastic cups (8 cm high and 6 cm diameter) with one tomato fruit that was approximately 4–5 cm long. The lid of the cups had a hole (3 cm *φ*) covered with fabric to allow ventilation. Four treatments were considered (15 replications), for each species, mirid age, and tomato ripeness: (a) fruit only (N); (b) fruit with water (W) (supplied through an Eppendorf vial with moist cotton wool); (c) fruit with water and alternative food (FW) (*Ephestia kuehniella* eggs and *Artemia* spp. cysts Entofood®, Koppert Biological Systems, Berkel en Rodenrijs, The Netherlands), supplied on a sticky paper strip (2.0 cm <sup>×</sup> 0.8 cm); (d) fruit with alternative food (F) (Entofood®) but no water. Tomato fruit (Figure 1, cv. Suntasty) were obtained from an organic commercial tomato greenhouse. Before the experiment, the fruit were washed with abundant water and individually inspected for any marks or possible feeding punctures. Fruit with any defect were discarded. The unripe condition was considered as fully grown green fruit and fully ripe was considered as fully red fruit. For each bioassay, a control (i.e., only fruit) was made to ensure that punctures were caused by mirids. Insects were allowed to feed for 24 h, after which punctures (injuries resulting from fruit feeding) were counted under a stereomicroscope with a magnification of 50×. Injury was considered as a puncture surrounded by a small whitish or yellowish halo [12,23]. Replicates in which death occurred, or nymphs molted into adults, were discarded.

### *2.4. Location on Tomato Plant versus Fruit*

The tomato fruit were the same cultivar as those in the fruit-damage bioassay (see Section 2.3, above). Females of both species were placed individually in 600 mL transparent plastic cups covered with fabric to allow ventilation. The cups contained a young tomato plant (cv. San Pedro) held in water in a 10 mL glass bottle, and an unripe (green) tomato fruit. The unripe tomato was elected after analyzing the data from the tomato-damage bioassay, in which the most injured tomato was unripe. Females were placed individually for 24 h before the experiment in empty test tubes (starved). After release into the plastic cups, the position of the females was recorded at 1 h, 2 h, 6 h, and 24 h as being on the young plant, on the fruit, or elsewhere in the cup (as proposed by McGregor et al. [24]). Observations were conducted in a controlled chamber (Fitoclima CP500, Aralab Lda., Albarraque, Sintra, Portugal) at a temperature of 25 ◦C and photoperiod of 14 h. At the

end, fruit were inspected for feeding punctures. Females were used only once. In total, 20 replicates were made for each species. *Insects* **2022**, *12*, x FOR PEER REVIEW 4 of 15

**Figure 1.** Plant damage by *Nesidiocoris tenuis* on tomato plants: (**a**) necrotic rings on leaves and shoots; (**b**) shoot wilting; (**c**,**d**) detail of necrotic rings on leaves. **Figure 1.** Plant damage by *Nesidiocoris tenuis* on tomato plants: (**a**) necrotic rings on leaves and shoots; (**b**) shoot wilting; (**c**,**d**) detail of necrotic rings on leaves.

#### *2.4. Location on Tomato Plant versus Fruit 2.5. Data Analysis*

The tomato fruit were the same cultivar as those in the fruit-damage bioassay (see Section 2.3, above). Females of both species were placed individually in 600 mL transparent plastic cups covered with fabric to allow ventilation. The cups contained a young tomato plant (cv. San Pedro) held in water in a 10 mL glass bottle, and an unripe (green) tomato fruit. The unripe tomato was elected after analyzing the data from the tomato-Phytophagy in semi-field conditions. Both *D. cerastii* and *N. tenuis* populations at the end of the experiment and plant damage (necrotic rings or necrotic patches) were compared between species using one-way ANOVA with species as an independent variable. Differences in flower abortion among treatments were compared with Pearson's χ 2 tests. These statistical analyses were performed with IBM SPSS statistics v.26 (IBM, Armonk, NY, USA).

damage bioassay, in which the most injured tomato was unripe. Females were placed individually for 24 h before the experiment in empty test tubes (starved). After release into the plastic cups, the position of the females was recorded at 1 h, 2 h, 6 h, and 24 h as being on the young plant, on the fruit, or elsewhere in the cup (as proposed by McGregor et al. [24]). Observations were conducted in a controlled chamber (Fitoclima CP500, Aralab Lda., Albarraque, Sintra, Portugal) at a temperature of 25 °C and photoperiod of 14 h. At the end, fruit were inspected for feeding punctures. Females were used only once. In total, 20 replicates were made for each species. *2.5. Data Analysis*  Phytophagy in semi-field conditions*.* Both *D. cerastii* and *N. tenuis* populations at the end of the experiment and plant damage (necrotic rings or necrotic patches) were compared between species using one-way ANOVA with species as an independent variable. Differences in flower abortion among treatments were compared with Pearson's χ2 tests. These statistical analyses were performed with IBM SPSS statistics v.26 (IBM, Armonk, NY, USA). Fruit damage bioassay. Classification tree methods were used to understand the relative importance of the variables used (i.e., the absence or presence of water and/or food, developmental stage *of* the mirid, tomato ripeness, and species) on the number of feeding punctures. Statistical analyses were performed using R software version 4.1.0 implemented in RStudio version 1.4.1106. Conditional inference trees were made using the "*ctree*" function (R package party, http://cran.r-project.org/web/packages/party/index.html, accessed on 9 June 2021), which bases node splitting on statistical tests, providing a *p*-value for the significance of splitting [25]. The importance of the variables was measured using the random forest algorithm [26] and computations were performed in the randomForest package with 1001 trees (ntree = 1001). The random forest algorithm combines many classification trees to produce more accurate classifications and has measures of variable importance and measures of similarity of data points as by-products of its calculations [27]. Data were analysed together (i.e., considering counts for both species) and separately for each species. All preliminary analysis considered the modalities in two groups: (i) without food (N and W); and (ii) with food (F and FW). We grouped them into absence (A) and presence (P) of food and/or water, respectively, in the presented output.

Fruit damage bioassay. Classification tree methods were used to understand the relative importance of the variables used (i.e., the absence or presence of water and/or food,

Position on tomato plant versus fruit. The position of insects was compared between species for each observation time using the Fisher's exact test and z-test with Bonferroni correction method. For the comparison among locations within each species and each observation time, the non-parametric χ 2 test for one sample was used. These statistical analyses were performed with IBM SPSS statistics v.26 (IBM, Armonk, NY, USA).

#### **3. Results**

#### *3.1. Phytophagy in Semi-Field Conditions*

The average temperature during the assay was 24.3 ◦C, with a minimum of 10.4 ◦C and a maximum of 45.4 ◦C. The relative humidity was ca. 50%. At the end of the experimental period, *N. tenuis* had a larger average population (137.4 ± 30.3 individuals/cage) than *D. cerastii* (68.4 ± 14.4 individuals/cage); this difference, however, was not significant (F = 4.182; df = 1; *p* = 0.075).

*Nesidiocoris tenuis* produced both necrotic patches (ca. 10% of total damage) on leaves and stems and necrotic rings (ca. 90% of total damage) (Figure 1), whereas *D. cerastii* only produced necrotic patches (Figure 2). Necrotic rings caused by *N. tenuis* occasionally led to withering of young shoots or leaves (Figure 1b), while this was not observed with *D. cerastii* patches. Plant damage numbers were significantly different between species (F = 17.114; df = 1; *p* = 0.003), and *N. tenuis* produced more necrotic injuries on the plants than *D. cerastii* (Figure 3). *Insects* **2022**, *12*, x FOR PEER REVIEW 6 of 15

**Figure 2.** Plant damage by *Dicyphus cerastii* on tomato plants: (**a**) puncture scars on expanded leaves, (**b**–**d**) detail of necrotic patches on leaves. **Figure 2.** Plant damage by *Dicyphus cerastii* on tomato plants: (**a**) puncture scars on expanded leaves, (**b**–**d**) detail of necrotic patches on leaves.

Both species also produced puncture scars on leaves as a result of feeding on young stem/leaf tissues (Figure 2a).

Control cages had significantly lower flower abortion than both *D. cerastii* (χ <sup>2</sup> = 12.047; df = 1; *p* = 0.001) and *N. tenuis* cages (χ <sup>2</sup> = 16.395; df = 1; *p* < 0.001), whereas the two mirid species displayed similar levels of flower abortion (χ <sup>2</sup> = 0.670; df = 1; *p* = 0.413) (Figure 4).

**Figure 3.** Number of plant injuries (necrotic rings or patches) per cage, by *Dicyphus cerastii* and *Nesidiocoris tenuis* on tomato plants (two plants/cage). Bars topped by different letters represent

means that are significantly different (ANOVA, *p* < 0.05).

(**b**–**d**) detail of necrotic patches on leaves.

**Figure 3.** Number of plant injuries (necrotic rings or patches) per cage, by *Dicyphus cerastii* and *Nesidiocoris tenuis* on tomato plants (two plants/cage). Bars topped by different letters represent means that are significantly different (ANOVA, *p* < 0.05). **Figure 3.** Number of plant injuries (necrotic rings or patches) per cage, by *Dicyphus cerastii* and *Nesidiocoris tenuis* on tomato plants (two plants/cage). Bars topped by different letters represent means that are significantly different (ANOVA, *p* < 0.05).

**Figure 2.** Plant damage by *Dicyphus cerastii* on tomato plants: (**a**) puncture scars on expanded leaves,

**Figure 4.** Percentage (+SE) of flower abortion (missing flowers/total number of flowers × 100) on control, *Dicyphus cerastii*, and *Nesidiocoris tenuis* tomato plants. Bars topped by different letters represent significantly different percentages (χ2, *p* < 0.05). **Figure 4.** Percentage (+SE) of flower abortion (missing flowers/total number of flowers × 100) on control, *Dicyphus cerastii*, and *Nesidiocoris tenuis* tomato plants. Bars topped by different letters represent significantly different percentages (χ 2 , *p* < 0.05).

#### *3.2. Fruit Damage 3.2. Fruit Damage*

Both species fed on tomato fruit and produced punctures that appeared as damaged epicarp/mesocarp cells. It was often possible to observe that the damaged area extended beyond the puncture point following stylet movement inside the fruit. Punctures were structurally similar (Figure 5), but the pattern of each species was different. *Dicyphus cerastii* punctures tended to be aggregated, forming clearly visible patches in cases of highly damaged fruit. *Nesidiocoris tenuis* punctures appeared less aggregated compared to *D. cerastii*. Punctures on fruit did not heal, as punctures on green fruit did not disappear even when fruit changed color during maturation. Punctures produced by females and nymphs appeared similar, for both species. We observed that, occasionally, females of both species laid eggs on fruit. It is possible that an amount of the punctures may have been egg laying attempts or the result of probing. Both species fed on tomato fruit and produced punctures that appeared as damaged epicarp/mesocarp cells. It was often possible to observe that the damaged area extended beyond the puncture point following stylet movement inside the fruit. Punctures were structurally similar (Figure 5), but the pattern of each species was different. *Dicyphus cerastii* punctures tended to be aggregated, forming clearly visible patches in cases of highly damaged fruit. *Nesidiocoris tenuis* punctures appeared less aggregated compared to *D. cerastii*. Punctures on fruit did not heal, as punctures on green fruit did not disappear even when fruit changed color during maturation. Punctures produced by females and nymphs appeared similar, for both species. We observed that, occasionally, females of both species laid eggs on fruit. It is possible that an amount of the punctures may have been egg laying attempts or the result of probing.

A high variability was observed in the number of punctures inflicted in all treatments and in both species, with low numbers or even no punctures to high numbers of punctures in the same treatment, which generated great variance in the data for both species studied. When analyzing the dataset considering both species or for each species separately, there was no difference between the two treatments without food (N and W) or between the

the one with the highest number of feeding punctures, followed by the presence/absence of food. Food presence (F, FW) or absence (N, W) was only significant in the case of green fruit, and species was only significant for females in the presence of food in the case of green fruit, and for females in the ripe fruit. While the most important factor for *D. cerastii* was also the tomato ripeness, for *N. tenuis* the most important factor was the individual's stage of development (i.e., whether it is a nymph or an adult). Food was the second factor

for *D. cerastii* but only the third, and more distant, for *N. tenuis* (Figure 6 and 7).

**Figure 5.** Feeding punctures in tomato fruit: (**a**) unripe fruit; (**b**) ripe fruit. **Figure 5.** Feeding punctures in tomato fruit: (**a**) unripe fruit; (**b**) ripe fruit.

A high variability was observed in the number of punctures inflicted in all treatments and in both species, with low numbers or even no punctures to high numbers of punctures in the same treatment, which generated great variance in the data for both species studied. When analyzing the dataset considering both species or for each species separately, there was no difference between the two treatments without food (N and W) or between the two treatments with food (F and FW). Considering both species, the most important variable was the tomato ripening stage (tomato\_age), with the unripe (green) tomato being the one with the highest number of feeding punctures, followed by the presence/absence of food. Food presence (F, FW) or absence (N, W) was only significant in the case of green fruit, and species was only significant for females in the presence of food in the case of green fruit, and for females in the ripe fruit. While the most important factor for *D. cerastii* was also the tomato ripeness, for *N. tenuis* the most important factor was the individual's stage of development (i.e., whether it is a nymph or an adult). Food was the second factor for *D. cerastii* but only the third, and more distant, for *N. tenuis* (Figures 6 and 7).

#### *3.3. Location on Tomato Plant vs. Fruit*

The locations of *N. tenuis* and *D. cerastii* were only significantly different at the first observation (1 h) (Fisher's exact test value = 6.423, *p* = 0.033), with the former more present on the young plant than the latter and the inverse regarding the cup (α = 0.05) (Figure 8). Both *N. tenuis* and *D. cerastii* were mainly found on the young plant. However, in the case of *D. cerastii*, differences among locations were only verified at 2 h, 6 h, and 24 h, with females more present on the young plant (or on the young plant or on the cup walls) than on the fruit. In the case of *N. tenuis*, the females were more often found on the young plant than on fruit or cup walls, except at 24 h; at this time there were no differences between the young plant and fruit and no female was observed on cup walls.

In this bioassay, damage on plants was not quantified since, in some treatments, they would not be identified, especially with *D. cerastii*. Feeding punctures on fruit were found (although they were not counted) in all cases when the females of both species were sighted on the fruit. Furthermore, in both species, at least one fruit with feeding punctures was found, although the female was never observed in that same fruit during the bioassay (one case in *N. tenuis* and two cases in *D. cerastii*).

**Figure 6.** Variable importance plot from the random forest model (randomForest). The variables are ordered top-to-bottom as most-to-least important for an increase in feeding punctures (counts) on tomato fruit. (**a**) Using all datasets with data from both species, *Dicyphus cerastii* and *Nesidiocoris tenuis*; (**b**) database containing data collected only for *D. cerastii*; (**c**) database containing data collected only for *N. tenuis*. %incMSE: increase in mean square error of predictions as a result of the

variable being permuted; inNodePurity: importance of each predictor variable**.** 

**Figure 5.** Feeding punctures in tomato fruit: (**a**) unripe fruit; (**b**) ripe fruit.

**Figure 6.** Variable importance plot from the random forest model (randomForest). The variables are ordered top-to-bottom as most-to-least important for an increase in feeding punctures (counts) on tomato fruit. (**a**) Using all datasets with data from both species, *Dicyphus cerastii* and *Nesidiocoris tenuis*; (**b**) database containing data collected only for *D. cerastii*; (**c**) database containing data col-**Figure 6.** Variable importance plot from the random forest model (randomForest). The variables are ordered top-to-bottom as most-to-least important for an increase in feeding punctures (counts) on tomato fruit. (**a**) Using all datasets with data from both species, *Dicyphus cerastii* and *Nesidiocoris tenuis*; (**b**) database containing data collected only for *D. cerastii*; (**c**) database containing data collected only for *N. tenuis*. %incMSE: increase in mean square error of predictions as a result of the variable being permuted; inNodePurity: importance of each predictor variable.

variable being permuted; inNodePurity: importance of each predictor variable**.** 

lected only for *N. tenuis*. %incMSE: increase in mean square error of predictions as a result of the

**Figure 7.** Classification trees from the conditional inference trees (ctree) model. For each internal node, input variable and P values are provided, the boxplot of the number of feed punctures is displayed for each end node. Numbers in boxes above the variable indicate the node number. Number above boxes (n) indicates number of fruit. (**a**) Using all datasets with data from both species, *Dicyphus cerastii* and *Nesidiocoris tenuis*; (**b**) database containing data collected only for *D. cerastii*; (**c**) database containing data collected only for *N. tenuis*. Dc: *Dicyphus cerastii*; Nt: *Nesidiocoris tenuis*; A: absence of food; P: presence of food. *3.3. Location on Tomato Plant vs. Fruit*  **Figure 7.** Classification trees from the conditional inference trees (ctree) model. For each internal node, input variable and P values are provided, the boxplot of the number of feed punctures is displayed for each end node. Numbers in boxes above the variable indicate the node number. Number above boxes (n) indicates number of fruit. (**a**) Using all datasets with data from both species, *Dicyphus cerastii* and *Nesidiocoris tenuis*; (**b**) database containing data collected only for *D. cerastii*; (**c**) database containing data collected only for *N. tenuis*. Dc: *Dicyphus cerastii*; Nt: *Nesidiocoris tenuis*; A: absence of food; P: presence of food. than on fruit or cup walls, except at 24 h; at this time there were no differences between the young plant and fruit and no female was observed on cup walls. In this bioassay, damage on plants was not quantified since, in some treatments, they would not be identified, especially with *D. cerastii*. Feeding punctures on fruit were found (although they were not counted) in all cases when the females of both species were sighted on the fruit. Furthermore, in both species, at least one fruit with feeding punctures was found, although the female was never observed in that same fruit during the bioassay (one case in *N. tenuis* and two cases in *D. cerastii*).

The locations of *N. tenuis* and *D. cerastii* were only significantly different at the first

**Figure 8.** Number of *Dicyphus cerastii* (Dc) and *Nesidiocoris tenuis* (Nt) females on the young tomato plant (green), tomato fruit (red), or cup wall (grey) after 1 h, 2 h, 6 h, and 24 h. Bars topped by different letters for each observation time represent significant differences between species (Fisher's exact test, *p* < 0.05); different letters within the same column indicate differences among location for each species (χ2, *p* < 0.05). **4. Discussion Figure 8.** Number of *Dicyphus cerastii* (Dc) and *Nesidiocoris tenuis* (Nt) females on the young tomato plant (green), tomato fruit (red), or cup wall (grey) after 1 h, 2 h, 6 h, and 24 h. Bars topped by different letters for each observation time represent significant differences between species (Fisher's exact test, *p* < 0.05); different letters within the same column indicate differences among location for each species (χ 2 , *p* < 0.05).

plant and mirid species (e.g., [6]).

observed in *D. cerastii* plants.

The occurrence of plant damage production by zoophytophagous mirids is influenced by factors such as prey scarcity (e.g., [11,21,28]), water stress (e.g., [8]), and also host

phine mirids [6,17,18] since this species has the particularity to feed on vascular tissues and to aggregate on feeding sites [29]. Our study corroborates this, as plants with *N. tenuis* suffered more damage than those with *D. cerastii*. Moreover, the higher number of necrotic rings produced by *N. tenuis* were also more severe for the plant compared to the necrotic patches observed for *D. cerastii*. We observed that, in some *N. tenuis* infested plants, the apical shoots or leaflets withered because of necrotic rings in stems, whereas this was not

Besides the damage to vegetative parts of the plant, dicyphines are also reported to damage reproductive organs, such as flowers and fruit [6]. *N. tenuis* is also recognized for causing flower and fruit abortion on tomato plants [30], and in a study by Sanchez et al. [31] the percentage of flower abortion by *N. tenuis* reached up to 50% during population peaks. Flower and fruit abortion on tomato plants has also been reported for *M. pygmaeus* [32]. However, to our knowledge, this type of damage has not been previously described for *Dicyphus* spp. In our experimental conditions, we found that the percentage of flower abortion was not different between both mirid species. As flower abortion is particularly important on cluster tomato cultivars, *D. cerastii* may have a similar impact to *N. tenuis* on such cultivars, despite the lower damage to vegetative tissues produced by *D. cerastii*.

#### **4. Discussion**

The occurrence of plant damage production by zoophytophagous mirids is influenced by factors such as prey scarcity (e.g., [11,21,28]), water stress (e.g., [8]), and also host plant and mirid species (e.g., [6]).

Phytophagous behaviour in *N. tenuis* is regarded as more severe than in other dicyphine mirids [6,17,18] since this species has the particularity to feed on vascular tissues and to aggregate on feeding sites [29]. Our study corroborates this, as plants with *N. tenuis* suffered more damage than those with *D. cerastii*. Moreover, the higher number of necrotic rings produced by *N. tenuis* were also more severe for the plant compared to the necrotic patches observed for *D. cerastii*. We observed that, in some *N. tenuis* infested plants, the apical shoots or leaflets withered because of necrotic rings in stems, whereas this was not observed in *D. cerastii* plants.

Besides the damage to vegetative parts of the plant, dicyphines are also reported to damage reproductive organs, such as flowers and fruit [6]. *N. tenuis* is also recognized for causing flower and fruit abortion on tomato plants [30], and in a study by Sanchez et al. [31] the percentage of flower abortion by *N. tenuis* reached up to 50% during population peaks. Flower and fruit abortion on tomato plants has also been reported for *M. pygmaeus* [32]. However, to our knowledge, this type of damage has not been previously described for *Dicyphus* spp. In our experimental conditions, we found that the percentage of flower abortion was not different between both mirid species. As flower abortion is particularly important on cluster tomato cultivars, *D. cerastii* may have a similar impact to *N. tenuis* on such cultivars, despite the lower damage to vegetative tissues produced by *D. cerastii*.

In this study, the presence of water did not influence fruit puncture level. Therefore, both mirid species and both development stages could obtain the water they needed from green or ripe tomato fruit, at least when water was not provided. Water provision has been reported as one reason for phytophagy on heteropteran predators (e.g., [9]). However, as puncture numbers did not differ when water was supplied for both *N. tenuis* and *D. cerastii*, these mirids may have looked for other resources when they fed on the fruit.

Among the nutrients obtained from phytophagy, carbohydrates may have a particular ecological function since they have been reported to influence both predation and reproduction in dicyphines. This was demonstrated for *N. tenuis*, which was able to reduce the amount of prey feeding needed to establish itself on tomato plants [5], and increased its progeny [33], in the presence of sucrose dispensers. In another study, *N. tenuis* reduced its phytophagy when provisioned with sucrose dispensers [34].

In our study, when considering both species combined, tomato ripeness was the most important factor, with green fruit suffering more punctures than mature ones. This difference may be due to distinct nutritional profiles between unripe and ripe fruit. Sugar concentration, among other nutrients, may be higher in ripe tomato fruit [35,36], so it is possible that mirids may obtain more nutritional value per feeding puncture on ripe fruit than on green ones. On green fruit, the main effect was the presence of prey, which reduced fruit damage. There were differences between the species as *N. tenuis* females produced more damage than those of *D. cerastii*. On ripe fruit the most important factor was mirid age, with females producing most damage and, in these fruit, food did not significantly reduce damage. However, and once again, females of *N. tenuis* produced more damage than those of *D. cerastii*.

Considering *N. tenuis*, the most important factor on fruit damage was age, with females damaging more fruit than nymphs. Differently, in studies with whole plants, *N. tenuis* nymphs showed higher carbohydrate content [34] and spent more time feeding on the apical part of the plant [20], compared to adults. A similar trend was found for nymphs of the neotropical mirids, such as *Macrolophus basicornis* (Stål), *Engytatus varians* (Distant) and *Campyloneuropsis infumatus* (Carvalho), that also produced fruit punctures, whereas females did not [23]. Even though we observed less fruit damage by *N. tenuis* nymphs than females, our results indicated that nymphs were less influenced by factors such as fruit ripeness or presence of prey, suggesting that *N. tenuis* nymphs may be less prone to change

their phytophagous behavior than adults. Following mirid age, fruit ripeness was the next important factor for *N. tenuis*, with green fruit sustaining more damage. The presence of prey reduced the amount of damage on green fruit, whereas on ripe fruit it did not produce differences. This further suggests that green fruit may be a less valuable nutritional source for *N. tenuis*.

Tomato ripeness was the most relevant factor to explain fruit punctures by *D. cerastii*. The presence of prey was also important for fruit damage reduction in green fruit. Differently to *N. tenuis,* mirid age was not important in this species, which may suggest that *D. cerastii* may not be as dissimilar in phytophagy between adults and nymphs as *N. tenuis*.

Plant damage by zoophytophagous mirids has been associated with prey scarcity [11,21]. However, in our study the presence of food did not affect puncture level on ripe or green fruit with *N. tenuis* nymphs. Similarly, McGregor et al. [24] reported that the presence of food did not influence the level of fruit feeding by *Dicyphus hesperus* on mature tomato fruit, and Lucas and Alomar [37] reported that in whole caged plants the presence of *E. kuehniella* eggs did not prevent fruit injury by *D. tamaninii*.

Plant damage by zoophytophagous mirids may be determined by a complex combination of factors, besides prey abundance. Different species may have distinct preference or behavior that produce different types and levels of damage. Under the same conditions *M. caliginosus* (in fact, *M. pygmaeus*, C. Castañé, pers. comm.) did not produce fruit damage, whereas *D. tamaninii* did [37]. A different dicyphine, the nearctic *D. hesperus* preferred to feed on tomato leaves producing negligible damage on fruit [24]. Host plant and cultivar may also determine phytophagy, as was demonstrated for *N. tenuis*, which varied its phytophagy among different tomato cultivars [38]. The health of the host plant may also shape phytophagy by dicyphines. *Macrolophus pygmaeus* was reported to increase the number and produce more evident fruit damage on tomato plants infected with Pepino mosaic virus (PepMV) [12], but the same did not occur with *N. tenuis* [14]. Defence-activated plants may also be less susceptible to mirid phytophagous behavior. This was demonstrated for *N. tenuis*, which produced less plant damage on tomato plants inoculated with the endophytic *Fusarium solani* K strain, a fungal isolate that confers tomato resistance to foliar and root fungal pathogens [39].

Other factors may explain differences in phytophagy, such as genetic variation within species [19,20]. In fact, for the same treatments, we observed high variability in puncture numbers inflicted on fruit. As the large majority of the individuals used to initiate, and all the ones used to refresh the rearings, came from nearby locations (less than 45 km of linear distance), it is likely that the geographic origin was not a key determining factor in the high variability in feeding puncture number. This suggests that other factors, other than those considered in our study, may be driving fruit feeding in both *N. tenuis* and *D. cerastii*, and genetically determined behaviors should probably be considered in future research.

The fact that most *N. tenuis* and *D. cerastii* females were found on the young plants rather than elsewhere in the cup, and that there was an increased frequency of *D. cerastii* females found on young plant over time, may be related to the search for a better oviposition site [24]. Despite this, we could observe a slight increase over time of *N. tenuis* females occurring on fruit, which became the same as that for young plants at 24 h, suggesting a potential risk to fruit by this species. Furthermore, although few females of both species were observed on fruit compared to plant parts, feeding punctures were observed on fruit where females were not seen throughout the observations, showing that the female was at some moment on the fruit and fed on it. Finally, as *D. cerastii* preferred to be on young tomato plants than tomato fruit for plant feeding over time, the potential for damage to tomato fruit by this zoophytophagous mirid may be lower when compared to *N. tenuis.*

In the field, in commercial, protected tomato crops, necrotic rings, shoot, and flower cluster withering, and also punctures on fruit, are common when *N. tenuis* is present at high densities. This has repercussions on tomato production, leading growers to use a tolerance threshold and apply control measures. In the case of *D. cerastii,* necrotic tissues and punctures on fruit have been observed in the field by our team in commercial greenhouses when this species is present in high population densities. In order to fully assess how *D. cerastii* may affect tomato production (both in quantity and quality), further research is needed in semi-field conditions and commercial greenhouses, to establish safe population density thresholds. Since the damage caused by *D. cerastii* was apparently different from the necrotic rings of *N. tenuis,* histological studies are needed to characterize the necrotic patch damage reported here. Furthermore, it is important to understand if the feeding behaviour of *D. cerastii* induced the production of volatile defence compounds in the damaged plant, as reported for *N. tenuis*, *M. pygmaeus* and *Dicyphus maroccanus*, Wagner (syn. *D. bolivari* Lindberg) [40,41], with the consequent attraction of other biological control agents [41].

#### **5. Conclusions**

Overall *D. cerastii* damage was less severe than *N. tenuis*, as it did not cause necrotic rings and was more likely to seek out parts of the plant than the fruit. Despite this, it fed on fruit and caused flower abortion. Therefore, as was already known for *N. tenuis* and *M pygmaeus*, *D. cerastii* has the potential to cause an economic impact on tomato fruit production, particularly for cluster tomato cultivars, since its damage is related to the parts of the plant responsible for fruit production. We suggest that decision making regarding its use as a biological control agent should be made through field evaluation considering different cultivars.

We also found that fruit damage was highly variable within treatments, indicating that there may be differences in phytophagy on both species and individual levels. Therefore, in the future, selection of less phytophagous populations/strains combined with adequate management strategies may also benefit from the predatory behavior of dicyphine mirids with lower negative impact on tomato production.

**Author Contributions:** Conceptualization, P.S., G.A.-D., and E.F.; method, P.S., G.A.-D., and E.F.; validation, P.S., G.A.-D., E.B.d.S., and E.F.; formal analysis, P.S., G.A.-D., and E.F.; investigation, P.S., G.A.-D., and E.B.d.S.; resources, G.A.-D.; data curation, P.S. and G.A.-D.; writing—original draft preparation, P.S., G.A.-D., and E.F.; writing—review and editing, P.S., G.A.-D., E.B.d.S., and E.F.; visualization, P.S., G.A.-D., and E.F.; supervision, E.F.; project administration, E.F.; funding acquisition, E.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by national funds under the research project Umbert-ECO (PTDC/ASP-PLA/29110/2017), by the PhD grant to G.A.-D. (SFRH/BD/118834/2016), both funded by FCT-Fundação para a Ciência e Tecnologia, I.P., by the LEAF exploratory project TomBugBite, under FCT financial support for LEAF Research Centre (UIDP/04129/2020+UIDB/04129/2020) and by CEF Research Centre (FCT UIDB/00239/2020).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to give special thanks to Vera Zina (CEF-Forest Research Centre, Instituto Superior de Agronomia), Luiz Felipe Silveira (Western Carolina University, USA), and Filipe Madeira (MORE-Laboratório Colaborativo Montanhas de Investigação) for their support in statistical analysis; Horticilha Agro-Indústria S.A. for providing the fruit used in biological assays, the two anonymous reviewers for their constructive comments, and the editor for his review of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

