Genetic Variability, Population Differentiation, and Correlations for Thermal Tolerance Indices in the Minute Wasp, Trichogramma cacoeciae

Simple Summary Augmentative biological control relies on the more or less frequent/abundant releases of biological control agents (BCAs) that have to be adapted to their short-term local environment including (micro-)climatic conditions. Thermal biology of BCAs is thus a key component for their success. The extent to which thermal tolerance indices may be relevant predictors of the field efficiency is however still poorly documented. Within this frame, we investigated the intraspecific variability for the ability to move at low temperatures in the minute wasp, Trichogramma cacoeciae. We collected, molecularly characterized, and compared for their thermal tolerance indices numerous strains originating from three contrasting geographic areas. Our findings evidenced both a geographic differentiation between strains for one of the thermal tolerance indices and a positive correlation between two of them, demonstrating the existence of an intraspecific variability. Abstract Temperature is a main driver of the ecology and evolution of ectotherms. In particular, the ability to move at sub-lethal low temperatures can be described through three thermal tolerance indices—critical thermal minimum (CTmin), chill coma temperature (CCT), and activity recovery (AR). Although these indices have proven relevant for inter-specific comparisons, little is known about their intraspecific variability as well as possible genetic correlations between them. We thus investigated these two topics (intraspecific variability and genetic correlations between thermal tolerance indices) using the minute wasp, Trichogramma cacoeciae. Strains from T. cacoeciae were sampled across three geographic regions in France—two bioclimatic zones along a sharp altitudinal cline in a Mediterranean context (meso-Mediterranean at low elevations and supra-Mediterranean at higher elevations) and a more northwestern area characterized by continental or mountainous climates. Our results evidenced a significant effect of both the longitude and the severity of the cold during winter months on CCT. Results were however counter-intuitive since the strains from the two bioclimatic zones characterized by more severe winters (northwestern area and supra-Mediterranean) exhibited opposite patterns. In addition, a strong positive correlation was observed between CCT and CTmin. Neither strain differentiation nor the covariations between traits seem to be linked with the molecular diversity observed on the part of the mitochondrial marker COI.


