Ecological phylogeography and coalescent models suggest a linear population expansion of Anastrepha fraterculus (Diptera: Tephritidae) in southern South America

Abstract

This work is a first approach to an integrated view of the genetics, ecology and dispersion patterns of Anastrepha fraterculus in southern South America. We studied the association of genetic variation with geographical patterns and environmental variables to provide insight into the crucial factors that drive the structure and dynamics of fly populations. Data from a 417 bp mitochondrial COII gene fragment from seven Argentinian populations and one South Brazilian population (from five ecoregions grouped in three biomes) were used to identify population clusters using a model-based Bayesian phylogeographical and ecological clustering approach. The sequences were also analysed under a coalescent model to evaluate historical demographic changes. We identified 19 different haplotypes and two clusters differing in all the environmental covariables. The assumption of neutral evolution and constant population size was rejected, and the population growth parameters suggested a linear population expansion starting 2500 years before present. The most likely ancestral location is Posadas, from where A. fraterculus would have expanded southwards and westwards in Argentina. This result is consistent with Holocene changes and anthropic factors related to the expansion of the Tupí–Guaraní culture, 3000–1500 years before present.

Anastrepha fraterculus, Bayesian analysis, biogeography, COII, Insecta, population genetics, Quaternary climate change, Tephritidae

INTRODUCTION

True fruit flies belong to Tephritidae, one of the largest families of Diptera, with > 4500 species worldwide, including some of the world’s most economically important fruit pests (White & Elson-Harris, 1992; Aluja et al., 2012). Fruit flies use ripening fruits as feeding and breeding substrates. Egg-laying holes and larval development inside the fruit speed up the fruit decay process, severely reducing profits. Economic losses attributable to infestation of fruit flies can be remarkably high (Allwood et al., 2001; Sarwar, 2015; Passos et al., 2018).

The genus Anastrepha (Schiner) is endemic to the Neotropical realm and has been recorded from the southern USA to northern Argentina (White & Elson-Harris, 1992). It is considered the economically most important native fruit fly genus in Latin America (Castañeda et al., 2010), with ~296 species (Norrbom et al., 1999) and 21 species groups (Norrbom et al., 1998; Norrbom & Korytkowski, 2011). The fraterculus group, with 34 species (Prezotto et al., 2019), is as widespread as the whole genus and infests a remarkably diverse group of hosts (Aluja, 1994; Norrbom et al., 1999).

The South American fruit fly, Anastrepha fraterculus (Wiedemann), is widely distributed throughout tropical and subtropical regions, from the southern USA to north and central Argentina (Cladera et al., 2014), and is highly polyphagous, infesting ~170 host species (Hernández-Ortiz et al., 2019), ≥ 29 in Argentina (Oroño et al., 2005). This species together with the invasive Mediterranean fruit fly, Ceratitis capitata (Wiedemann), are the most economically important tephritid species in Argentina (Alberti et al., 2008).

Many studies report that individuals identified as A. fraterculus collected from different regions or hosts (see references listed by Cáceres et al., 2009) show clear-cut differences in terms of several lines of evidence; these include adult and egg morphology, mating compatibility, hybridization, host preference, classical and molecular cytogenetics, and isozymes (Cladera et al., 2014). Also, phylogenetic inferences based on molecular markers (Rocha & Selivon, 2002; Prezotto et al., 2019) and DNA sequencing (Steck & Sheppard, 1993; McPheron et al., 1999) have been reported. On these grounds, many authors support the view that the Latin name A. fraterculus is applied to a cryptic complex of different biological species (AF complex) (Morgante et al., 1980; Hernández-Ortiz et al., 2004, 2012; Cladera et al., 2014; Dias et al., 2016). The presence of multiple gene pools and the lack of monophyly in this nominal species were corroborated by Smith-Caldas et al. (2001).

Although the number of putative species within the AF complex is not yet known, Hernández-Ortiz et al. (2004, 2012, 2015) identified at least eight mostly allopatric morphotypes: ‘Mexican’ (México, Guatemala and Panamá), ‘Venezuelan’ (lowlands of Venezuela), ‘Andean’ (highlands of Venezuela and Colombia), ‘Peruvian’ (lowlands of Perú and Ecuador; A. sp. 4), ‘Ecuadorian’ (highlands of Perú and Ecuador) and three Brazilian morphotypes, ‘Brazilian-1’ (Argentina and Brazil), ‘Brazilian-2’ and ‘Brazilian-3’ (corresponding to A. sp. 1, A. sp. 2 and A. sp. 3 in the studies by Selivon et al., 2002, 2005). In a recent study, the genetic divergence of ITS1 sequences was correlated with morphological differences for six of these morphotypes (Prezotto et al., 2019).

The existence of various biological entities within the AF complex might be related to its wide continental distribution, encompassing diverse ecological conditions and geographical areas where they occur (Prezotto et al., 2019). It is possible that A. fraterculus populations became structured in response to adaptation to different biomes. It is well known that morphological variability within the AF complex appears along with different ecological conditions and contrasting host preferences (Hernández-Ortiz et al., 2012).

In many phylogeographical studies (Avise, 2004) mitochondrial genetic markers have proved informative for population genetic studies owing to their strictly maternal inheritance, lack of recombination, relatively high mutation rate and the availability of very efficient polymerase chain reaction primers (Gallo-Franco et al., 2017). These markers have been used successfully in tephritids for resolving interspecific relationships (Barr & McPheron, 2006; Boykin et al., 2006; Vaníčková et al., 2015) and intraspecific structure (McPheron et al., 1994; Alberti et al., 2008; Lanzavecchia et al., 2008), developing species identification, diagnostic methods (Barr & McPheron, 2006) and pathway analysis of populations (Lanzavecchia et al., 2008; Ruiz-Arce et al., 2012).

