Adaptation of Drosophila subobscura chromosomal inversions to climatic variables: the Balkan natural population of Avala

The adaptive value of chromosomal inversions continues raising relevant questions in evolutionary biology. In many species of the Drosophila genus, different inversions have been recognized to be related to thermal adaptation, but it is necessary to determine to which specific climatic variables the inversions are adaptive. With this aim, the behavior of thermal adapted inversions of Drosophila subobscura regarding climatic variables was studied in the natural population of Avala (Serbia) during the 2014–2017 period. The results obtained were compared with those previously reported in the Font Groga (Barcelona, Spain) population, which presents different climatic and environmental conditions. In both populations, it was observed that most thermal adapted inversions were significantly associated with the first, second or both principal components, which were related with maximum, minimum and mean temperatures. Moreover, a significant increase over years (2004–2017) for the minimum temperature was detected. In parallel, a significant variation over time in Avala was only observed for the frequencies of ‘warm’ and ‘non-thermal’ adapted inversions of the U chromosome. However, stability in the chromosomal inversion polymorphism was observed for the 2014–2017 period which might result from the temporal span of the study and/or selective process acting on the population.


Introduction
Chromosomal inversions were detected one century ago by Sturtevant (1921), and Dobzhansky observed that their genetic content presented adaptive properties (for a revision see Lewontin et al. 1981). The presence of inversions is a karyotypic characteristic of a large number of species across the animal kingdom, and this structural chromosomal change is a key element in evolutionary processes such as adaptation to diverse environmental conditions or speciation (Sperlich and Pfriem 1986;Hoffmann et al. 2004;Feuk et al. 2005;Ayala and Coluzzi 2005;Hoffmann and Rieseberg 2008;McAllister et al. 2008;Kirkpatrick 2010;Nie et al. 2012;Fuller et al. 2018;Kapun and Flatt 2019). One of the best studied cases is the adaptive potential of chromosomal inversion polymorphism of Drosophila subobscura. For European researchers, this species is somehow the equivalent of D. pseudoobscura for North American evolutionary scientists. Nevertheless, D. subobscura has the advantage that presents inversions in all chromosomes of its karyotype with the exception of the dot (for a revision see Kunze-Mühl and Sperlich 1955;Krimbas 1992Krimbas , 1993. Interestingly in the Palearctic populations of this species, it was detected that the frequencies of inversions (or their combinations, the so-called arrangements) presented clines correlated with latitude (Krimbas and Loukas 1980;Krimbas 1992Krimbas , 1993. Although different explanations were possible, the discovery of latitudinal cline variation of chromosomal inversion frequencies in the North and South American colonizing 1 3 populations of D. subobscura which exhibit the same sign than that found in the Old World, strongly supported the adaptive role of inversions (Prevosti et al. 1988Ayala et al. 1989). To elucidate the adaptive meaning of D. subobscura inversions in Palearctic populations, Menozzi and Krimbas (1992) analyzed their possible relationship with many environmental variables. It must be remembered that D. subobscura, as a species, can only survive and reproduce under certain climatic conditions. This fact was clearly demonstrated by noting the distribution of D. subobscura in the Palearctic region and in those colonized areas of America (North and South). In these three zones, a parallel and gradual variation of climatic conditions (western maritime, Mediterranean and semiarid/arid) can be observed (The Times 1972;Prevosti et al. 1988Prevosti et al. , 1989. Even in a particular location and in the same season, the proportion over years of D. subobscura with regard to other species of the Drosophila genus depends on the weather conditions (Argemí et al. 1999(Argemí et al. , 2003Galludo et al. 2020).
Changes in the inversion chromosomal polymorphism over time in natural populations of this species regarding global warming have been extensively studied (Orengo and Prevosti 1996;Rodriguez-Trelles and Rodriguez 1998;Solé et al. 2002;Balanyà et al. 2004Balanyà et al. , 2006Balanyà et al. , 2009Rego et al. 2010;Dolgova et al. 2010;Zivanovic et al. 2012Zivanovic et al. , 2015Orengo et al. 2016). From early studies, inversions were classified as 'cold', 'warm' or 'non-thermal' adapted (Menozzi and Krimbas 1992;Rodriguez-Trelles et al. 1996;Rego et al. 2010;Zivanovic et al. 2016). However, the adaptive response of particular inversions in face of different climatic variables deserves to be better analyzed. Recently, a detailed study was carried out in the Font Groga (Barcelona) populations, where the relation of different inversions with regard to climatic variables (mean, maximum and minimum temperatures, humidity and rainfall) was analyzed (Galludo et al. 2018). The most relevant result of that research was that 'cold' inversions (A st , J st and E st ) increased in frequency when Tmax (maximum temperature) and Tmin (minimum temperature) decreased, whereas the 'warm' inversions (E 1+2+9 , O 3+4+1 and O 3+4+8 ), augmented in frequency with an increase in Tmax and Tmin. Furthermore, inversions seem to adapt not only to temperature, but also to humidity (Hm) and rainfall (Rf). For these reasons, it would be evolutionary relevant to compare the behavior of thermal adapted inversions in another different population of D. subobscura with regard to climatic conditions, season, different habitat and historic factors (another refugium of last glaciations).
The present research is articulated in two objectives regarding the thermal adapted chromosomal inversions of D. subobscura. The first objective is the main and aims to deepen the knowledge of the relationship between particular chromosomal inversions and some climatic variables. For this purpose, the chromosomal polymorphism of the Balkan D. subobscura population of Mount Avala (Serbia) was analyzed during four consecutive years. This population was previously studied (Zivanovic and Mestres 2010;Zivanovic et al. 2014aZivanovic et al. , 2015Zivanovic et al. , 2016 and could be very useful because it shows a different chromosomal polymorphism composition of that from Front Groga (Galludo et al. 2018;Zivanovic et al. 2016). Furthermore, both populations seem to differ regarding the climate and vegetation. Finally, the Iberian Peninsula and the Balkans are considered biodiversity refugia during the Pleistocene glaciations adding historical processes in their differentiation (Taberlet et al. 1998;Hewitt 2000;Tzedakis 2004;Provan and Bennett 2008;Wielstra and Arntzen 2020). Thus, Avala and Font Groga D. subobscura populations would be differentiated due to different historic and adaptive process. For these reasons, to compare both populations regarding the effects of climatic variables on the inversion frequencies could improve the knowledge of the role of chromosomal inversions on adaptation. To achieve this aim, our research was organized as follows: first, the possible differences in chromosomal inversions composition and abundance between Avala and Font Groga were analyzed. Also, a comparison between both populations with regard to some climatic variables was carried out. Third, a qualitative contrast between the vegetation of these two localities was made. Once detected the possible differentiation between both populations, a complete analysis to elucidate the effect of some climatic variables on different inversions in Avala population was computed. The results were compared with those from Font Groga. This comparison should allow to obtain more information on a relevant question: will the chromosomal inversions (those that are the same and those that are not) in different biological and physical conditions react in a similar adaptive way with respect to variations in some climatic variables?
The other objective is secondary, but related to first one, taking the advantage that Avala population was previously studied several years ago (Zivanovic and Mestres 2010;Zivanovic et al. 2015). It seeks to answer the following question: will thermally adapted inversions increase in frequency as our planet's global warming augments, or is there a maximum threshold for these frequencies? For this purpose, the chromosomal inversion data from the last year of the present research (2017) was compared with those available from Avala recorded in 2004 and 2011 (Zivanovic and Mestres 2010;Zivanovic et al. 2015) to analyze the changes over time in 'cold', 'warm' or 'non-thermal' adapted inversion frequencies. In parallel, the behavior of some climatic variables in Avala for the period 2004-2017 was analyzed to explore the potential global warming effect.

Drosophila subobscura collections and preparation of chromosomes
Flies were collected from Mount Avala population (Serbia, 44° 41′ 25″ N 20° 30′ 51″ E), located at 18 km south of Belgrade. They were sampled from a forest with polydominant communities of Fagetum submontanum mixtum that is approximately 450 m a.s.l. Collections were obtained strictly in the same place in 2014 (from the 10th to the 22nd of June), 2015 (from the 1th to the 14th of July), 2016 (from the 10th to the 18th of June) and 2017 (from the 9th to the 19th of June). In all samples, flies were netted from 4 to 8 p.m. using 40 fermenting apple baits.
Wild males and sons of wild females were crossed in individual vials with virgin females from the Kussnacht . To obtain the chromosomal preparations third instar larvae from the F 1 were dissected and polytene chromosomes were stained and squashed in aceto-orcein solution. To properly identify the inversions, the chromosomal maps of Kunze-Mühl and Müller (1958) and Krimbas (1992Krimbas ( , 1993 were used. To achieve the karyotypes with a probability higher than 0.99, at least eight larvae were analyzed from the progeny of each cross. All crosses were carried out at 18 °C, 60% relative humidity and 12 h/12 h light/dark cycle. Finally, chromosomal inversions were classified as 'C' (cold adapted), 'W' (warm adapted) and 'N' (non-thermal adapted) following the criterions of Menozzi and Krimbas (1992), Rego et al. (2010) and Zivanovic et al. (2016). For the aim of comparison, the Font Groga (Barcelona, 41° 25′ 54″ N 2° 07′ 20″ E) population is situated on the foothills of the Tibidabo mountain (on the edge of Barcelona city at 415 m a.s.l.). Flies were collected in a forest of pines (Pinus pinea) and ilexes (Quercus ilex) with Mediterranean brushwood, during the autumn of five consecutive years (2011-2015) (Galludo et al. 2018).

Climatic data and statistical analyses
The Serbian Republic Hydrometeorological Service provided the meteorological data from Avala (mean, maximum and minimum temperatures; dawn, noon, dusk and mean humidity; rainfall). Temperature, humidity and rainfall were measured in centigrade degrees, percentage and millimeters of precipitation, respectively.
All statistical computations were carried out using the basic and vegan packages of R language (R Development Core Team 2014). To compare the chromosomal composition between different years and departure of the observed frequencies of chromosomal karyotypes from expectations, Fisher's exact test was used (statistically significant P < 0.05). The corresponding P values were obtained using a bootstrap procedure (100,000 runs). In all cases of multiple comparisons, the FDR correction (Benjamini and Hochberg 1995) was applied, and it was reported as significant for P < 0.05. The index of free recombination (IFR) was computed according to Carson (1955). Using all D. subobscura chromosomes (A, E, J, U and O), a comparison between the annual samples from Avala (2014-2017) and Font Groga (2011Groga ( -2015 populations was carried out as follows. With this set of populations using the Bhattacharyya distance (Bhattacharyya 1946), a Principal Coordinate Analysis was performed according to Balanyà et al. (2006) and Mestres et al. (2009). Moreover, a cluster was computed with the same data using the GEVA-Ward procedure, because it is considered excellent for chromosomal inversion data (Irigoien et al. 2010;Zivanovic et al. 2016). To measure how faithfully the cluster preserved the pairwise distances between the original data, the Pearson cophenetic correlation was computed.
To compare the climatic variables (mean, maximum and minimum temperatures; dawn, noon, dusk and mean humidity; rainfall) between Avala and Font Groga during the same period of years (2014-2017) and using the values from May (the month before trapping the flies) a Mann-Whitney U test was computed. Furthermore, a Ward cluster was computed.
The relationship between chromosomal inversions and climatic variables (mean, maximum and minimum temperatures, mean humidity and rainfall) was computed only for those inversions observed in all Avala samples. For this reason, the infrequent inversions J 3+4 , U 1 , U 1+2+3 , O 7 , O 22 , O 3+4+2 , O 3+4+7 and O 3+4+17 were excluded of the analysis. Moreover, we gathered together the meteorological information for March, April and May during the studied period (2014-2017), because we assumed that the climate of these previous months before collecting the D. subobscura individuals could considerably influence the chromosomal composition of the samples. The mean of the values obtained for the climatic variables during these three months was used in the calculations (Supplementary Table S1). Using the climatic variables, a Principal Component Analysis (PCA) was carried out. In this way, the principal components synthesize the information provided by the climatic variables. For each chromosome and inversion, we considered a Poisson regression model taking as regressors these principal components. For the individual significance of each component, the FDR correction (Benjamini and Hochberg 1995) was applied and it was reported as significant for P < 0.05. Moreover, we assessed the ratio of variation in abundance for the chromosomal inversions when all climatic variables were kept constant and only one variable was either increased or decreased in only one unit.
Finally, in the Avala population the possible global warming effect was measured using data for mean, maximum and minimum temperatures, humidity and rainfall recorded in the period 2004-2017 for the average of months of March, April and May, because they were the months before the D. subobscura collections. For each climatic variable, a temporal series was computed. Differences between frequencies of 'cold', 'warm' and 'non-thermal' adapted inversions for the samples of 2004, 2011 and the period 2014-2017, were studied using the Fisher exact test, as previously explained. The thermal adaptation of the whole karyotype was measured using the CTI index .

Climate and vegetation comparisons between Avala and Font Groga
Regarding the kind of climate, both populations differ: Avala is classified as western maritime, whereas Font Groga is Mediterranean (The Times 1972; Prevosti et al. 1988). This translates in different kind of natural landscapes: Avala belonging to a forests and grasslands of temperate continental regions and Font Groga to a typical Mediterranean one. These features and historic factors during the Pleistocene conditioned the vegetation of both localities. In the case of Avala, it is composed of submontanum beech wood where Fagus moesiaca is the dominant tree together with three other linden species (Tilia argentea, Tilia cordata, Tilia platyphylos) and brushwood plants (Rubus hirtus, Pulmonaria officinalis, Asperula odorata, Carex silvatica, Lamium galeobdolon and others). In the Font Groga site, the vegetation is mainly composed of pines (Pinus pinea) and ilexes (Quercus ilex) with brushwood (Arbutus unedo, Ruscus aculeatus, Erica arborea, Hedera helix, Rubus ulmifolius, Smilax aspera, Laurus nobilis and others). Using the U Mann-Whitney test and applying the FDR correction, differences were not detected in the climatic variables when compared Avala and Font Groga populations for the period 2014-2017 (Tmax, P = 0.6650, adjusted P = 0.8817; Tmin, P = 0.7715, adjusted P = 0.8817; Tmean, P = 1.0, adjusted P = 1.0; Rainfall, P = 0.0304, adjusted P = 0.1215; Dawn Hum, P = 0.0294, adjusted P = 0.1215; Noon Hum, P = 0.0591, adjusted P = 0.1575; Dusk Hum, P = 0.1939, adjusted P = 0.3879; Mean Hum, P = 0.4705, adjusted P = 0.7528). However, considering the P values without FDR correction a certain trend of differentiation was observed for rainfall and two measures of humidity (dawn and noon). Finally, the cluster analysis generated three groups: one containing only the Avala sample of 2014, another with the remaining collections from Avala, and the last group with all Font Groga samples ( Supplementary Fig. S1). The cophenetic correlation coefficient was 0.954, indicating that the tree is a trustworthy reproduction from the distances used. For the studied climatic variables, the cluster result seems to support a certain differentiation between these two localities.

The chromosomal inversion polymorphism of Avala population
Focusing in the chromosomal polymorphism, the frequencies for each inversion for the years 2014-2017 in Avala population are presented in Table 1. It is worth noting that the J 3+4 inversion, generally associated with arid climates (Krimbas 1992(Krimbas , 1993, was detected in Avala for the first time in 2015. Likely, infrequent inversions U 1 , O 7 and O 22 are a product of very rare recombination events. The U 1+2+3 inversion was not reported in Avala until 2017 and its frequency was negligible. It is usually found in Southern Italy, Israel and Iran (Krimbas 1992). The O 3+4+2 and O 3+4+17 were observed for the first time in Avala in 2011, both at very low frequencies. The first one is found at low frequencies in different populations of the Palearctic region (Krimbas 1992;Galludo et al. 2018;Madrenas et al. 2020), but in high frequencies in American colonizing populations (Prevosti et al. 1988;Balanyà et al. 2003). The second is distributed in the Iberian Peninsula at very low frequency (Solé et al. 2002;Mestres et al. 1998;Galludo et al. 2018). The distribution of karyotypes in the samples of Avala population is presented in Table 2. There were no deviations from the Hardy-Weinberg equilibrium for any chromosome in any of the years analyzed (Supplementary Table S2). Regarding the IFR values, they are very similar in the 2014-2017 samples, indicating an equivalent amount of chromosomal polymorphism. The values estimated correspond to those of populations located in the central area of D. subobscura distribution (Krimbas 1992).

Comparison between the chromosomal polymorphisms of Avala and Font Groga
Qualitative differences can be observed between the chromosomal polymorphisms of Avala and Font Groga regarding the type and abundance of inversions. The multivariate analyses (Principal Coordinate Analysis and cluster) demonstrated the existence of a differentiation (Fig. 1). In the PCoA, first, second and third axes explained 67.75%, 20.52% and 5.44%, respectively. The samples from both populations are clearly separated in the graphical representation (Fig. 1a). A similar result can be observed in the cluster analysis (Fig. 1b), where two groups are defined, one for the Avala samples and the other for those from the Font Groga. The cophenetic correlation coefficient was 0.902, which is interpreted as that the tree accurately describes the genetic distances between the samples used.

Relations between inversions and climatic variables
The first three principal components (PC1, PC2 and PC3) obtained using the climatic variables account for 57.52%, 38.49% and 3.99% of the variance, respectively. Thus, using the three components we conserve the 100% of the explained variability and no information is lost. PC1 is mainly related with temperature (Tmean, Tmax and Tmin) and a certain

Changes over time of the Avala chromosomal polymorphism
The climatic variables studied during 2004-2017 presented annual variations ( Supplementary Fig. 2). Although all temperatures increased over time, only Tmin was significant (P = 0.037, with a slope of 0.092). Large fluctuations were observed for humidity with a trend to decrease over time, although it was not significant. Finally, the rainfall showed the expected erratic distribution.  All comparisons between different years and chromosomes (A, J, U, E and O), but considering 'cold', 'warm' and 'non-thermal' adapted inversions are presented in Supplementary Table S3. It is worth to remember that there are not 'non-thermal' adapted inversions for the A chromosome. Significant values were observed only for the 'warm' and 'non-thermal' adapted inversions of the U chromosome when comparing the years 2004 with 2011 and 2004 with 2017. The variations of 'cold', 'warm' and 'non-thermal' adapted inversions over time are shown in Fig. 2. The 'cold' inversions presented a certain pattern of variation in different directions depending on the chromosome during the interval 2014-2016, to return to levels similar to those of 2004 (Fig. 2a). Regarding the 'warm' inversions, their behavior is such as that of the 'cold' inversions, with the clear exception of those from the U chromosome, which frequency was increasing from 2004 to 2011 and then seemed to suffer only small fluctuations (Fig. 2b). Finally, the changes of the 'non-thermal' inversion over time were insignificant, with the exception of those from the U chromosome showing the opposite behavior of the 'warm' inversions for the same chromosome (Fig. 2c). In summary, it seemed that the    increase of U chromosome 'warm' inversions would be at the expense of the 'non-thermal' for the same chromosome.
The CTI values for Avala in 2014, 2015, 2016 and 2017 were 0.383, 0.262, 0.302 and 0.371, respectively. In general, these values are rather similar, but a little bit lower, that those reported in the same population in 2004 (0.374) and 2011 (0.426) ). All possible statistical comparisons between these CTI values were computed and no evidence of significance was observed (Supplementary  Table S4).

Discussion
There is still a gap between the information provided by the molecular studies of inversions and their phenotypes regarding the adaptation to environmental conditions. The in situ hybridization and genomics studies have been able to identify many structural genes located inside chromosomal inversions of the Drosophila genus (for instance, Clark et al. 2007;Laayouni et al. 2007;Guillén and Ruiz 2012;Fabian et al. 2012;Orengo et al. 2017;Reis et al. 2018;Karageorgiou et al. 2019;. Furthermore, molecular studies have allowed understanding the recombination and linkage disequilibrium for genes located inside and outside of inversions (for example, Schaeffer et al. 2003;Laayouni et al. 2003;Schaeffer and Anderson 2005;Pegueroles et al. 2010Pegueroles et al. , 2013Pegueroles et al. , 2016Stevison et al. 2011;Smukowski et al. 2015). However, we have not completely identified all genetic elements (structural genes or regulatory sequences) located inside inversions or outside but close to them. Neither and even more difficult, we have not properly understood their possible interactions at functional level. Due to all these shortcomings in our knowledge and also from an evolutionary point of view, it is still fundamental to study the chromosomal inversions as whole genetic units although the understanding of their molecular content is still incomplete (Pennisi 2017).
Although Font Groga and Avala are populations located in the central area of the D. subobscura distribution (Krimbas 1992(Krimbas , 1993, they differ regarding historic, climatic and biotic factors. For this reason, their inversion chromosomal polymorphisms are distinct in term of composition and frequencies. This species presents a high capacity of dispersion (Begon 1976;Serra et al. 1987;Ayala et al. 1989) and gene flow (Araúz et al. 2009;Mestres et al. 2009;Pegueroles et al. 2013) that would tend to unify the chromosomal composition between populations. Thus, given the particular environmental conditions of each population, natural selection would likely have a great effect on adaptive inversions generating the differentiation. However, exactly the same thermal adapted inversions ('cold adapted': A st , A 1 , J st , U st , E st and O st ; 'warm adapted': A 2 , J 1 , U 1+2 , U 1+8+2 , E 1+2+9 , E 1+2+9+12 , O 3+4 , O 3+4+1 and O 3+4+8 ) were present in both populations, Font Groga (Galludo et al. 2018) and Avala. In the latter, most of them (9 out of 15, after the FDR correction) were significant with PC1, PC2 or both, being these principal components related mainly with temperature. These results differed from that from Font Groga (Galludo et al. 2018), where 14 out of 15 inversion were significant, but it is worth noting that the principal components in that population were not equivalent to those from Avala. For instance, the PC1 from Font Groga depended strongly of temperatures, but also of humidity and rainfall. Furthermore, in Avala population, the relation of each inversion with the variation of one unit for a particular environmental variable while maintaining constant the others, indicated that most inversions (with the exceptions of the E 1+2+9+12 and O 3+4+6 inversions) do not change appreciably. This result could be a consequence of annual fluctuations of climatic variables in Avala during the period 2014-2017, as no clear directional change was present for any of the studied variables, no changes were detected for the thermal adapted inversions.
However, it was relevant to study the change of chromosomal inversions over time in Avala population (2004Avala population ( -2017. From a qualitative point of view, several inversions considered warm and arid were detected recently for the first time and in low frequencies in Avala (J 3+4 in 2015 and U 1+2+3 2017). It could be explained by a tendency towards decreased humidity in that region, and those inversions could have arrived from the Eastern Mediterranean region. Instead, the quantitative study of inversion chromosomal frequencies produced the main results. Thus, for all polymorphic chromosomes of D. subobscura, the 'cold' adapted inversions remained rather constant for the period 2004-2017. In a global warming scenario, a decrease in these frequencies would be expected. A quite similar result was observed for the 'warm' inversions with a notable exception, the U chromosome, showed a significant increase from 2004 to 2011 samples mainly due to the increase in frequency of U 1+2 inversion (Zivanovic et al. 2015). This inversion was most abundant in the Caucasus and Iran regions (Krimbas 1993). However, from 2011 no significant variation in frequency has been observed for 'warm' inversions. The paradox was that the increase of U 'warm' inversions was at the expense of 'non-thermal' adapted inversion, not at the 'cold' adapted as it could be expected. For these reasons, the 'non-thermal' adapted inversions remained quite stable over time with the exception of the U chromosome, with a remarkable decrease in the U 1+2+6 frequency from 2004 to 2011.
It would seem that the chromosomal inversion frequencies of D. subobscura from Avala have reached a plateau of stability. This result could be an artifact due to the few years studied, and because inversion polymorphism data were not available for each one of the 13 years the involved in the analysis. Perhaps longer periods of time (about 1 3 20-30 years) would be necessary to detect notable changes in inversion frequencies, as was carried out in several previous researches with this species (Orengo and Prevosti 1996;Solé et al. 2002;Balanyà et al. 2004Balanyà et al. , 2006Orengo et al. 2016;Zivanovic et al. 2019). However, in intermediate periods of time (about 15 years), changes in this polymorphism were detected and in the direction expected by global warming (Rodriguez-Trelles and Rodriguez 1998; Zivanovic and Mestres 2011;Zivanovic et al. 2014b). Other hypotheses are possible, for instance the occurrence of episodes of intense cold or heat waves in any of the years in which the study is framed (Schar and Jendritzky 2004;Founda and Giannakopoulos 2009;Barriopedro et al. 2011;Stefanon et al. 2012). It was observed that D. subobscura inversion polymorphism was able to respond to these changes (Rodríguez-Trelles et al. 2013). However, these climatic episodes could generate particular results that could be wrongly interpreted from the perspective of climate change. This distortion due to particular climatic conditions could also be mitigated if the studied period is large enough.
On the topic of inversion frequencies variation according to global warming, there are still relevant questions to answer: is there a limit to the frequency increase of warm adapted inversions in a population? If there is a threshold, are D. subobscura populations close to it? Previous researches would indicate that the accumulation of 'warm' adapted inversions is additive, with interactions representing a secondary role, if any (Zivanovic et al 2016;Arenas et al. 2018). The thermal adaptation QTL mapping in D. melanogaster would support this observation (Morgan and Mackay 2006;Norry et al. 2008;Borda et al. 2018). Usually, all characters on which directional selection acts have a limit, because other fitness-related characters are adversely affected (for a revision see Lynch and Walsh 1998; Sgrò and Hoffmann 2004). For D. subobscura, the 'warm' adapted inversions could contain genes or have genes in linkage disequilibrium located near them, which would be dragged by the directional thermal selection generating serious alterations of fitness. Moreover, it is worth to remember the complex organization of chromosomal inversions in this species. From a historic point of view and due to its large number of heterozygotes for inversions collected in natural populations it was considered that D. subobscura could carry systems of balanced lethals (Krimbas 1992(Krimbas , 1993 for a complete revision). Although inversions with heterotic effect carrying lethal genes have been described in this species (Mestres et al. 2001), it was generally observed that the frequency of homozygotes and heterozygotes for inversions fits panmictic expectations. Moreover, some arrangements constituted by non-overlapped inversions located far away in the same chromosome were not broken by recombination, likely due to linkage disequilibrium (Sperlich and Feuerbach 1969;Sperlich and Feuerbach-Mravlag 1974;Mestres et al. 1998). If one of the inversions in the arrangement contained the thermal adapted genes, all inversions constituting this arrangement would be selected together. This selection would have consequences over other fitness traits, and for this reason trade-offs would likely be necessary. Fitness trade-offs involving inversions have been recently described in D. melanogaster (Durmaz et al. 2018).
Is there any evidence indicating that any D. subobscura population has already reached a 'warm' adapted inversions accumulation threshold? In a recent study in the D. subobscura isolated Atlantic population of Madeira , the frequencies of thermal adapted inversions and the CTI remained quite stable over time, although there was a clear indication of global warming in the island (significant increase of Tmin, Tmax and Tmean). The frequencies of the 'cold' adapted inversions A st (0.228) and E st (0.257) were not negligible, indicating that they could contain genes adaptive to other prominent environmental conditions for the species survival or the inversion system of the species does not allow to continue with the directional selection. Avala is not an isolated population and could change also due to gene flow, receiving climatic adapted inversions. In this sense, it will be interesting to follow up the evolution in frequency of J 3+4 , U 1+ 2+3 and O 3+4+2 inversions, considered adaptive to warm and dry conditions. Moreover, the CTI values estimated in this study (0.330 in average) are not particularly high, and in this sense the composition of 'warm' adapted inversions could increase, because the maximum value (0.958) was from the Italian population of Etna . Probably, the stability detected in Avala could be a transient situation and longer periods of time should be studied in the future. Furthermore, the inversion composition of Avala must be adaptive, not only regarding the thermal variables, but also to other environmental conditions. In conclusion, it would seem that depending on the environmental conditions, each D. subobscura population could accumulate until a certain maximum amount of thermal adapted inversions without altering other fitness components. Thus, the threshold for thermal adapted inversions accumulation would probably be different in each population depending on the particular environmental conditions of its ecosystem.
Finally, it is worth noting that natural ecosystems are complex in structure. The climatic variables as temperatures, humidity and rainfall act changing the communities and D. subobscura is an element of them. To study the direct effect of these variables on the chromosomal inversion polymorphism is possibly a reduction and simplification of the biological problem. Probably, inversions contain genes that could confer adaptive advantage in front a combination of conditions, both physical and biological present in the ecosystems (Galludo et al. 2018;Kapun and Flatt 2019). The interactions of genes located in the chromosomal structure of D. subobscura karyotype would be responsible for adaptations to a wide range of situations that take place in the ecosystems. All these topics deserve future studies because the adaptive value of inversions to different climatic conditions has been described in many species Ananina et al. 2004;Etges and Levitan 2008;Takahashi and Takano-Shimizu 2011;van Heerwaarden et al. 2012;Ayala et al. 2014;Berg et al. 2017;Kapun and Flatt 2019).