Comparative analysis of adaptive and neutral markers of Drosophila mediopunctata populations dispersed among forest fragments

Abstract Comparison of adaptive and neutral genetic markers is a valuable approach to characterize the evolutionary consequences of populations living in environments threatened by anthropogenic disturbances, such as forest fragmentation. Shifts in allele frequencies, low genetic variability, and a small effective population size can be considered clear signs of forest fragmentation effects (due to genetic drift) over natural populations, while adaptive responses correlate with environmental variables. Brazilian Atlantic Forest had its landscape drastically reduced and fragmented. Now, several forest remnants are isolated from each other by urban and crop areas. We sampled Drosophila mediopunctata populations from eight forest remnants dispersed on two adjacent geomorphological regions, which are physiognomic and climatically quite distinct. Microsatellite data of inversion‐free chromosomes (neutral genetic marker) indicate low structuration among populations suggesting that they were panmictic and greatly influenced by gene flow. Moreover, significant differences in chromosomal inversion frequencies (adaptive genetic marker) among populations and their correlations with climatic and geographical variables indicate that genetic divergence among populations could be an adaptive response to their environment. Nonetheless, we observed a significant difference in inversion frequencies of a population in two consecutive years that may be associated with edge and demographic effects. Also, it may be reflecting seasonal changes of inversion frequencies influenced by great temperature variation due to edge effects. Moreover, the forest fragment size does not affect genetic variation of neutral markers. Our data indicate that despite oscillations in chromosomal inversion frequencies, D. mediopunctata populations from Brazilian Atlantic Forest and their divergence may be driven by adaptive factors to local differences, perhaps because it is a small flying insect easily carried by the wind increasing its migration rates.