Previous studies on morphometry (Hernández-Ortiz et al., 2012), mating behaviour (Petit-Marty et al., 2004; Vera et al., 2006), ethology (Petit-Marty et al., 2004), cytogenetics (Basso et al., 2003) and mitochondrial COII (Alberti et al., 1999, 2002) and COI (Dias et al., 2016) DNA sequences conducted in Argentina suggest that A. fraterculus populations from this country belong to a single biological species (Alberti et al., 2008) and, according to Rull et al. (2012), match the Brazilian-1 morphotype. Based on the mitochondrial COII gene sequence clustering and analysis of genetic and geographical distances, Alberti et al. (2008) reported lack of association between haplotypic variability and geographical distribution. In contrast, populations of A. fraterculus seem to be homogeneous within smaller areas (Ludueña et al., 2010).

Currently, eco-genetic models, developed to understand population structure and the forces modulating genetic variability distribution, represent a practical application of theoretical evolutionary genetics that can contribute to solve problems of pest management in agricultural settings (Ruiz-Arce et al., 2019). However, despite their global importance, only a few publications have seriously examined the long-term population dynamics in fruit flies (Aluja et al., 2012). An in-depth knowledge of these dynamics and the history of the target species is important for pest management because it improves the ability to predict how populations will react to environmental changes (Porretta et al., 2007; Aluja et al., 2012; Wei et al., 2015).

The current population structure and distribution of natural populations of pest species result from the interaction between their ability to invade based on environmental requirements and geographical variation of environmental features and historical events. In this context, research on fine-scale population structure and genetic variation patterns at larger spatial and temporal scales can increase the understanding of interactions between micro-evolutionary processes and biotic and abiotic conditions. A substantial contribution is provided by phylogeographical approaches that focus on the processes governing the distribution of genealogical lineages within species across the geographical landscape (Arbogast & Kenagy, 2001; da Silva et al., 2018).

Within this framework, simulation modelling is a promising approach to develop robust, spatially explicit evolutionary landscape genetic hypotheses, offering an opportunity to understand the range of conditions and causes that can generate complex patterns. This makes it possible to validate those hypothesis by using the processes identified as important for the empirical patterns as input for the simulations (Balkenhol et al., 2015).

The present study applies a phylogeographical and ecological approach to analyse the distribution of variability of a mitochondrial COII gene fragment in samples from Argentinian and South Brazilian A. fraterculus populations. Our main goal was to improve the understanding of the micro-evolutionary processes related to genetic diversity and population structure, in addition to the spatial and temporal scales of divergence and migration patterns of this species in the southernmost portion of its distribution. We intended to assess the influence of environmental conditions on the demographic history of A. fraterculus by using a Bayesian ecological clustering approach and simulations based on the theory of coalescence. The hypotheses to be tested were that the distribution of A. fraterculus populations is affected by environmental factors and that colonization of the southernmost populations is relatively recent. The predictions emanating from these hypotheses are that: (1) A. fraterculus underwent a significant population expansion in its southern range; (2) the relationships between populations differ from the expectations based on the isolation-by-distance model; and (3) environmental differences contribute to explain geographical patters of genetic variation.

MATERIAL AND METHODS

Population sampling and environmental variables