Introduction
The process of local adaptation is a central theme in evolutionary ecology and, among selective abiotic factors, temperature is thought to be of major importance, in particular for ectotherms [1,2]. Within this frame, environmental gradients represent ideal spatial organizations to investigate population differentiation related to thermal adaptation [3,4]. In particular, both latitudinal and altitudinal gradients usually offer a progressive variation of the mean ambient temperature along variable geographic distances and, thus, experimentally tractable case-studies [5,6]. Yet, some drawbacks-e.g., co-variation between temperature and other abiotic factors, and distinction between direct and indirect causes-have to be considered cautiously.
Low temperatures are a significant source of stress for many insect species and are supposed to limit fitness and ultimately the distribution of species [7,8]. One particular aspect of thermal adaption is the ability of individuals to maintain a locomotor activity at low temperatures, which can be described using three thermal indices: (i) the critical thermal minimum (CTmin), i.e., the temperature at which the individuals lose locomotor activity, (ii) the chill coma temperature (CCT), i.e., the temperature at which the individuals lose all movement, and (iii) the activity recovery (AR), i.e., the temperature at which animals recover locomotor activity. Several authors emphasize the relevance of these thermal tolerance indices to experimentally investigate inter-and, to a lesser extent, intraspecific variabilities [9][10][11][12][13][14][15]. At the intra-specific level, the extent to which these indices actually reflect local adaptations and how they covary remain however poorly documented apart from some "model organisms" such as Drosophila species.
Besides academic purposes, eco-evolutionary processes affecting thermal biology are also of prime interest for more applied purposes, in particular in the frame of biological control-the use of living organisms (called biological control agents or BCAs) against noxious ones [16]. Incidentally, it is noteworthy that the "ins and outs" differ according to the biological control strategy. Because classical biological control aims at the longterm establishment of an exotic BCA, the adequacy of this latter towards the thermal environment in the future area of introduction is considered as a key parameter of its establishment success. Climate matching approaches are in such a case actually relevant (see, for instance, [17][18][19]) and some intra-specific genetic variations may indeed explain contrasted results [20]. However, because augmentative biological control (in particular the inundative sub-strategy) only aims at a transitory BCA presence and pest control, the thermal adequacy of BCAs has to be reflected and then optimized with regard to some particular (micro)climatic conditions. More precisely, those conditions may (i) only represent part of the seasonal variations (open field crops), (ii) be disconnected from the actual climate and artificially controlled (greenhouses productions), (iii) be modulated by the architecture and physiology of the plants (see, for instance, [21][22][23]), and/or (iv) be far from the BCA's "comfort zone". With regard to this last point, a relevant challenge is thus to choose or even genetically improve biological control agents that could be used at low temperatures early in the season in order to prevent pest outbreaks. To our knowledge, such attempts remain rare [24][25][26]. For macro-organisms used as biocontrol agents (e.g., predatory arthropods or insect parasitoids), one of the limits to this challenge is probably the absence of relevant and experimentally tractable proxies for evaluating the overall performances at low temperatures. We hypothesized that the thermal tolerance indices (critical thermal minimum, chill coma temperature, and activity recovery) may play the role of such proxies.
Given this context (both fundamental and applied research), we investigated the intra-specific differentiation for thermal tolerance indices as well as the correlation between these traits in the minute parasitic wasp, Trichogramma cacoeciae Marchal (Hymenoptera, Trichogrammatidae). As with other Trichogramma species, T. cacoeciae is an oophagous parasitoid whose biology and ecology can be summarized as follows [27]: The female wasp lays each of its eggs inside a host egg whose embryo is rapidly killed. The offspring's pre-imaginal development then occurs within the host egg until metamorphosis and the emergence as an adult. In T. cacoeciae, reproduction is ensured through thelytoky (diploid females producing diploid daughters without fecundation by males); this parthenogenesis not being induced by heritable symbionts [28]. Because asexuality is traditionally considered as a desirable feature for biological control [29], T. cacoeciae has often been evaluated or even used as a biocontrol agent (see, for instance, [30][31][32]). Moreover, T. cacoeciae is known to be a widespread species and has been sampled repeatedly in Europe [33]. At a very local scale (south-eastern France), T. cacoeciae is also the most abundant species along an altitudinal gradient [34]. However, the relative contributions of recurrent recolonizations through active and/or passive dispersal and local adaptations to contrasted climates remain poorly investigated. It is worth noting that, at this scale, no use of T. cacoeciae for augmentation biological control is documented so that natural eco-evolutionary processes are not disturbed by external inputs of individuals. To test for the presence of local adaptation, we constituted a pool of T. cacoeciae strains collected at both macro-(latitudinal cline) and micro-(altitudinal gradient) scales. The strains were characterized genetically using part of the cytochrome oxidase I [35], and compared phenotypically using the three thermal tolerance indices: CTmin, CCT, and AR. We also investigated correlations between these three thermal tolerance indices in order to point out possible negative and positive correlations between traits that are, respectively, indicative of evolutionary constraints or strategies.