| INTRODUC TI ON
Changes of natural landscape in different ecosystems around the world caused by forest fragmentation processes may be very harmful to biodiversity (Millette & Keyghobadi, 2015;Newmark & McNeally, 2018) and may affect the genetic structure of populations (Radespiel & Bruford, 2014;Rhoads, Williams, & Krane, 2017;Rosche et al., 2018). An appropriate approach to infer the nature of genetic responses in stressful environments is to compare population parameters such as number of migrants per generation, fixation index, and allele frequencies shifts using adaptive and non-adaptive genetic markers (Hoffmann & Willi, 2008;Merilä & Hendry, 2014;Stojanova et al., 2018).
Microsatellite loci are codominant multi-allelic genetic markers, which allow assessing both temporal and spatial genetic structure of natural populations (Gredler, Hish, & Noor, 2015;Hartvig et al., 2018;Silva, Machado, & Mateus, 2015). They are also a marker commonly used in conservation genetics to estimate the loss of genetic variability and to infer the demographic history of populations, assuming they are neutral or nearly neutral even if located in a coding region (Ellegren, 2004;Lombaert et al., 2018;Stamenković-Radak et al., 2012;Takezaki, 2017).
Forest fragmentation and deforestation are believed to make environmental conditions more heterogeneous, with pronounced changes in biotic and abiotic conditions (Keyghobadi, 2007).
Stochastic shifts in frequencies of genetic markers and high divergence among populations are expected to be the major responses of populations from a fragmented landscape (Milligan et al., 2018;Schippers et al., 2015).
Currently, the Atlantic Rainforest biome in Brazil is extremely fragmented (Joly, Metzger, & Tabarelli, 2014). About 80% of its forest remnants encompass areas smaller than 50 ha (Ribeiro, Metzger, Martensen, Ponzoni, & Hirota, 2009). Many of its forests remnants are scattered amid pastures, agricultural fields and growing urban landscape, especially in the states of São Paulo, Rio de Janeiro and Minas Gerais (Joly et al., 2014). There are more than a hundred forest fragment remnants in the city of Campinas, which is 61 km north of the Tropic of Capricorn (Cielo-Filho & Martins, 2016;Cielo-Filho, Gneri, & Martins, 2007). Although all forest remnants from this region can be classified as Seasonal Semi-deciduous Forests, var-ious studies have shown that they are heterogeneous (Cielo-Filho & Martins, 2016;Cielo-Filho et al., 2007;Salis, Shepherd, & Joly, 1995). Furthermore, they are located over an ecotone in the transition area between two geomorphological regions (GMRs): Peripheral Depression and Atlantic Plateau (Joly et al., 2014;Ross, 2013).
Marked differences are observed between these two GMRs.
Considering climatic variables such as average monthly temperature and annual sum of monthly precipitation, they differ in approximately 1°C and more than 70 mm, respectively (Table 1; see also Alvares, Stape, Sentelhas, Moraes Gonçalves, & Sparovek, 2013). In relation to their geological properties, Peripheral Depression unit is characterized by a flat topography and lands of magmatic sedimentary rocks with crystalline rocks; and Atlantic Plateau is characterized by orogenic belts, a continuous range of mountains with deep valleys and channels, with several soil types including cambisols, lithic, podzolic and podzolic yellow-red and red-yellow oxisol, and rocky outcrops (Ross, 2013).
Chromosomes II, IV, and X are polymorphic for inversions (with 17, two and four gene arrangements, respectively) while chromosomes III and V are inversion-free (Ananina et al., 2002;Brianti, Ananina, & Klaczko, 2013). Thus, to avoid biases on the genetic structure estimates by hitchhiking effect due to selection on chromosome inversions (Kennington & Hoffmann, 2013;Santos et al., 2016), one may use microsatellite loci located in the inversion-free chromosomes (Cavasini, Batista, & Klaczko, 2015;Laborda, Gazaffi, Garcia, & Souza, 2012). On the other hand, adaptive responses can be inferred from chromosome II inversion polymorphism, which is associated with variation in other traits that affect fitness, such as size and shape of the wing and genitalia and polychromatism (Andrade, Vieira, Ananina, & Klaczko, 2009;Bitner-Mathé, Peixoto, & Klaczko, 1995;Hatadani & Klaczko, 2008;Hatadani, Baptista, Souza, & Klaczko, 2004). Furthermore, chromosome II inversion polymorphism shows correlation with climatic variables (temperature and precipitation) congruent with an altitudinal cline and seasonal cycling variation described for the natural population from Parque Nacional do Itatiaia, Rio de Janeiro State (Ananina et al., 2004;Batista et al., 2012;Klaczko, 2006). Ananina et al. (2004) studying various populations of D. mediopunctata in a geographic transect found a puzzling result. In spite of their geographic distance, two natural populations showed stark differences in chromosome II inversion frequencies. Santa Genebra population, located in Peripheral Depression, showed clear differences in chromosome II inversion frequencies when compared to Japi population. These two natural populations are distant about 50 km. On the other hand, Japi had frequencies similar to Itatiaia, which is about 250 km distant. Interestingly, both Japi and Itatiaia are located on the same GMR, Atlantic Plateau.
Two simple alternative scenarios (adaptive vs. neutral) can be advanced to explain this finding and other chromosomal inversion frequency geographical differences. In the adaptive case, we expect populations to be panmictic with low genetic structure (high gene flow) and geographical differentiation of chromosomal inversion frequencies correlated to each GMR. Under the neutral case, the observed difference is a casual oscillation caused by forest fragmentation (genetic drift), and we expect to observe populations highly isolated (great genetic structure and low gene flow), great level of linkage disequilibrium in neutral markers and stochastic shifts in chromosomal inversion frequencies (with no pattern).
To test them, we sampled populations in the Campinas Area (Campinas city and Western nearby Capivari town) and in the Eastern Area (Itatiaia and Teresópolis in Rio de Janeiro State; and Juiz de Fora in Minas Gerais State; see Figure 2). Furthermore, we contrasted two kinds of genetic markers-microsatellite markers and chromosomal inversions-to unravel different evolutionary forces-which may shape the observed genetic variation-using estimates such as number of migrants per generation (N m ), F-statistics (F ST , F IS ), and allele frequencies shifts.