Sampled populations belong to five ecoregions grouped in three biomes (Table 1; Fig. 1). A detailed description of ecoregions in Argentina is provided by Morello et al. (2012), and distribution maps of world terrestrial ecoregions and biomes are given by Olson et al. (2001) and Dinerstein et al. (2017). Climatic and normalized difference vegetation index (NDVI) data were obtained respectively from WorldClim v.2 (https://www.worldclim.org/data/worldclim21.html) and NASA Earth Observations (NEO) (https://neo.sci.gsfc.nasa.gov/view.php?datasetId=MOD_NDVI_M). Data related to bioclimatic variables and solar radiation correspond to the mean for the period 1970–2000, represented as GeoTiff files. NDVI data corresponded to the year 2001. In the case of solar radiation and NDVI data, 12 files were downloaded, one for each month of the year. In all cases, the spatial resolution of GeoTiff files is 30′ (~1 km²). Six environmental variables were included as covariates in the DNA analysis described below: (1) mean temperature of warmest quarter (TWQ); (2) mean temperature of coldest quarter (tcq); (3) precipitation of wettest quarter (PPWQ); (4) precipitation of driest quarter (ppdq); (5) average solar radiation (SRAD); and (6) average NDVI (NDVI). Temperature variation throughout the year is summarized by TWQ and tcq, and the same holds for PPWQ and ppdq in reference to rainfall. The SRAD and NDVI values were averaged over months to obtain a single value. Environmental variable values at each sampled location (Table 1) were extracted from the GeoTiff files with the function extract of the package raster (Hijmans, 2019) of the software R v.3.6.0 (R Core Team, 2019).

Table 1.

Geographical coordinates and environmental data for the sampled populations of Anastrepha fraterculus

Sample	Province/state	Country	Biome*	Ecoregion†	N	Host plant	Longitude	Latitude	(1) TWQ (°C)	(2) tcq (°C)	(3) PPWQ (mm)	(4) ppdq (mm)	(5) SRAD (kJ m−2 day−1)	(6) NDVI
YU	Jujuy	Argentina	1	Yungas	3	Y	−64.45	−23.63	26.71	16.30	484.75	20.50	16615.29	203.54
HM	Tucumán		1	Yungas	3	P	−65.33	−26.80	24.60	13.58	462.00	40.75	16455.56	216.73
PO	Misiones		1	Paraná forest	9	Y	−55.87	−27.38	26.27	16.75	521.25	342.25	16104.58	155.56
CA	Buenos Aires		8	Pampa	3	G	−58.63	−34.65	23.40	11.34	329.00	174.75	16040.00	130.22
MR			8	Pampa	2	P, L, M	−58.37	−34.83	22.90	10.89	324.25	178.50	16324.35	156.99
ME	San Luis		7	Dry Chaco	3	Y	−65.03	−32.35	22.43	10.20	289.00	55.00	17092.21	161.14
CO	Entre Ríos		8	Espinal	3	Y	−58.15	−31.03	25.03	13.24	432.75	195.00	16832.17	196.39
PE	Rio Grande do Sul	Brazil	7	Uruguayan savanna	2	Y	−52.35	−31.77	24.23	14.45	394.00	278.25	14394.21	174.31

N = sample size. Sampling sites: CA, Castelar; CO, Concordia; HM, Horco Molle; ME: Merlo; MR, Ministro Rivadavia; PE, Pelotas; PO, Posadas; YU, Yuto. Host plant acronyms: G, green guava (Feijoa sellowiana) (Berg.); L, plum (Prunus salicina) Lindl.; M, medlar (Eriobotrya japonica) (Thunb.); P, peach (Prunus persica) L.; Y, yellow guava (Psidium guayaba) L. Environmental variables: (1) mean temperature of warmest quarter (TWQ); (2) mean temperature of coldest quarter (tcq); (3) precipitation of wettest quarter (PPWQ); (4) precipitation of driest quarter (ppdq); (5) average solar radiation (SRAD); and (6) average normalized difference vegetation index (NDVI).

*Biome numbers according to Dinerstein et al. (2017) (see Fig. 1): 1, tropical and subtropical moist broadleaf forests; 7, tropical and subtropical grasslands, savannas and shrublands; 8, temperate grasslands, savannas and shrublands.

†Descriptions of ecoregions are provided by Morello et al. (2012) (Argentinian ecosystems) and https://www.worldwildlife.org/ecoregions/nt0710 (Uruguayan Savanna).

Figure 1.

Partial map of South America, representing the distribution of different biomes and the locations sampled for Anastrepha fraterculus. A, biomes are represented in different colours and with numerical codes according to Dinerstein et al. (2017): (1) tropical and subtropical moist broadleaf forests; (2) tropical and subtropical dry broadleaf forests; (4) temperate broadleaf and mixed forests; (7) tropical and subtropical grasslands, savannas and shrublands; (8) temperate grasslands, savannas and shrublands; (9) flooded grasslands and savannas; (10) montane grasslands and shrublands; (12) mediterranean forests, woodlands and scrub; and (13) deserts and xeric shrublands. B, distribution of clusters identified in the phylogeographical and ecological analysis. The light blue and yellow areas represent the probable ranges of clusters 1 and 2, respectively. The arrowhead indicates the most probable root population. Population labels are coloured according to their corresponding biomes (1, 7 and 8). Sampling sites: CA, Castelar; CO, Concordia; HM, Horco Molle; ME, Merlo; MR, Ministro Rivadavia; PE, Pelotas; PO, Posadas; YU, Yuto.

The recorded environmental variables were chosen because they are considered meaningful for modelling A. fraterculus distribution. Temperature has been shown to affect A. fraterculus life-history traits (Cardoso et al., 2002). Montoya et al. (2008) reported effects of rainfall and soil moisture on survival of both immature and adult stages of Anastrepha ludens (Loew) and Anastrepha obliqua (Macquart), and Toledo et al. (2014) observed that these environmental factors also affect the ability of an entomopathogenic nematode (Heterorhabditis bacteriophora Poinar) to infect A. ludens. The SRAD has been claimed to affect insect development and the ability to survive (Querci & Romei, 1945a, b; Battisti et al., 2013). The NDVI is a simple coefficient (without dimensions), widely used to quantify healthy vegetation, estimated from spectral reflectance at the visible (VIS) and near-infrared (NIR) bands according to the expression: NDVI = (NIR − VIS)/(NIR + VIS) (Tucker et al., 1983; Townshend & Justice, 1986). This coefficient has been helpful in a niche modelling framework to produce predictions of the annual and seasonal distributions of suitability for Ceratitis capitata, the medfly, and identifying regions that represent potential risk areas for invasions (Szyniszewska & Tatem, 2014; Szyniszewska et al., 2016). It has also been used to assess the areas of poor growth and health of host plants attributable to Hessian fly [Mayetiola destructor (Say)] infestation (Bhattarai et al., 2019) and proved to be one of the environmental variables that could overshadow spatial distance as a driver of bat fly (Diptera: Streblidae) species composition (Eriksson et al., 2020).

DNA fragment analysis

The analysis was based on a 417 bp fragment of mitochondrial DNA corresponding to the COII gene described by Alberti et al. (2008). The sample involved 28 individuals from seven Argentinian populations and one Brazilian population (Table 1; Fig. 1). The corresponding COII gene fragment sequences were obtained from GenBank (Supporting Information, Table S1).

All analyses were conducted with different packages of the software R v.3.6.0 (R Core Team, 2019). Sequence files in fasta format were imported and converted into DNAbin objects. Diversity among sequences was quantified by the mean number of pairwise mismatches (Π), nucleotide diversity (π) (Nei & Li, 1979) and haplotype diversity (h). The value of Π was estimated with the function dist.dna of the package ape (Paradis & Schliep, 2018), whereas π and h were obtained with the functions nucleotideDiversity and heterozygosity of the package stratag (Archer et al., 2016). Pairwise between-population diversities (d_XY) and divergences (d_A) were estimated by applying the evolutionary model of Kimura (1980), with the function nucleotideDivergence of the same package.

DNAbin object was converted into a genind object using the function DNAbin2genind of the package adegenet (Jombart, 2008; Jombart & Ahmed, 2011), which retains only polymorphic sites. From this format, Nei’s genetic diversities (H) of each population were estimated with the function poppr of the package poppr (Kamvar et al., 2014, 2015). Pairwise Nei’s (1978) genetic distances between populations (D) were estimated with the function dist.genpop of the package adegenet. A phenogram with bootstrap support (100 pseudoreplicates) was obtained with the function aboot of poppr. The hypothesis of isolation by distance was evaluated by comparing pairwise population geographical distance with nucleotide divergence matrices by Mantel test (with 2000 permutations) using the package ade4 (Dray et al., 2007).

Population structure was quantified by Φ statistics with the package stratag using the function statPhist (applying the model of Kimura 1980, Φ _K80) and by analysis of molecular variance with the function poppr.amova of the package poppr (Φ _AMOV).

Demographic changes were evaluated graphically by a mismatch distribution plot produced with the function MMD of the package pegas (Paradis, 2010) and by calculating the raggedness index (rg) (Harpending, 1994), the neutrality Tajima’s D (TD) (Tajima, 1989) and Fu’s F_S (Fu, 1997) tests. The statistic rg was calculated with an ad hoc script, whereas TD and F_S were obtained with the functions tajimasD and fusFS, respectively, of stratag. The statistical significance of rg, TD and F_S was determined by comparison with the neutral expected distribution obtained from coalescent model simulations generated with the function coal_model of the package coala (Staab & Metzler, 2016). A total of 1000 sets of 28 DNA sequences corresponding to the fragment analysed were simulated (28 sequences, one locus, 417 bp), with a mutation rate equivalent to the observed Π estimate. From each set of simulated DNA sequences, estimates of TD, F_S and rg were obtained (TD_sim, F_Ssim and rg_sim, respectively). In each case, the observed coefficients (TD, F_S and rg) were considered significant if they fell outside the 95% confidence interval generated by the coalescent simulations.

To determine the most probable time-dependent demographic model, a coalescent analysis of the COII fragment sequences was conducted by Monte Carlo simulations of Markov chains (MCMC) with the package coalescentMCMC (Paradis, 2015). The method implemented in this package uses neighbourhood rearrangement, following Kuhner et al. (1995), and allows for testing different models. In our study, we compared four demographic models: (1) constant Θ(4Nµ), where Θ is the effective population size times mutation rate per site, N is the population size and μ is the mutation rate; (2) exponential; (3) step; and (4) linear growth. For each model, 15 000 simulated (phylogenetic) trees were obtained, from which the last 10 000 were kept (burnin = 5000).

The number of estimated demographic parameters depended on the model. The first one (constant model) assumed that Θ remained without changes through time, and therefore, the analysis provided only the estimation of Θ.

The exponential model assumed that Θ varied through time according to the expression:

$${\displaystyle \begin{array}{l}{\displaystyle {\mathrm{\Theta }}_{\left(t\right)}={\mathrm{\Theta }}_0{\mathrm{e}}^{\mathrm{\rho }t}}\end{array}}$$

where Θ ₀ and Θ _(t) are the values of Θ at present and at t time units before present, respectively, and ρ is the growth rate. In this case, two parameters were estimated, Θ ₀ and ρ.

The step model assumed two constant values of Θ, denoted as Θ ₀ and Θ ₁, as follows:

$${\displaystyle \begin{array}{l}{\displaystyle {\mathrm{\Theta }}_{\left(t\right)}=\{\begin{array}{l}{\mathrm{\Theta }}_{0}t\le \mathrm{\tau }\\ {\mathrm{\Theta }}_{1}t>\mathrm{\tau }\end{array}}\end{array}}$$

where τ is a point in time to be estimated. In this case, three parameters were estimated: Θ ₀, Θ ₁ and τ.

Finally, in the linear model Θ varied through time according to the expression:

$${\displaystyle \begin{array}{l}{\displaystyle {\mathrm{\Theta }}_{\left(t\right)}={\mathrm{\Theta }}_0+\frac{({\mathrm{\Theta }}_{\mathrm{T}\mathrm{M}\mathrm{R}\mathrm{C}\mathrm{A}}-{\mathrm{\Theta }}_0)}{T_{\mathrm{M}\mathrm{R}\mathrm{C}\mathrm{A}}}t}\end{array}}$$

where Θ _TMRCA is the value of Θ at the time to the most recent ancestor, denoted as T_MRCA. In this case, the estimated parameters were Θ ₀, Θ _TMRCA and T_MRCA.

Models were compared according to their deviance information criterion (DIC) (Spiegelhalter et al., 2002), defined as:

$${\displaystyle \begin{array}{l}{\displaystyle \mathrm{D}\mathrm{I}\mathrm{C}={\displaystyle \overline{\mathrm{D}\mathrm{e}\mathrm{v}}}-\mathrm{D}\mathrm{e}{\mathrm{v}}_{({\displaystyle \overline{\mathrm{\Theta }}})}}\end{array}}$$

Where $\overline{\mathrm{D}\mathrm{e}\mathrm{v}}$ is the deviance averaged over the simulated trees and $\mathrm{D}\mathrm{e}{\mathrm{v}}_{(\overline{}\theta )}$ the deviance for the averaged estimated parameters (see vignettes provided by Paradis, 2015).

To reveal possible migration events, we analysed the distribution of the COII fragment sequence variability by means of a Bayesian phylogeographical and ecological clustering approach, including geographical and environmental variables as potentially explanatory factors, using the package BPEC (Manolopoulou & Hille, 2016). For the resulting clusters, population structure was evaluated by an AMOVA, as described above.

Sequences were read from a nexus format file with the function bpec.loadSeq to obtain the list of sequences. The covariates included in the analysis were the geographical coordinates and the six environmental variables indicated in Table 1. The phylogeographical analysis was conducted by the function bpec.mcmc, with two main arguments, the list of sequences and the data frame with geographical and environmental variables. The conditions for this function were 100 000 iterations, saving 10 000 posterior samples, with a parsimony relaxation parameter, ds = 0. To avoid bias owing to a scale effect, environmental variables were scaled to mean = 0 and variance = 1. The algorithm simultaneously draws inferences about the genealogy (haplotype network) and the population clustering and identifies locations with high posterior probability of being ancestral. The haplotype network was obtained with the function bpec.treePlot, and the probable geographical distribution of the clusters identified was plotted with the function bpec.contourPlot.

RESULTS

Gene diversity and population structure

The mean nucleotide frequencies in the sequence analysed have been reported previously by Alberti et al. (2008). Out of the 417 bp, 34 (8.15%) sites were variable, with two (31 sites) or three alleles (three sites), in the total sample of 28 individuals. The transition-to-transversion ratio (13:27) fitted the expectation (1:2) if all substitutions were equally probable (χ ² = 0.01, P = 0.91). Among the 71 alleles, 19 were present in single copies (singletons). The mean number of pairwise mismatches between sequences was Π = 4.87.

The total number of haplotypes in the whole sample was 19. All populations apart from PE had more than one haplotype, yielding high haplotype diversity within populations (Table 2). Nucleotide diversity within each population (π _i) ranged from zero to 0.02, with a mean value $\overline{\pi }=0.008$⁠. The global nucleotide diversity (π = 0.012) was slightly higher than the average within-population $\overline{\pi }$⁠, indicating low nucleotide divergence among populations (d_A). Values of pairwise d_A were, in most cases, close to zero (Table 3). In congruence with the low estimates for d_A, both statistics of population structure were low and non-significant: Φ _K80 = −0.08 (P = 0.87) and Φ _AMOV = −0.01 (P = 0.51).

Table 2.

Haplotype diversity, nucleotide diversity and Nei’s genetic diversity within populations of the COII mitochondrial DNA fragment analysed in Anastrepha fraterculus

Population	Haplotype diversity	Nucleotide diversity	Nei’s diversity
CA	1.000	7.19 × 10−3	0.098
CO	1.000	1.28 × 10−2	0.157
HM	1.000	4.80 × 10−3	0.059
ME	1.000	6.39 × 10−3	0.078
MR	1.000	7.19 × 10−3	0.088
PE	0.000	0.000	0.000
PO	0.972	2.12 × 10−2	0.260
YU	0.667	1.60 × 10−3	0.020
Total	0.915	0.012	0.148

Nei’s diversities were estimated considering only variable sites. Sampling sites: CA, Castelar; CO, Concordia; HM, Horco Molle; ME, Merlo; MR, Ministro Rivadavia; PE, Pelotas; PO, Posadas; YU, Yuto.

Table 3.

Nucleotide diversity (above diagonal) and divergence (below diagonal) between populations estimated according to Kimura (1980)

	CA	CO	HM	ME	MR	PE	PO	YU
CA	–	9.99 × 10−3	6.00 × 10−3	6.00 × 10−3	7.19 × 10−3	3.60 × 10−3	1.37 × 10−2	3.60 × 10−3
CO	< 10–19	–	8.26 × 10−3	1.04 × 10−2	9.99 × 10−3	6.39 × 10−3	1.98 × 10−2	7.19 × 10−3
HM	< 10–19	−5.33 × 10−4	–	6.39 × 10−3	6.00 × 10−3	2.40 × 10−3	1.57 × 10−2	3.20 × 10−3
ME	−7.99 × 10−4	7.99 × 10−4	7.99 × 10−4	–	6.79 × 10−3	4.00 × 10−3	1.56 × 10−2	4.80 × 10−3
MR	−8.67 × 10–19	< 10–19	< 10–19	−8.67 × 10–19	–	3.60 × 10−3	1.59 × 10−2	4.40 × 10−3
PE	< 10–19	< 10–19	< 10–19	7.99 × 10−4	−4.34 × 10–19	–	1.36 × 10−2	7.99 × 10−4
PO	−4.66 × 10−4	2.82 × 10−3	2.73 × 10−3	1.84 × 10−3	1.67 × 10−3	3.00 × 10−3	–	1.33 × 10−2
YU	−7.99 × 10−4	8.67 × 10−19	< 10−19	7.99 × 10−4	< 10−19	< 10−19	1.93 × 10−3	–

Sampling sites: CA, Castelar; CO, Concordia; HM, Horco Molle; ME, Merlo; MR, Ministro Rivadavia; PE, Pelotas; PO, Posadas; YU, Yuto.

Genetic distances between populations are best represented by restricting the analysis to variable sites. The unweighted pair group method with arithmetic mean (UPGMA) tree based on Nei’s genetic distances (Fig. 2) showed that populations were not grouped according to geographical distances. The correlation between geographical and genetic distances was negative and non-significant (r = −0.47, P = 0.995), indicating that genetic differentiation between populations did not fit the isolation-by-distance model.

Figure 2.

Phenogram obtained by unweighted pair group method with arithmetic mean (UPGMA) from Nei’s genetic distances between populations of Anastrepha fraterculus. The numbers on the nodes represent the bootstrap support based on 100 pseudoreplicates. See legend to Figure 1 for the names of the sampling sites.

Demographic change analysis

Raggedness was borderline significantly lower than the expectation for a simulated stable population (rg = 0.022, P = 0.098). The mismatch distribution plot (Fig. 3) differed from the stable expectation, showing a high peak around two and a small secondary peak around 16. Both Tajima’s and Fu’s neutrality tests were negative and significant (TD = −1.57, P = 0.035; F_S = −7.87, P < 0.001).

Figure 3.

Mismatch distribution plot obtained from 28 sequences of a 417 bp mitochondrial COII fragment in Anastrepha fraterculus.

Taken as a whole, the results of TD, F_S, rg and the mismatch distribution plot suggested that the studied sample did not fit the assumption of neutral evolution and constant population size. To gain a deeper insight into possible demographic changes, we estimated Θ and related parameters from coalescent trees simulated by MCMC procedures (Fig. 4). The three plots in Figure 4 compare the changes of Θ through time with the constant model, where Θ = 0.090 and DIC = −0.28. Assuming the linear model, Θ grew from 0.108 to 0.976 from T_MRCA = 0.019 units of time before present, with DIC = −42.78. The exponential model, in contrast, suggested that Θ has been reduced with a rate ρ = 1, and DIC = −0.46. Finally, parameters estimated for the step model suggested that Θ increased from 0.026 to 0.664 in a moment situated 0.024 time units before present, with DIC = −13.86.

Figure 4.

Expected variation of Θ (effective population size times mutation rate per site) through time for the linear, exponential and step population growing models. In all three plots, the predicted line is compared with the constant model (horizontal dashed line). In the linear model plot, the vertical dashed line represents the time to the most recent common ancestor (T_MRCA).

The two best-supported models, linear and step, were consistent in suggesting population expansion starting ~0.02 time units before present. The main difference between these models was the dynamics of the change, steep or gradual. According to DIC estimates, the best-fitting model was the linear one. This result was consistent with the negative and significant TD and F_S values, compatible with a population expansion. Assuming a divergence rate of 2.3% (Brower, 1994), T_MRCA represented ~20 000 generations. Given that A. fraterculus is estimated to have eight generations per year (Pereira dos Santos et al., 2017), the results suggested that populations of this species from Argentina and southern Brazil have been expanding during the last 2500 years.

Phylogeography and ecological clustering

The ecological clustering identified two clusters, represented in yellow and light blue in Figure 5, and the most probable number of migrations was one (P = 1). The most likely root node is extinct (haplotype 39 in the left-hand network shown in Fig. 5). Cluster 1 (light blue) included five haplotypes found in Posadas: P1, P3, P6, P9 and P10. Cluster 2 (yellow) involved 12 haplotypes. Haplotype P2 belongs to cluster 1 or 2 with the same probability (Fig. 5). One haplotype from Concordia (Co8) could not be assigned by BPEC to any of the clusters. The AMOVA results showed non-significant structure at the ecological clustering level (P < 0.05).

Figure 5.

Haplotype network of the mitochondrial COII fragment of Anastrepha fraterculus obtained by the Bayesian phylogeographical analysis conducted with the software package BPEC. On the left side, nodes are represented by pies, in which yellow or light blue indicates the probability of membership of the identified clusters (see Fig. 1), whereas pink indicates missing intermediate haplotypes. On the right side, the nodes are labelled with the haplotype names given by Alberti et al. (2008).

The posterior distribution of the covariates showed highly significant differences between the two clusters for all environmental variables (t varied from ten to 431, P < 10⁻¹⁵; Fig. 6). Considering geographical locations and environmental variables, the most likely ancestral location was Posadas (P = 0.56; cluster 1; Fig. 1). The evidence from the phylogeographical and ecological clustering suggested that the populations analysed in this study from central northwest Argentina and south Brazil might have been derived from the area of Posadas and expanded southwards and westwards, resulting in the populations included in cluster 2.

Figure 6.

Boxplots representing the posterior distribution of the environmental variables (mean and SD) in the two clusters identified by the ecological clustering. See footnote to Table 1 for the full names of the environmental variables.

DISCUSSION

Anastrepha fraterculus is considered a complex of synmorphic species with at least eight morphotypes (Hernández-Ortiz et al., 2004, 2012, 2015), three of them reported in sympatry in Brazil (Selivon et al., 2002; Selivon & Perondini, 2007; Vaníčková et al., 2015). In order to differentiate synmorphic species and to analyse genetic structure, discriminate geographical collections and measure variation within and between populations of fruit flies and other insects, molecular markers are useful (Ruiz-Arce et al., 2012; Zhao et al., 2017). Previous studies based on the phylogenies and geographical range of samples, concluded that COII revealed more variability than COI and would be a better marker to distinguish cryptic species in the A. fraterculus complex (Ludueña et al., 2010; Vaníčková et al., 2015). Previous studies, based on COI and COII (Alberti et al., 2008; Ludueña et al., 2010; Dias et al., 2016), morphological traits (Hernández-Ortiz et al., 2012) and mating compatibility (Petit-Marty et al., 2004; Vera et al., 2006), have shown populations from Argentina and south Brazil to be highly related and to belong to the same morphotype, known as Brazilian-1 (Rull et al., 2012). However, Prezotto et al. (2019), using ITS1 and morphometric data, revealed differences in Argentinian samples of such magnitude that it raises questions about the occurrence of a single entity in this country.

In the present work, the percentage of polymorphic sites and global haplotypic diversity (h) were high when compared with other tephritid species (Ruiz-Arce et al., 2015; Dias et al., 2016; Gallo-Franco et al., 2017). Regardless of the molecular marker selected, low genetic variability has often been observed in Anastrepha, and the underlying causes for this trend remain unknown (Passos et al., 2018; Prezotto et al., 2019). In our study, the nucleotide diversity (π) within populations and divergence among populations (d_A) were low, but within the range recorded using COI + ND6 mitochondrial DNA regions in other species of the genus, such as Anastrepha striata Schiner (π = 0.0005; Gallo-Franco et al., 2017), A. obliqua (π = 0.012–0.016; Ruiz-Arce et al., 2012; Aguirre‐Ramirez et al., 2017) and A. ludens (π = 0.02; Ruiz-Arce et al., 2015). Differences in genetic diversity between species might be associated with dietary habits or sample strategies (Gallo-Franco et al., 2017; Zhao et al., 2017). Low values of π and high values of h, as observed in our case, suggest the existence of small populations that have undergone recent population growth (Rosetti & Remis, 2012).

In congruence with the low values of d_A, both statistics of population structure (Φ _K80 and Φ _AMOV) were low both between populations and between ecological clusters. The UPGMA tree based on Nei’s genetic distances did not group populations according to their geographical provenance, in agreement with our working hypothesis. Lack of association between genetic and geographical distances has also been observed in other tephiritid species, where genetic differentiation among populations does not fit the isolation-by-distance model expected for neutral markers (Karsten et al., 2013; Gallo-Franco et al., 2017; Passos et al., 2018). Some of the hypotheses regarding this topic are the existence of high gene flow among populations, probably associated with commercial activities (Hill, 2008; Alberti et al., 2008; Ruiz-Arce et al., 2015), the overlapping distribution of populations (Karsten et al., 2013) and the absence of physical barriers, which allows constant unrestricted movement of A. fraterculus across the natural regions under study.

Mild population structuring is commonly observed in flying insect species, specifically in those that are good dispersers (da Silva et al., 2018). However, different morphotypes have been found in Brazil occurring in sympatry (Selivon et al., 2002; Selivon & Perondini, 2007; Vaníčková et al., 2015), and cryptic population structure has been described in Argentinian populations (Oroño et al., 2013; Rodriguez et al., 2019), suggesting that variability in resource characteristics within a site might contribute to genetic discontinuities. Regarding this issue, our results indicate that A. fraterculus was able to thrive along different ecoregions in a short period of time, indicating good dispersal ability and high gene flow between populations, preventing geographical differentiation. This feature is also consistent with one of our predictions, a recent invasion by a limited number of founders, followed by population growth (Gallo-Franco et al., 2017) and expansion (Malacrida et al., 2007).

Our results suggest a recent colonization by A. fraterculus in Argentina, from the north-eastern region towards the west and south. Recently, new mathematical models based on coalescent theory have enhanced the analysis of the demographic history of natural populations, and estimates of genetic diversity and divergence time allow dispersal and demographic processes to be inferred (Rosetti & Remis, 2012). In other tephritid species, estimates of Fu’s F_S and Tajima’s D were negative, and the raggedness index was also consistent with the hypothesis of population expansion (Ruiz-Arce et al., 2012; Karsten et al., 2013; Low et al., 2014; Gallo-Franco et al., 2017; Passos et al., 2018). Our findings are consistent with these studies in showing a genetic signature of demographic expansion.

Our haplotype network shows a star-like shape, with the presence of a single common haplotype and a number of rare haplotypes connected to it by successional mutation steps. This pattern is consistent with the demographic analysis. The short divergence time would not have enabled mutations to occur and accumulate, despite the ecological or geographical isolation of populations (Avise, 2000). Similar results were found in Simulium tani Takaoka & Davies (Low et al., 2014) and A. striata (Gallo-Franco et al., 2017), in agreement with the observation that although most haplotypes are not shared among locations, they are still found at very low frequencies. The conclusion from the demographic analysis, suggesting a fast expansion for our populations, might explain the low genetic differentiation and the lack of isolation by distance.

Knowledge and understanding of the effect of past events is important to assess possible factors that shape population dynamics. Fly population distribution and genetic variability are typically associated with environmental variables, such as temperature, precipitation and soil moisture (Montoya et al., 2008; Battisti et al., 2013; Ma et al., 2020). Therefore, studying the effect of both global and local environmental patterns can provide significant insight into the crucial factors that drive the ecology and dynamics of insect populations at different scales (Aluja et al., 2012). In addition, in the context of climatic change, it is useful to investigate rates of adaptation in pest populations and to predict how their ability to invade might be altered (Larson et al., 2019; Ma et al., 2020).

In the present study, the estimated population expansion time, assuming a divergence rate of 2.3% (Brower, 1994) and eight generations per year (Pereira dos Santos et al., 2017), was 2500 years before present (BP). Moreover, our analysis located the most likely ancestral population in an area around Posadas. The sample from this site showed the highest variability, which is expected to occur in the centre from which the expansion originated (Ruiz-Arce et al., 2015). This result is consistent with a recent study by Giardini et al. (2020), based on chromosomal polymorphism, which supports one possible route of invasion of A. fraterculus sp. 1 from south-east Brazil, followed by a wide dispersion of this pest throughout the Argentinian territory. Additionally, our work provides valuable information about the timing of the demographic events.

Several studies in subtropical north-east Argentina (Misiones) and south-east Brazil, based on pollen (Behling, 2002), vegetation and fire-history records (Gessert et al., 2011), sediment/palaeosol sequences (Zech et al., 2009b), precipitation changes (Smith & Mayle, 2018) and palaeoclimate and archaeological analyses (Iriarte et al., 2017) reveal a humid climate period in the Late Holocene, starting 3000 years BP. This scenario, when climate became wetter again owing to the restrengthening of the South American monsoon (Zech et al., 2009a, b), enabled the expansion of humid forests at the expense of more drought-tolerant savanna/grassland ecosystems (Gessert et al., 2011; Iriarte et al., 2017). This context might have contributed to A. fraterculus population expansion. In fact, climatic conditions becoming favourable, with a warmer and more humid climate, have been reported to be an important driver in genetic dynamics of insect populations (Meeyen et al., 2014; Gallo-Franco et al., 2017).

Associated with this forest expansion, multiple sources suggest that the Tupí–Guaraní-speaking people split and began to spread southwards from Amazonia to south-east South America ~3000 years BP (Iriarte et al., 2017). It is likely that increasing precipitation through the mid- to late Holocene would have led to the expansion of gallery forests and increased the upper Paraguay River and Paraná River levels, allowing the Tupí–Guaraní-speaking people to spread and establish villages. Increases in vegetation score and colonization by the Tupí–Guaraní-speaking people occurred later as one progresses further south, dating to 2500–2000 years BP in the Upper Paraná, ~2000 years BP in the Middle Paraná (27°S, were Posadas is located) and ~2000–1500 years BP down the Paraná (Iriarte et al., 2017). We suggest that the demographic growth of this culture might have been a driver of A. fraterculus expansion by human-mediated dispersal.

In our work, we introduced environmental conditions when performing BPEC; in this way, we identified two clusters that showed significant differences for all environmental variables. Cluster 1 comprised Posadas and Yuto, both of which belong to the Amazon domain, and they correspond to the regions of Argentina with the highest temperature and humidity records. Cluster 2 included the remaining localities, which are part of the Chaco domain, characterized by lower rainfalls. The haplotype network derived from these analyses was fairly consistent with the nested design reported by Alberti et al. (2008) and suggested that a single migration started from the ancestral location near Posadas. We used inductive inference with a landscape genetics approach to associate the patterns of genetic variation observed in A. fraterculus populations. In empirical studies, the true driving processes are never known, only inferred. One of the greatest potential strengths of simulation modelling is in adapting the results of empirical studies to the consideration of genetic connectivity over larger spatial extents or responses to future scenarios of landscape change (Balkenhol et al., 2015). In particular, we found evidence supporting the hypotheses of temperature, humidity and vegetation status as factors that might be responsible for creating genetic variation.

Genetic, environmental and demographic data provide basic information for interpreting spatial dynamics, which could be useful for the development of population models (Soemargono et al., 2011) and the design of efficient sampling programmes for population estimation and pest management (Soemargono et al., 2011; Cladera et al., 2014). In this context, it should be noted that the definition of populations and potential management units might require not only geographical data, but also other spatial, ecological and temporal data (Ruiz-Arce et al., 2012). There is abundant literature suggesting that management should be undertaken at a broad scale, rather than on a fine scale (Aluja et al., 2012; Karsten et al., 2013); otherwise, any local effort will be ineffective, prone to foster insecticide resistance, and will render the programmes economically unsustainable (Aluja et al., 2012). In the case of the South American A. fraterculus populations analysed, all our results suggest that the identified clusters, although representing a single evolutionarily significant unit, can be treated as two separate management units. Moreover, human-mediated dispersal of A. fraterculus has been increasing with current globalization and single-crop techniques. Therefore, it might be important to screen routes of fruit movement and establish quarantine procedures to reduce the pathways of invasion and gene flow patterns (Karsten et al., 2013).

Considering the economic importance of this species, there is a need for an integrated and more precise knowledge of its intraspecific divergence over time and space. All the results presented in this study provide useful information and a first multidisciplinary approach to gain a comprehensive overview of the genetics, ecology and dispersion patterns of A. fraterculus from Argentina, and the results highlight the importance of including environmental covariates in the development of demographic and evolutionary models.

Our findings should be viewed with some caution because of two facts: they are based solely on mitochondrial DNA sequence data, dependent on the assumption of an intrinsic mutation rate and its underlying mode of inheritance and effective population size (Avise, 2000), and our sample size is relatively small. However, all the analyses were highly consistent with each other, indicating statistical robustness. Moreover, our results on the colonization route of A. fraterculus are in agreement with cytogenetic evidence (Giardini et al., 2020). Further studies integrating ecological, morphological and genomic perspectives, and a sampling over the entire distribution range of A. fraterculus are needed, and the results would help us to understand more fully the origin and evolution of the species and its geographical expansion.

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article at the publisher’s web-site:

Table S1.COII sequence names and GenBank accessions.

ACKNOWLEDGEMENTS

L.I.F. and P.G.-C. contributed equally to this work. This research was supported by funding from Universidad de Buenos Aires (UBA; UBACYT 20020170100270BA to J.C.V.) and Agencia Nacional de Promoción Científica y Tecnológica (ANPCYT; PIC 2018-02567 to Dr María Isabel Remis). We thank Dr M. Aluja and two anonymous reviewers for their helpful comments, which contributed to improve this manuscript. We are also indebted to Dr J. Cladera for the critical reading and style revision. The authors declare that there is no conflict of interest regarding the publication of this article.