Terminology
For this experiment, 40 strains of T. cacoeciae originating from south-eastern France were used (Table 1 and Figure 1). Here, the term "strain" refers to a group of wild-caught females (hereafter referred as to the "founders") collected in the same location, in the same micro-habitats (same plant and hosts' patch), at the same time and from which a perennial laboratory offspring has been obtained. These strains were representative of three geographic and climatic contexts, hereafter referred as to the "origins". Three origins were distinguished here: "southern meso-Mediterranean", "southern supra-Mediterranean" and "northwestern" (see details below). emergence as an adult. In T. cacoeciae, reproduction is ensured through thelytoky (diploid females producing diploid daughters without fecundation by males); this parthenogenesis not being induced by heritable symbionts [28]. Because asexuality is traditionally considered as a desirable feature for biological control [29], T. cacoeciae has often been evaluated or even used as a biocontrol agent (see, for instance, [30][31][32]). Moreover, T. cacoeciae is known to be a widespread species and has been sampled repeatedly in Europe [33]. At a very local scale (south-eastern France), T. cacoeciae is also the most abundant species along an altitudinal gradient [34]. However, the relative contributions of recurrent recolonizations through active and/or passive dispersal and local adaptations to contrasted climates remain poorly investigated. It is worth noting that, at this scale, no use of T. cacoeciae for augmentation biological control is documented so that natural eco-evolutionary processes are not disturbed by external inputs of individuals. To test for the presence of local adaptation, we constituted a pool of T. cacoeciae strains collected at both macro-(latitudinal cline) and micro-(altitudinal gradient) scales. The strains were characterized genetically using part of the cytochrome oxidase I [35], and compared phenotypically using the three thermal tolerance indices: CTmin, CCT, and AR. We also investigated correlations between these three thermal tolerance indices in order to point out possible negative and positive correlations between traits that are, respectively, indicative of evolutionary constraints or strategies.

Terminology
For this experiment, 40 strains of T. cacoeciae originating from south-eastern France were used (Table 1 and Figure 1). Here, the term "strain" refers to a group of wild-caught females (hereafter referred as to the "founders") collected in the same location, in the same micro-habitats (same plant and hosts' patch), at the same time and from which a perennial laboratory offspring has been obtained. These strains were representative of three geographic and climatic contexts, hereafter referred as to the "origins". Three origins were distinguished here: "southern meso-Mediterranean", "southern supra-Mediterranean" and "northwestern" (see details below).

Field Collection and Laboratory Rearing
Thirty strains were initiated from founders (about ten/strain, all these individuals being stored after their reproduction for the molecular characterization, see Section 2.3 § below) collected between 2016 and 2017 in a narrow geographic area (latitudes between 43.56 • N and 43.80 • N, longitudes between 6.84 • E and 7.19 • E) in south-eastern France (Department of "Alpes-Maritimes", hereafter "southern"). This area is characterized by a marked altitudinal cline from the sea level to an altitude of 1500 m in less than 30 km. Based on climatic and ecological features [36], this altitudinal cline encompasses two bioclimatic zones, the meso-Mediterranean and the supra-Mediterranean, respectively represented here by 14 southern-meso strains and 16 southern-supra strains (see Table 1). The diversity and abundance of Trichogramma within this area is detailed elsewhere [34] but it is noteworthy that:

•
Although there was an intense sampling effort [37] locations, and more than 9000 patches of the sentinel eggs of Ephestia kuehniella Zeller (Lepidoptera: Pyradidae) exposed], Trichogramma individuals were difficult to collect. On average, only 2% of hosts' patches were actually parasitized.

•
Among the eight recovered Trichogramma species, T. cacoeciae was the most abundant one and T. evanescens was the only species present all along the altitudinal cline.
The last 10 strains (hereafter, "Northwestern strains") were sampled within a more northwestern area (latitudes between 43.79 • N and 48.46 • N, longitudes between 2.49 • E and 6.63 • E) characterized by more continental or even mountainous climates. All the strains of this last set were collected using sentinel eggs of E. kuehniella between 2015 and 2016 except the TCMZ strain which is an old reference strain (see Table 1). Because Northwestern strains were collected for the needs of a different project, less ecological information is available than for the southern ones. Moreover, the founders were not kept for molecular characterization (see § below), which was performed on the lab-reared offspring.
After sampling, all strains were then reared in standardized conditions within the Biological Resource Center "Egg Parasitoids Collection" (doi.org/10.15454/AY4LMT) using the substitution host E. kuehniella. Abiotic conditions were kept constant with a mean temperature of 18 • C ± 1 • C, a relative humidity of 75% ± 10%, and a photoperiod (light: dark) of 16 h:8 h (see for instance [26]). When the phenotyping for the thermal indices was initiated, most of the strains (main exception: TCMZ) had undergone between 15-30 generations in laboratory depending on their sampling date (see Table 1).

Molecular Characterization
This molecular characterization aims at discarding the existence of cryptic species that are frequent in Hymenopteran species [38][39][40][41] as well as the hypothesis of a phylogeographic process [37,42,43]-for instance, divergences at the molecular and phenotypic levels in an allopatric context followed by, after the expansion of some lineages, a secondary contact-that would better explain the observed patterns in thermal biology than processes of local adaptations. To achieve these goals, the mitochondrial gene cytochrome c oxidase subunit 1 (COI) was used for similar work (see [35,38]). It is worth noting that, T. cacoeciae being thelytokous, there is no genetic mixing including mito-nuclear introgression.
The molecular characterization was performed on founders (for the southern mesoand southern supra-Mediterranean strains) or on lab-reared offspring (three individuals) for Northwestern strains. Genomic DNA was extracted from recently killed individuals' tissues using the prepGEM ® Insect kit (ZYGEM, PIN0500). Individuals were placed individually in 15 µL of mix and incubated for 3 h at 75 • C and then for 5 min at 95 • C. DNA extracts were stored at −20 • C.
A fragment (about 700 bp) was amplified using the primers LCO1490 and HCO2198 [44] and 1 µL of DNA was used for the PCR reaction (performed in a total volume of 25 µL). The sequences obtained were then aligned and compared with a set of previously obtained COI haplotypes see Supplementary Materials Table S1 and [45]. A neighborjoining tree was inferred using the Kimura 2-parameters distance and 500 replicates for bootstraps. All the analyses were performed using the MEGA software.

Thermal Tolerance Indices Phenotypic Characterization
As Trichogramma pre-imaginal life is parasitic and occurs within the host, phenotyping was performed on the mobile adults (less than 1 mm).
A prototype was especially designed and built for this experiment. This device allowed us to control the temperature of the air within a small circular arena (3 × 1 cm) where a small group of Trichogramma had been enclosed. Briefly, the cooling and warming were provided by a Peltier module controlled by a software developed for this purpose. The temperature of the arena was monitored in real time while the behavior of the individuals was observed from the top through a transparent wall and video-recorded using a digital video camera (Nikon D800 digital camera equipped with a 60 mm 1:2.8 G ED AF-S Micro Nikkor lens), the light being provided from the bottom.
After preliminary tests, each phenotyping session was done as follows: First, a group of about 50 young adults (i.e., emerged within the last 24 h) were transferred without anesthesia in the arena. The stage of "young adult" was used as only adults are mobile and as age may alter some thermal indices (David et al., 1998). Second, the temperature within the arena was initially set at +17 • C, a temperature at which the wasps were kept for no more than five minutes. Then, a progressive decrease of the temperature was controlled at a rate of 0.3 • C/ min until, in general, −2 • C ("routine program" on Table 1) and, sometimes, −4 • C ("extended program" on Table 1) was applied. Once the minimum was reached, an increase of the temperature occurred at the same rate until +14 • C. Each T. cacoeciae strain was phenotyped once. All the video-sequences were then carefully checked and only the behavior of the individuals (about 15-20 individuals) located on the floor of the arena and well visible were taken into account to estimate: • Critical thermal minimum (CTmin): temperature at which the last individual lost its ability to walk; • Chill coma temperature (CCT): temperature at which the first individual toppled with no more movement; • Activity recovery (AR): temperature at which the first individual recovered its ability to walk.

Statistical Analyses
Statistical analyses were performed with the software R (version R-3.6.0) and its working environment Rstudio (Version 1.1.453).
The intraspecific differentiation was tested on each of the three thermal tolerance indices (CTmin, CCT, and AR) using two complementary approaches-an hypothesis test approach and a model comparison approach realized with linear mixed effects models (LME models). Two sets of explanatory variables were investigated separately:

•
The first set of variables (hereafter, the hypothesis test approach) included the "origin" (qualitative variable with three modalities: southern-meso, southern-supra, and northwestern) as the only fixed effect. The sampling location and the COI haplotype were concatenated into a single qualitative variable and used as a random effect. This was to prevent any pseudo-replication linked to the use of strains deriving from the same "natural population" (see Field Collection and Laboratory Rearing).

•
The second set of variables (hereafter, the model comparison approach) considered the two climatic variables (Tmean and Tmini-see previous section) as well as three geographic ones (altitude, latitude, and longitude) as fixed effects. These quantitative variables thus substituted the qualitative variable "origin" used in the first analysis.
As for the first analysis, the concatenated information about the sampling location and the COI haplotype was used as a first random effect (intercept). Insofar as several strains are linked to the same meteorological station, this information was used as a second random effect (intercept).
In both cases (the hypothesis test approach and the model comparison approach), the analysis started with the selection of random effects based on in the first case on the Akaike information criteria (AIC) and in the second case on the conditional Akaike information criteria (cAIC) [46,47].
For each trait, all the models, i.e., the null one, the model with "origin" as the sole fixed factor in hypothesis test approach, and the 21 models obtained with the model comparison approach-were compared (models leading to singularity being discarded) and those minimizing the information criterion were conserved. If fixed effects were involved, they were tested using likelihood ratio tests (LRT) or Fisher's exact test (see Table 3). For each thermal index (respectively, CTmin, CCT, and AR) and each approach (hypothesis test and, when relevant, model comparisons-see Materials and Methods), the model selection was realized using the Information Criteria (see Supplementary Materials Table S2). The selected fixed factors were tested using likelihood ratio test (LRT) and/or Fisher's exact test (f -value).
For significant qualitative variables, post-hoc Tukey tests were also performed. Several R packages were used (in particular cAIC4, multcomp, effects, and ggplot) for the statistical analysis or graphic representations.
In order to explore possible correlations between thermal tolerance indices, pairwise Spearman's nonparametric rank correlation tests were performed to test possible correlations between the three thermal tolerance indices. Then, a principal component analysis (PCA) was carried out after the normalization of the data in order to graphically visualize the between-trait correlations as well as individual and mean performances (R packages: FactoMineR and Factoextra).

Molecular Characterization Using COI
The 40 strains were characterized using part of the COI. As shown in Figure 2, six valid haplotypes were obtained from the 40 strains, three ("Hap 006", "Hap 019" and "Hap 119") of them representing 94% of the strains and the other ones being observed only once.

Behavior of the Wasps at Cold Temperatures
When the temperature decreased, all strains lost, at some point, their ability to move their limbs in a coordinated manner and were therefore unable to walk. At colder temper atures, individuals started to topple on one of their sides with no more movement. Within a pool of individuals, there was, however, a large overlap between individual CTmin and CCT what explains the following results. Indeed, the median value for the last individua reaching CTmin was 0 °C with an interquartile between −1 °C and 0 °C. The median value for the first individual reaching CCT was +1°C with an interquartile between 0 °C and +1 °C. Then, when temperature increased, the mobility was restored. The median value for the AR was +9 °C with an interquartile of +8.50 °C and +9 °C. AR values have to be treated with caution insofar as, the CCT varying between strains, the strains had variable expo sures to cold temperatures before the warming. No mortality was observed during the

Behavior of the Wasps at Cold Temperatures
When the temperature decreased, all strains lost, at some point, their ability to move their limbs in a coordinated manner and were therefore unable to walk. At colder temperatures, individuals started to topple on one of their sides with no more movement. Within a pool of individuals, there was, however, a large overlap between individual CTmin and CCT what explains the following results. Indeed, the median value for the last individual reaching CTmin was 0 • C with an interquartile between −1 • C and 0 • C. The median value for the first individual reaching CCT was +1 • C with an interquartile between 0 • C and +1 • C. Then, when temperature increased, the mobility was restored. The median value for the AR was +9 • C with an interquartile of +8.50 • C and +9 • C. AR values have to be treated with caution insofar as, the CCT varying between strains, the strains had variable exposures to cold temperatures before the warming. No mortality was observed during the phenotyping. Figure 3 more precisely represents the distributions of each of the three thermal tolerance indices according to the three geographical origins (northwestern, southern meso-Mediterranean, and outhern supra-Mediterranean).

Genetic Differentiation between Geographic Origins
The hypothesis test approach ("origin" as the sole explanatory variable) did not allow detecting relevant factors for CTmin and activity recovery (Table 3 and Figure 3). By contrast, evidence for spatial differentiation was observed for the chill coma temperature, the fixed effect "origin" being significant. More precisely, a significant difference was observed between the "Northwestern" origin (median of −1.4 • C) and the "southern supra-Mediterranean" one (median of −0.4 • C) (z = 3.256; p = 0.003). No significant difference was observed for the two other comparisons, "northwestern" versus "southern meso-Mediterranean" (z = 2.090, p = 0.073) or "southern meso-Mediterranean" versus "southern supra-Mediterranean" (z = 1.260, p = 0.208). meso-Mediterranean, and outhern supra-Mediterranean).