| Sampling methods, geographic and climatic variables
We carried out 23 field trips to collect drosophilids from February 2005 to March 2011 in eight localities (Supporting Information Table S1). Five of them were near the boundary of the two GMRs  Eastern Area samples were used since they live in bigger areas of continuous and well-preserved forest, they probably suffer comparatively less effect of genetic drift than of natural selection.
Populations were sampled, when possible, in different seasons trying to avoid biases due to seasonal differences in gene frequencies.
The other populations showed significant differences in their inversion frequencies (with Bonferroni correction): CS (χ 2 = 12.09, df: 3; p = 0.007), CA (χ 2 = 44.21, df: 20; p = 0.0014), and IT (χ 2 = 58.58; df: Costa e Silva population was sampled only twice and showed a clear difference in chromosome inversion frequencies between two consecutive years. So, it is highly desirable to have other samples from this fragment. However, despite the Administration efforts, the place has not been suitable for collections since 2008 due to security reasons beyond our possibilities. Thus, each sample was treated as a different population: CS07 and CS08. In spite of the variations in chromosome inversion frequencies found in the other two populations (CA and IT), we decided to pool all the samples obtained for each population to carry out the statistical analyses since this seems to introduce a smaller a priori bias. Moreover, and most importantly, all the tests we used in the paper (see below) were robust. We repeated the tests using only cold-dry season (Fall-Winter) data, removing hot-rainy season (Spring-Summer) data; furthermore, we also removed the data from 2010 and 2011, that were atypical years (due to the La Niña phenomenon). In all cases, for the statistical tests, we did for chromosome inversion frequencies the results remained qualitatively the same.
We collected drosophilids according to the proceedings described by Batista et al. (2018). Then, we brought the flies to laboratory and sorted them according to external morphology (Freire-Maia & Pavan, 1949;Frota-Pessoa, 1954). We crossed wild-caught males individually with two or three virgin females from the homokaryotypic strain ITC-229-ET, routinely maintained in laboratory conditions (Carvalho, Peixoto, & Klaczko, 1989). We set up isofemale lines from wild females and used F 1 male genitalia for species identification (we compared male F 1 genitalia to drawings described by Frota-Pessoa, 1954).
We obtained the geographical variables (latitude, longitude, altitude- Table 1) using a GPS device. We used data summaries of the nearest meteorological station (average monthly mean temperature and annual sum of monthly precipitation) for each site available at www.agritempo.gov.br (date of access: January 10, 2018) and www.

| Cytological methods and chromosomal inversion frequencies
We prepared slides of polytene chromosomes following a protocol adapted from Ashburner (1989). First, we dissected 3rd instar larvae  We estimated chromosomal inversion frequencies using the "male" and "egg sample" methods (Ananina et al., 2004;Arnold, 1981). We inferred the wild male karyotype using up to eight F 1 karyotyped larvae of the cross between the male and ITC-229-ET females. For the egg sample method, we karyotyped one isofemale F 1 larva that we used for inferring the chromosomal inversion frequencies. We compared both estimated frequencies using a chi-squared test (Zar, 2013), with not a single significant test. So, we pooled them and used the pooled frequencies for statistical analyses. Ananina et al. (2004) pooled inversions DS and DP, because they have similar properties related to seasonality, temperature, precipitation, and altitude. Therefore, as they did, we grouped these inversions in our analysis. All other inversions (DV, DJ, DT, DR) were pooled as "Other" (OT-see in Bitner-Mathé et al., 1995;Ananina et al., 2004).

| Genetic and statistical analyses carried out using chromosomal inversion frequencies
We performed several analyses to characterize population divergence and their genetic structure revealed by inversion frequencies. First, we carried out an exploratory cluster analysis of the populations using Ward's algorithm method based on the minimum variance of Single Euclidian Distance of observed chromosomal inversion frequencies (Paradis, Claude, & Strimmer, 2004).
This is not a test nor a phylogeny, its goal was to obtain a visual representation of the data.
To determine geographical and geomorphological influence on populations, we performed an independent test using the multiple matrix regression with randomization (MMRR). MMRR is essentially a multiple linear regression method applied to matrices (Wang, 2013). It may be used to quantify the association between distance matrices (such as geographic and environmental distances) and a dependent variable, such as genetic distances to make sure the results are not due to a statistical artifact. It is a method to disentangle the relative effects of isolation by distance (IBD) and isolation by environment (IBE), which does not suffer the limitations of the Partial Mantel test (Guillot & Rousset, 2013).
Then, to evaluate cluster heterogeneity, we used a chi-squared test. We estimated hierarchical F ST of the inversion frequencies using Arlequin 3.5 (Excoffier & Lischer, 2010), this is similar to the analysis proposed by Ferrari and Taylor (1981). Populations were considered with low levels of population genetic structure when F ST < 0.05; moderate with 0.05 < F ST < 0.15; and great if F ST > 0.15 (Wright, 1951).
Finally, we tested the correlations between meteorological variables (average monthly mean temperature; annual sum of precipitation) versus genetic data (inversion frequencies after angular transformations Zar, 2013). We used the same procedure for geographical variables (latitude, longitude, altitude, logarithmic transformation of the square root of the area) versus genetic data.

| Microsatellite genotyping
We also used the samples of six populations that were still available for microsatellite genotyping in the two distinct geomorphological regions (GMRs): three on the Atlantic Plateau-Teresópolis (TE n = 32 F 1 females); Itatiaia (IT, n = 32 F 1 females); Colinas do Atibaia (CA, n = 30 F 1 females); and three on the Peripheral Depression-Capivari (CV, n = 32 F 1 females); Santa Genebra (SG, n = 32 females from the field); and Costa e Silva (CS07, n = 34 wild males).
We resolved PCR products in an ABI PRISM 3,500-XL automated sequencer (Applied Biosystems) using GeneScan 600 Liz (Applied Biosystems) as a molecular weight marker. We read genotypes using GeneMarker v.2.2 (SoftGenetics) with manual checking.

| Genetic and statistical analyses carried out using microsatelites
We used GenAlEx 6.5 (Peakall & Smouse, 2012) to infer the estimates of genetic diversity, mean number of alleles per locus (N A ), mean effective number of alleles (N E ) and allele frequencies. Also, we estimated expected heterozygosity of microsatellite loci with Arlequin 3.5. (Excoffier & Lischer, 2010). We inferred null alleles influence in our samples, linkage disequilibrium and tested the alleles from all populations for deviations from Hardy-Weinberg equilibrium using software Genepop on web (Rousset, 2008). We obtained expected heterozygosity parameter from equation H e =1−∑p i 2 , where p i 2 is the total expected frequency of the homokaryotypes for a chromosomal inversion.
We examined genetic structure population revealed by microsatellites using Wright's F-statistics (Wright, 1951)

estimated with
FSTAT (Goudet, 1995) and through an analysis of molecular variance (AMOVA-based on F-statistics) with Arlequin 3.5 (Excoffier & Lischer, 2010), with 1,000 randomization and significant estimates considering p < 0.05. Similar to above, with chromosome inversions, we used the same criteria for determining the levels of population genetic structure. Then, we estimated the effective number of migrants per generation (Nm) by private allele method, implemented on Genepop on web (Rousset, 2008).
We also carried out an MMRR test, as we did with chromosome inversions, using Nei's Genetic Distance for the microsatellite data.
Finally, we tested the correlations between meteorological variables (average monthly mean temperature; annual sum of precipitation) versus genetic data (frequencies of three most common alleles and expected heterozygosity [H e ] after angular transformations Zar, 2013). We used the same procedure for geographical variables (latitude, longitude, altitude, logarithmic transformation of the square root of the area) versus genetic data.

| Population genetic structure revealed by chromosomal inversions
We observed three distinct groups in the cluster analysis ( βENV = 0.454, p = 0.003). These estimates showed that even controlling the geographical distance, the environmental factor (geomorphological region-GMR) was still significant.
Pairwise genetic distances indicated that populations grouped in cluster 1 were less different than others populations within cluster 2 (Supporting Information Table S3). Their genetic distances varied between 0.004 and 0.013, while in the clade with clusters 2 and 3 the distances varied between 0.006 and 0.793. CV and JF were the most distant populations (their distance was 0.793); while SG and CS07 showed the lowest genetic distance (0.004).
Hierarchical F ST analysis of chromosomal inversion genetic distance revealed that more than 91% of the total variation was among individuals within populations with moderate levels of genetic structure among all populations (F ST = 0.088; p = 1 × 10 −5 ). AMOVA also showed that 5% of the total variation can be explained by the

| Clines and correlations of chromosomal inversions polymorphism
We examined the correlations between inversion frequencies and latitude, longitude, altitude, area, average monthly mean temperature, and annual sum of monthly precipitation of the collecting sites (Table 4) We also tested independently the correlations between frequencies of inversions in the second chromosome proximal region, using only wild male karyotypes, with the same meteorological and geographical variables mentioned above (data not shown). The general pattern was consistent with the distal inversions mentioned above.

| Population genetic structure revealed by microsatellites
We observed great variability of microsatellite markers (Table 5) We estimated inbreeding coefficient (F IS ) values for every locus and every population ( We found relatively low levels of population genetic structure (overall F ST estimate of microsatellite markers was F ST = 0.022; p < 0.0001) among these six populations. We estimated genetic differentiation between pairs of populations (pairwise F ST analysis of microsatellites genetic distance- Table 7). Only four out of 15 comparisons were not significant. The smallest difference was between CS07 and SG as well as between IT and CV (F ST = 0.003 ns ; ns: non-significant), while populations most genetically distant were CS07 and CV (F ST = 0.036; p < 0.00001). We observed no significant differences among hierarchical levels microsatellite-based AMOVA (Table 3).
Variation within populations explained about 98% of total variation. On the other hand, only 0.6% of the total variation is due to the variation between the different GMRs and is non-significant (F RT = 0.006 ns ).
Multiple matrix regression with randomization analysis for microsatellite data detected no effect of spatial variation (r 2 = 0.0017; p = 0.993). We observed that geographical (geodesic) and environmental (geomorphological region) distances were not associated with microsatellite both for geodesic distance (βGEOG) and ecologi-

| D ISCUSS I ON
Apart from some beneficial effects over bird guilds (Terraube et al., 2016), habitat loss and forest fragmentation are expected to TA B L E 3 Hierarchical analyses of molecular variance (AMOVA) with fixation indexes due to differences: between geomorphological regions (GMRs; F RT ); among populations within GMRs (F PR ); and among populations in the total (F ST ); and components of the total genetic variance (in percentage) for: between GMRs; within GMRs; and within all populations. AMOVAs estimated using: Inv-chromosomal inversion frequencies; and ssr-microsatellite loci generate small-scale environmental heterogeneity as well as to affect natural populations by drastic reduction of their population size and gene flow (Haddad et al., 2015). These demographic effects can be characterized using genetic markers to assess the genetic structure, effective population size, linkage disequilibrium, dispersion pattern, and changes in allele frequencies (Keyghobadi, 2007).
F-statistics analyses revealed excesses of homozygotes for microsatellites in most samples, with significant F IS . This, probably, is due to an artifact caused by the possible presence of null alleles (Supporting Information Table S2 shows null alleles frequencies per locus in all populations). Some population parameters such as F IS and expected homozygosity may be biased by the presence of null alleles (Chapuis & Estoup, 2007;Waples, 2018 Note. ns, non-significant; *p < 0.05; **p < 0.01; ***p < 0.001 CV: Capivari; SG: Santa Genebra; CS07: Costa e Silva; CA: Colinas do Atibaia; IT: Itatiaia; TE: Teresópolis. population size; therefore, this population does not seem to suffer the negative effects caused by human activity and urbanization (Kenig, Jelic, Kurbalija, Stamenkovic-Radak, & Andjelkovic, 2010).
Their results suggest that divergence in chromosome inversion polymorphisms among Serbian populations of D. subobscura may be an adaptive response to differences among environments they live in (Stamenković-Radak et al., 2012).
Divergence among D. mediopunctata populations may be more associated with climatic and geomorphological properties of both regions than to harmful effects of forest fragmentation. Similarly, divergences in species diversity for Ithomiinae butterflies sampled in the same GMRs may be correlated to geomorphological differences (Brown & Freitas, 2000;. Overall, populations can be grouped in three different clusters according to their genetic distance and their geomorphological location ( Figure 3). The only exception is the population from Costa e Silva, which grouped within clusters 1 and 2 in successive collecting years. The variation in frequencies observed between two consecutive years may be interpreted as result of migration from adjacent areas, since it is located near the border of the two GMR. However, it also could be caused by drastic changes in thermal regime asso- Variations in genetic polymorphisms concomitant with environmental gradients can be considered signs of local adaptation which may lead to population divergence (Wellenreuther et al., 2017). We could not find a consistent pattern of geographical variation for any microsatellite locus. In contrast, inversion polymorphism showed clinal variations congruent with previous findings (Ananina et al., 2004;Batista et al., 2012). This highlights that for a complete understanding of how the fragmentation process is affecting a species, studies should be carried out using different genetic markers to evaluate the joint effects of natural selection, migration and genetic drift.
In summary, our results suggest that, with few exceptions, differences in inversion frequencies of fragmented populations can be maintained according to their geomorphological origin in spite of the effects of gene flow and genetic drift.