Genetic Differentiation between Geographic Origins
The hypothesis test approach ("origin" as the sole explanatory variable) did not allow detecting relevant factors for CTmin and activity recovery (Table 3 and Figure 3). By contrast, evidence for spatial differentiation was observed for the chill coma temperature, the fixed effect "origin" being significant. More precisely, a significant difference was observed between the "Northwestern" origin (median of −1.4 °C) and the "southern supra-Mediterranean" one (median of −0.4°C) (z = 3.256; p = 0.003). No significant difference was observed for the two other comparisons, "northwestern" versus "southern meso-Mediterranean" (z = 2.090, p = 0.073) or "southern meso-Mediterranean" versus "southern supra-Mediterranean" (z = 1.260, p = 0.208). For each thermal tolerance index, the boxplot indicates the median, the interquartile, the minimums and maximums as well as any values considered extreme. For each index, the significance of the variable "origin" is provided in Table 3 (Hypothesis test approach). Results of post hoc tests are represented by different letters (a and b to indicate a significant difference and ns to indicate a non-significant difference). Boxplot in red, blue and green color respectively represent the three geographic origins: southern meso-Mediterranean strains (SM), southern supra-Mediterranean (SS), northwestern strains (N).
The model comparison approach gave different outcomes at the step of model selection according to the thermal index considered (see Supplementary Materials S2). For CTmin (linear mixed effects model), the values of cAIC were very close to each other (from 104.08 to 104.30) and all the models are equivalent to each other in information and to the null model. For the activity recovery, the values of AIC were close (from 172.3 to 181.2) and the best model was the null one, i.e., those without fixed effects. For the chill coma temperature (linear mixed effects model), the values of CAIC were far more variable (from 179.13 to 226.02), the best model being those involving three fixed effects: the daily mean For each thermal tolerance index, the boxplot indicates the median, the interquartile, the minimums and maximums as well as any values considered extreme. For each index, the significance of the variable "origin" is provided in Table 3 (Hypothesis test approach). Results of post hoc tests are represented by different letters (a and b to indicate a significant difference and ns to indicate a non-significant difference). Boxplot in red, blue and green color respectively represent the three geographic origins: southern meso-Mediterranean strains (SM), southern supra-Mediterranean (SS), northwestern strains (N).
The model comparison approach gave different outcomes at the step of model selection according to the thermal index considered (see Supplementary Materials Table S2). For CTmin (linear mixed effects model), the values of cAIC were very close to each other (from 104.08 to 104.30) and all the models are equivalent to each other in information and to the null model. For the activity recovery, the values of AIC were close (from 172.3 to 181.2) and the best model was the null one, i.e., those without fixed effects. For the chill coma temperature (linear mixed effects model), the values of CAIC were far more variable (from 179.13 to 226.02), the best model being those involving three fixed effects: the daily mean temperature observed during the three coldest months (Tmean), the daily minimal temperature observed during the same period (Tmini), and the longitude (cf. Table 3). Two of these (Tmean and longitude) were proven significant, the chill coma temperature being negatively correlated with the mean temperature during the coldest months (Tmean) and positively with the longitude (Table 3 and Figure 4). This thus confirms the results of the hypothesis test approach. temperature observed during the three coldest months (Tmean), the daily minimal temperature observed during the same period (Tmini), and the longitude (cf. Table 3). Two of these (Tmean and longitude) were proven significant, the chill coma temperature being negatively correlated with the mean temperature during the coldest months (Tmean) and positively with the longitude (Table 3 and Figure 4). This thus confirms the results of the hypothesis test approach. Figure 4. Influence of the climatic variable "Tmean" and the geographic variable "Longitude" on the chill coma temperature. The points, linear trend, and gray ribbon, respectively, represent the raw data and the predicted regression as well as the 95% confidence interval. The statistical significance of these two variables is indicated in the Table 3.

Pairwise and Multivariate Correlations between Thermal Tolerance Indices
Among the three pairwise correlations between CTmin, CCT, and AR, only those between CTmin and CCT were significant and highly positive (Spearman's rank correlation tests: rho = 0.43, p-value = 0.006). On the contrary, no correlation was found between AR and CCT (rho = −0.14, p-value = 0.38) or between AR and CTmin (rho = −0.11, p-value = 0.52).
To further investigate possible covariations between these three traits, a principal component analysis (PCA) was performed on normalized data, the first two axes explaining 83% of the variability ( Figure 5). Most of the variability on the first axis was explained by the positive correlation between CTmin and CCT while most of the variability on the second axis was explained by the AR. Within this phenotypic space, both the meso-and supra-Mediterranean strains appeared to be more variable than the Northwestern ones. In particular, 29% of the meso-Mediterranean strains and 44% of supra-Mediterranean ones exhibited high values on the first PCA axis (see dotted ellipse on Figure 5A). Such high values were not encountered in the northwestern strains. As shown in Figure 5B, differences in phenotypic performances did not depend on the COI haplotypes. Figure 4. Influence of the climatic variable "Tmean" and the geographic variable "Longitude" on the chill coma temperature. The points, linear trend, and gray ribbon, respectively, represent the raw data and the predicted regression as well as the 95% confidence interval. The statistical significance of these two variables is indicated in the Table 3.

Pairwise and Multivariate Correlations between Thermal Tolerance Indices
Among the three pairwise correlations between CTmin, CCT, and AR, only those between CTmin and CCT were significant and highly positive (Spearman's rank correlation tests: rho = 0.43, p-value = 0.006). On the contrary, no correlation was found between AR and CCT (rho = −0.14, p-value = 0.38) or between AR and CTmin (rho = −0.11, p-value = 0.52).
To further investigate possible covariations between these three traits, a principal component analysis (PCA) was performed on normalized data, the first two axes explaining 83% of the variability ( Figure 5). Most of the variability on the first axis was explained by the positive correlation between CTmin and CCT while most of the variability on the second axis was explained by the AR. Within this phenotypic space, both the meso-and supra-Mediterranean strains appeared to be more variable than the Northwestern ones. In particular, 29% of the meso-Mediterranean strains and 44% of supra-Mediterranean ones exhibited high values on the first PCA axis (see dotted ellipse on Figure 5A). Such high values were not encountered in the northwestern strains. As shown in Figure 5B, differences in phenotypic performances did not depend on the COI haplotypes.

Discussion
Our study evidences two main findings, the existence of a significant geographic differentiation for the chill coma (CC) and the positive correlation between the chill coma (CC) and critical thermal minimum (CTmin). Because the strains were reared during several generations in standardized conditions before their phenotyping, these two patterns are unlikely to be explained by (trans-generational) plastic responses. They thus rather reflect an actual intra-specific genetic variability. It is also noteworthy that the two patterns cannot be linked with the molecular clustering observed on COI. Indeed, based on surveys in other areas [45], the COI haplotype "Hap 006" (the most frequent one in cluster 1) incidentally appears widely distributed in Europe and also present in North Africa which suggests an ancient expansion from its original area. With such a wide geographic

Discussion
Our study evidences two main findings, the existence of a significant geographic differentiation for the chill coma (CC) and the positive correlation between the chill coma (CC) and critical thermal minimum (CTmin). Because the strains were reared during several generations in standardized conditions before their phenotyping, these two patterns are unlikely to be explained by (trans-generational) plastic responses. They thus rather reflect an actual intra-specific genetic variability. It is also noteworthy that the two patterns cannot be linked with the molecular clustering observed on COI. Indeed, based on surveys in other areas [45], the COI haplotype "Hap 006" (the most frequent one in cluster 1) incidentally appears widely distributed in Europe and also present in North Africa which suggests an ancient expansion from its original area. With such a wide geographic distribution, it is thus not surprising that contrasting life-history strategies may have been selected according to contrasting local conditions. So far, the COI haplotype "Hap 019" appears to be less cosmopolitan and, to date, only observed in the south-east of France. This does not however prevent this haplotype from including a wide range of phenotypes. Although the use remains, of course, of limited relevance to investigate population genomics, we can, however, discard the hypothesis that the phenotypic differentiation between bioclimatic zones is explained by different ancient lineages. We would thus rather explain the two patterns (significant differentiation on CC and positive correlation between CC and CTmin) as the consequence of local adaptions. As highlighted in Supplementary Materials Table S2, this study thus actually contributes to filling the deficit of knowledge about the intraspecific variability of thermal tolerance indices in insects and, in particular, in parasitoids.
Under temperate, mountainous, or even polar climates, the survival of ectotherms to low temperatures is a central question in understanding their ecology and evolution [48]. According to (i) the intensity and duration of the cold, (ii) the species' biology, and (iii) the stage of the development, various behavioral and/or physiological responses are observed including the research of local favorable micro-habitats, migration towards more favorable locations, "dormancy" (diapause or quiescence), and/or an improved tolerance [49,50].
With respect to this last point, several authors emphasize the relevance of the three thermal tolerance indices, critical thermal minimum [9], chill coma temperature, and activity recovery  as relevant proxies of the adaptation to low temperatures. As proxies, the mean values per geographic origins (meso-Mediterranean, supra-Mediterranean, and northwestern) observed for each thermal tolerance index are not ecologically meaningful per se but the ranking between the three origins (e.g., for CC: northwestern < meso-Mediterranean < supra-Mediterranean) is informative. This same reasoning was already held by Andersen and colleagues [9] and Le Laan and colleagues [52] at an interspecific level.
At the intraspecific level, several studies evidenced some strain differentiation according to latitude [53] and/or altitude. Nevertheless, the pervasiveness of such patterns and the scale at which they occur remain unclear. Our understanding of such processes is also limited by the lack of information about the potential genetic correlations between these thermal tolerance indices as well as between these indices and other facets of the adaptation to cold temperatures (e.g., lowest temperature for the development or propensity to diapause) although this aspect is mandatory to understanding possible constraints or synergies between traits [53]. In fact, the few references we have found about correlations between thermal tolerance indices deal with interspecific comparisons [9,10] or with environment-induced phenotypic variations [54]. The extrapolation of such results at the intra-specific level remains unreliable.
Within this frame, we thus addressed the question of the spatial scale at which strain differentiation and local adaptations may occur by considering both a macro-("southern" versus "northwestern" strains) and a micro-(within southern strain, "meso-" versus "supra-Mediterranean" strains) scales. Taken as whole, our results evidenced a phenotypic differentiation at the macro-scale with a significant effect of both the longitude and the severity the cold during winter. This led to a significantly lower median for the CCT in northwestern strains than in the southern supra-Mediterranean ones. Because environmental conditions were standardized during the phenotyping and, before, during several generations of lab-rearing, this differentiation can only be explained by genetic differences. This result was quite counter-intuitive insofar as we expected similar trends for these two origins for which similar cold climates are observed during winter. The reasons why the micro-process (genetic differentiation along the altitude cline) differ from the macro-process (genetic differentiation along the latitudinal/longitudinal one) are not understood. Discrepancies between the effect of, on one side, altitude, and, on the other, latitude, were however already observed for other traits and species [55][56][57][58]. Moreover, the strong correlation between the CCT and the CTmin suggests either independent but convergent selections on these two traits, or a genetic (pleiotropy) or physiological linkage [59,60].
A future direction of this study is to complete the phenotyping of these strains using other traits involved in the adaptation to cold temperatures. Based on the literature on Trichogramma species, two complementary facets are thus planned to be investigatedoverwintering strategies in field conditions [61][62][63], and propensity to diapause in laboratory conditions [64,65]. The objective will, hence, be to investigate if global strategies for low temperature exist and, if so, what could be the possible trade-off with other parts of the thermal range, the "comfort zone" or sub-lethal warm temperatures [46,66,67]. In addition to their relevance for fundamental research, such results might open major perspectives for the genetic improvement of Trichogramma strains used as biological control agents.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/insects12111013/s1 (the link will be generated): Table S1: Summary of reference strains used for the molecular characterization of the Trichogramma cacoeciae used in this study. Haplotypes corresponding to the strains used in this study are indicated in bold. The columns "reproduction" and "amino-acid sequence" indicate, respectively, the mode of reproduction of the strain (A = arrhenotoky, T = thelytoky) and the variation in the corresponding amino-acid sequence (label "varprot") and thus putative pseudogene or technical artefacts. Upper part: strains available at the Biological Resource Center "Egg Parasitoids Collection" (doi.org/10.15454/AY4LMT) (source: Warot 2018). Identification was formerly realized by B. Pintureau. Lower part: Genbank references. Table S2: Summary of the statistical models evaluated in the study. According to the selection of the random effects, the nature of the model was either a linear mixed effects model (LME) or a linear model (LM). For each of the traits, the two approaches are presented, i.e., the test hypothesis ("Origin" as the sole fixed effect) and the model comparison. In this second case, the null model and the 31 possible combinations obtained from the two meteorological variables (t_mean and t_mini) and the geographic ones (alt, lat, and lon) are provided. The model minimizing the Information criteria (cAIC for LME and AIC for LM) is indicated in bold.