What Can Integrated Analysis of Morphological and Genetic Data Still Reveal about the Anastrepha fraterculus (Diptera: Tephritidae) Cryptic Species Complex?

The South American fruit fly Anastrepha fraterculus (Wiedemann) is a complex of cryptic species, the so-called “Anastrepha fraterculus complex”, for which eight morphotypes are currently recognized. A previous analysis of ITS1 in samples of the Anastrepha fraterculus complex, while revealing high distinctiveness among samples from different localities of South America, Central America, and Mexico, no direct association was made between sequence type and morphotype. In the present report, a correlated analysis of morphometry and ITS1 data involved individuals belonging to the same population samples. Although showing a low level of intra-populational nucleotide variability, the ITS1 analysis indicated numerous inter-population sequence type variants. Morphotypes identified by morphometric analysis based on female wing shape were highly concordant with ITS1 genetic data. The correlation of genetic divergence and morphological differences among the tested samples gives strong evidence of a robust dataset, thereby indicating the existence of various taxonomic species within the A. fraterculus complex. However, the data revealed genetic and morphological variations in some regions, suggesting that further analysis is still required for some geographic regions.


Introduction
The genus Anastrepha Schiner is endemic to the Neotropics and has been found to occur from the southern United States to northern Argentina and the Caribbean Islands [1]. The genus comprises about 250 valid species with 21 species groups that differ in terms of morphology [2][3][4]. Among the 34 species of the "fraterculus group", some are considered economically damaging, such as Anastrepha fraterculus (Wiedemann). In the 1940s, variations in wing pattern along its wide geographic distribution, similar to that of the genus, lead to the assumption that cryptic species might be involved [5][6][7]. Since then, samples from the nominal Anastrepha fraterculus populations have been analyzed for other biological parameters and significant data on the existence of different species have been accumulated over time.
The first non-morphological evidence indicating that different species could be involved in the Anastrepha fraterculus complex came from a karyotypic analyses of specimens collected in Brazil and

Biological Material
Specimens of the AF complex from 23 samples were obtained from 20 localities spread throughout the Neotropics: Mexico, Guatemala, Colombia, Ecuador, Peru, Argentina, and Brazil (Table 1 and Figure 1). Nineteen of these samples were derived from field-infested fruits and brought to the laboratory and kept until adult emergence, whereas two samples belonged to laboratory colonies, and two others were derived from McPhail traps. While all samples were screened for ITS1 sequences, 19 also underwent female wing-shape morphometry analysis. For the ITS1 analysis, adult females preserved in ethanol were used, except for those collected in Brazil, which emerged directly from the collected infested fruits. Adult flies from Brazil were identified using a combined analysis of biological characteristics [19,20], whereas flies from other countries were classified by morphometry, as previously proposed [27,28].

Biological Material
Specimens of the AF complex from 23 samples were obtained from 20 localities spread throughout the Neotropics: Mexico, Guatemala, Colombia, Ecuador, Peru, Argentina, and Brazil (Table 1 and Figure 1). Nineteen of these samples were derived from field-infested fruits and brought to the laboratory and kept until adult emergence, whereas two samples belonged to laboratory colonies, and two others were derived from McPhail traps. While all samples were screened for ITS1 sequences, 19 also underwent female wing-shape morphometry analysis. For the ITS1 analysis, adult females preserved in ethanol were used, except for those collected in Brazil, which emerged directly from the collected infested fruits. Adult flies from Brazil were identified using a combined analysis of biological characteristics [19,20], whereas flies from other countries were classified by morphometry, as previously proposed [27,28].

Analysis of ITS1 Sequences
DNA for the ITS1 analysis was extracted from 10 individuals of each sample collected in Brazil and from 5 individuals from each sample for the other localities. DNA was extracted from the isolated thorax of each individual in accordance with the Jowett protocol [36]. For PCR amplification, new primers based on the ITS1 sequence of A. fraterculus (s.l.) (GenBank #AY775552) were used, i.e., 18SF-TAACTCGCATTGATTAAGTCCC and 5.8SR-GATATGCGTTCAAATGTCGATG, which bind to the 3 end of 18S and the 5 end of 5.8S genes, respectively. PCR reactions consisted of one cycle (2 min at 94 • C), 31 cycles (30 s at 94 • C, 30 s at 60 • C, and 2 min at 68 • C) and a 5 min extension at 72 • C. Amplified products were electrophoresed in 0.8% agarose gel and visualized after ethidium bromide staining. Oneshot Escherichia coli cells were employed to generate five clones from each fragment by applying Topo TA Cloning (Invitrogen, Carlsbad, CA, USA). Plasmid DNA was extracted using the Perfectprep Plasmid Mini Kit (Eppendorf, Hamburg, Germany). The BigDye reaction kit in an ABI Prism Automatic Sequencer (Applied Biosystems, Foster City, CA, USA) was employed for sequencing at the Institute of Chemistry, University of São Paulo, São Paulo, Brazil.
Amplified fragments of around 900 bp contained complete ITS1 and partial sequences of the flanking genes. After comparison with the ITS1 sequences of A. fraterculus (l.s.) (GenBank #AY775552), the regions corresponding to 18S and 5.8S flanking genes were trimmed off. Amplification and sequencing of A or T homopolymeric regions are prone to the production of artifacts that are interpretable as either intra-individual polymorphisms or PCR slippage. Analysis of the ITS1 electropherograms of five clones from each sampled specimen with the Electropherogram Quality Analysis web tool [37] indicated the absence of intra-individual polymorphisms. Since 2 to 10 individuals were analyzed per sample, sequences differing, for example, in a homopolymeric region by only one extra nucleotide that was present in every individual of a given sample or in individuals from distinct samples were interpreted as a real nucleotide substitution and not as a PCR artifact. It would be highly unlikely that PCR errors had occurred in so many individuals at the same position in the sequences. The amplified ITS1 sequences were deposited in GenBank and are listed in Table 2. Some of the samples analyzed herein (Tables 1 and 2) were collected in geographic areas that overlap with samples analyzed by Sutton et al. [32], i.e., from Mexico/Guatemala, the lowlands of Peru/Ecuador, the highlands of Colombia, Peru, Bolivia, and south-eastern Brazil. For the sake of comparisons, the ITS1 sequences of the morphotypes identified herein were named accordingly as sequence types already described-types TI/TIa, TII, TIII(A,B,C), TIV-available in the GenBank [32].
Insects 2019, 10, 408 6 of 23 Furthermore, three previously described sequences from Paraguay (HQ829865, HQ829868, HQ829876) were also analyzed [38]. The sequence set was aligned through Clustal Omega software [39], based on a previously proposed alignment [32], and genetic diversity was estimated by MEGA 6.0 [40]. Overall similarity among sequences was tested by UPGMA (unweighted pair-group method with arithmetic means) [41], the same methodology employed a priori to compare the ITS1 sequence types [32]. For the phylogenetic inferences the maximum likelihood (ML) [42] and fast minimum evolution (FME) provided at the NCBI site were also applied. The ML inference was made based on the best substitution model determined by the tool "Search the Best Model" implemented in Mega 6.0. According to the Akaike information criterion (corrected), the best model was the GTR (general time reversible) model. Regarding ML, the initial trees for the heuristic search were obtained automatically by applying the BioNJ method to a matrix of pairwise distances estimated using the maximum composite likelihood method [43] and then selecting the topology with the superior log likelihood value. The FME was based on pairwise distances inferred according to Jukes-Cantor methodology. Gamma distribution was used to model evolutionary rate differences among sites. The percentage of replicate trees in which the associated sequences clustered together was measured by bootstrap tests (1000 replicates). The analyses were conducted in MEGA 6.0.

Morphometric Assessment
Landmark-based geometric morphometry was applied to quantify and distinguish wing morphological features in 239 female specimens from 19 population samples of the AF complex (Table 1). Morphological studies were carried out using female wings since among the usual features employed in Anastrepha taxonomy (wings, ovipositor, mesonotum) it is the most appropriate to precisely recognize homologous landmarks [20,29].
For permanent preparations, the right-wings of 10-12 specimens per sample were removed from their bases, washed with distilled water, and later dehydrated in a gradual alcohol series (50%, 70%, 100%), with 20 min for each step. Then, they were placed in xylene for 2-3 min and mounted in a dorsal view between the slide and coverslip with Canada Balsam. Images of each wing were digitalized at an identical magnification with a high-resolution 5050 Olympus camera fitted to a Zeiss Stemi 508 stereomicroscope.
The variations in female wing-shape were assessed from digital images by recording 18 homologous landmarks for each specimen using the TPS software package [44][45][46]. Landmarks were defined by the intersection or termination of wing veins, as previously described [29]: (1)  A generalized Procrustes analysis using the MorphoJ program was conducted by overlapping landmark configurations that remove variations caused by differences in translation, orientation, and size into a common coordinate system [47]. A Canonical variate analysis was performed with the matrix of Procrustes coordinates for 19 population samples. A cluster analysis of phenotypic similarities using the matrix of squared Mahalanobis distances (SMD) involving pairwise comparisons by the UPGMA method (unweighted pair-group method with arithmetic means) was implemented in the Statistica software program (StatSoft, Inc., Tulsa, OK, USA, 2006). A morphotype centroid-size analysis was done with a Kruskal-Wallis test, followed by Dunn's multiple comparisons with Bonferroni correction through the PMCMR package in the program R [48]. Reference samples and permanently mounted slides were deposited at the Entomological Collection of the Instituto de Ecología A.C., Xalapa, México, and the wing preparations from the Brazilian samples were deposited in the Laboratory of Evolutionary A generalized Procrustes analysis using the MorphoJ program was conducted by overlapping landmark configurations that remove variations caused by differences in translation, orientation, and size into a common coordinate system [47]. A Canonical variate analysis was performed with the matrix of Procrustes coordinates for 19 population samples. A cluster analysis of phenotypic similarities using the matrix of squared Mahalanobis distances (SMD) involving pairwise comparisons by the UPGMA method (unweighted pair-group method with arithmetic means) was implemented in the Statistica software program (StatSoft, Inc., Tulsa, OK, USA, 2006). A morphotype centroid-size analysis was done with a Kruskal-Wallis test, followed by Dunn's multiple comparisons with Bonferroni correction through the PMCMR package in the program R [48]. Reference samples and permanently mounted slides were deposited at the Entomological Collection of the Instituto de Ecología A.C., Xalapa, México, and the wing preparations from the Brazilian samples were deposited in the Laboratory of Evolutionary Studies on Fruit Flies, Department of Genetics and Evolutionary Biology, Institute of Bioscience, University of São Paulo, São Paulo, Brazil.

Analyses of ITS1 Sequences
The 32 ITS1 sequences obtained varied from 530 to 583 nucleotide (nt) positions, were rich in AT (ca. 70-84%),and showed complex structures. As described in detail for other samples of the AF complex [32], the sequences had poly (A)/poly (T) portions, highly variable regions interspersed with conserved subsequences, and scattered SNPs (single nucleotide polymorphism). Although most samples had a single sequence type, in some, 2 to 3 sequences were detected with 1 or 2 nucleotide substitutions. Since this represents a similarity of 99% or more, it was assumed that the sequences were alike, and therefore, the sample was represented by a single sequence, i.e., the predominant one ( Table 3). Classification of the representative sequences obtained from distinct AF morphotypes compared with representative sequence types found previously [32] are described below and shown in Figure 3. However, inconsistencies were found in data related to the collection sites, sequence types, and the nucleotide sequences in the GenBank. For example, in [32] (Table 2), collections 43-44 from Colombia (San Andrés Island) are listed as having the TI/TIa sequence type, whereas in the text [32] (p. 186), they are stated as being type TII. A comparative analysis of the sequences deposited in the GenBank showed they have 99% similarity with sequence type TII, which is found in Mexico, Central America, and Venezuela. Another case involving sequences from Colombia (collections 17-20, 25, 40, 41) is listed as being type TIIIC. An analysis of sequences from GenBank indicated the presence of SNPs and one unique repeat ATT, as described in [32] (p. 188), thereby proving that they are of type TIV. Thus, the types of ITS1 sequences obtained were determined by comparisons with a set of 29 sequences available in GenBank for which the respective accession codes could be unambiguously correlated with their collection sites (localities and altitude) [32]. These sequences are listed in Table 2.

Analyses of ITS1 Sequences
The 32 ITS1 sequences obtained varied from 530 to 583 nucleotide (nt) positions, were rich in AT (ca. 70-84%),and showed complex structures. As described in detail for other samples of the AF complex [32], the sequences had poly (A)/poly (T) portions, highly variable regions interspersed with conserved subsequences, and scattered SNPs (single nucleotide polymorphism). Although most samples had a single sequence type, in some, 2 to 3 sequences were detected with 1 or 2 nucleotide substitutions. Since this represents a similarity of 99% or more, it was assumed that the sequences were alike, and therefore, the sample was represented by a single sequence, i.e., the predominant one ( Table 3). Classification of the representative sequences obtained from distinct AF morphotypes compared with representative sequence types found previously [32] are described below and shown in Figure 3. However, inconsistencies were found in data related to the collection sites, sequence types, and the nucleotide sequences in the GenBank. For example, in [32] (Table 2), collections 43-44 from Colombia (San Andrés Island) are listed as having the TI/TIa sequence type, whereas in the text [32] (p. 186), they are stated as being type TII. A comparative analysis of the sequences deposited in the GenBank showed they have 99% similarity with sequence type TII, which is found in Mexico, Central America, and Venezuela. Another case involving sequences from Colombia (collections 17-20, 25, 40, 41) is listed as being type TIIIC. An analysis of sequences from GenBank indicated the presence of SNPs and one unique repeat ATT, as described in [32] (p. 188), thereby proving that they are of type TIV. Thus, the types of ITS1 sequences obtained were determined by comparisons with a set of 29 sequences available in GenBank for which the respective accession codes could be unambiguously correlated with their collection sites (localities and altitude) [32]. These sequences are listed in Table 2.

Characterization of the ITS1 Sequence Types
The representative ITS1 sequences obtained were aligned with previously distinguished representative sequence types-types TI, TIa, TII, TIIIABC, TIV [32]-from sites overlapping those of the present analysis ( Figure 3). ITS1 sequences have an interrupted poly(A) region and bases 1-109 and 618-669 (the numbers refer to the nucleotide positions shown in the alignment of Figure 3), respectively, at their 5 and 3 ends. In the intervening regions, there are two hypervariable subregions, HVRI and HVRII, with bases 380-420 and 539-582, respectively, interspersed with conserved subsequences and scattered SNPs (single nucleotide polymorphism). However, as the alignment of the 5´Interrupted Poly (A) Region (bases 1-109) is clearly problematic due to the high variability among samples, these were removed in the characterization of the ITS1 sequence of the distinct morphotypes.

Sequence Type TI
This type, encountered in samples from south-eastern Brazil, Paraguay, Bolivia, and the highlands of Peru, has a unique TCACATATA expansion (bases 301-310), an ATT extension (bases 368-376), and SNPs at bases 353, 403, and 417. A variant of this sequence, TIa, besides the TI markers, has a unique ATT extension (bases 580-582) inside Highly Variable Region II in samples from Peru (Pe96VRreg, Pe63AyaHua) and Bolivia (BoSCr). TIa was also found in the ArHmol sample from Argentina. Other samples of Brazilian morphotypes had the basic TI markers but showed differences characterizing other TI variants. The sequences from the Brazilian-1 morphotype had characteristic CCC (bases 129-131) and CCCC (bases 177-180) expansions, as well as a few other SNPs; sequences of Brazilian-2 morphotype had an additional G at positions 125 and 232, an additional A at position 668, and other nucleotide substitutions. The sequences found in the three samples of Brazilian-3 morphotype had a unique TAT extension at positions 359-361 and other SNPs. Hence, the sequences were considered new variants of TI associated specifically with the three Brazilian morphotypes and were named TIb (Brazilian-1), TIc (Brazilian-2), and TId (Brazilian-3). The sequence from Tucumán (ArTuc, Argentina) was considered a variant of TI but was distinct from the other variants (TIa to TId) due to its unique extensions, (T 9 (bases 291-298), C 4 (bases 302-305), CACA (bases 493-496), and GGAAA (bases 602-606)), as well as several SNPs, and as such, it was named TIe. Genetic distances between these sequence types were low, the highest being those of ArHmol (mean of 0.014) and ArTuc (mean of 0.142) relative to the other samples (Table S1). The average intragroup genetic distance between sequences was 0.039 (Table S2). Although the Mexican morphotype is characterized by sequence type TII, the high level of divergence in certain populations, due to the presence of numerous unique markers indicates the presence of variants, TIIa (MxJic), TIIb (MxQro), and TIIc (MxTap).

Sequence Types TIII A, B, C
Although mutually very similar, TIII A was characterized by an expansion in the HVRI CAATATATA (bases 380-388) by insertions of an A (bases 391 and 569) and an AT (positions 419-420) and possible A expansions at base 605. This sequence type was found in the lowlands of Peru and Ecuador (<900 m) but showed little divergence between the sequences PeCaj, PePiu, and PeLmo. The average distance between TIIIA sequences was 0.002 (Table S3); consequently, the sequence type TIIIA characterized the Peruvian morphotype.

Sequence Type TIV
This type has unique SNPs at positions 119, 151, 343, and 421 and a unique repeat expansion (ATT) n (n = 2 and 3) at positions 539-541. Every sequence from the highlands of Colombia (>1500 m) was type IV, and the average distance between samples was 0.006 (Table S4). The sample from Ibagué (CoIba) obtained differed by possessing several additional SNPs other than the four common ones of sequence type TIV, as well as the addition of G, T, or A at bases 125, 131, 137,139,162, 166, 232, 243, and 445. Furthermore, the CoIba sequence also showed the unique repeat expansion (ATT) 3 characteristic of sequence type TIV. The data indicate that Andean morphotypes are characterized by having sequence type TIV, with at least one variant, CoIba, which could be named TIVa.

Relationships of Sequence Types among Samples of the AF Complex
To achieve an improved comparison between the sequences obtained with those previously described [32], the similarities among the entire set of 55 ITS1 sequences (Table 2) were inferred by a UPGMA methodology. The results are shown in Figure 4. UPGMA revealed four main clusters and placed four sequences as single terminal branches: one from Argentina (ArTuc), one from Colombia (CoIba), and three from Mexico (MxTap, MxQro, MxJic). Each of the main cluster aggregates' specific sequence types corresponded to a distinct morphotype of the AF complex: One cluster with sequence types TIV and TIVa corresponded to the Andean morphotype, another with the TIIIA sequence corresponded to the Peruvian morphotype, and one cluster with sequence type TII was found in the Mexican morphotype, as well as a large cluster with sequence type TI and its variants. The latter was found in samples from the highlands of Peru and Bolivia as well as Paraguay, Argentina, and south and south-eastern Brazil. This cluster was subdivided into a branch with the ArHmol sequence from Argentina (type TIa), a sub-cluster from Brazil with Brazilian-1 morphotype (TIb), a branch from Peru with the sequence Pe64AyaPon (type TI), another sub-cluster from Brazil comprising morphotype Brazilian-2 with its sequence (TIc), and a larger sub-cluster comprising sequence types TI and TIa from Peru, Bolivia, and Paraguay, sequences of morphotype Brazilian-3 from Brazil (TId), and sequences from Argentina. The sequence Br09ParCha (type TI) included in this sub-cluster is the only one available in GenBank among all the studied samples from Brazil [32].
Additionally, the genetic relationships among the sequences were inferred by a maximum likelihood phylogenetic analysis by employing the entire set of ITS1 sequences. The ML tree showed that the sequences were also grouped in four clades, each of which congregated the same set of sequences as those included in the clusters of the UPGMA analysis ( Figure 5). The clade with samples from the Mexican morphotype was the most basal. There was also a clade with sequences from the Peruvian morphotype, one with sequences of the Andean morphotype, and the last one with sequences of the three Brazilian morphotypes and the sequences TI/TIa from Peru, Bolivia, Paraguay, and Argentina. Within this clade, the sequences of the Brazilian-1 and Brazilian-2 samples were separated into distinct, although internal, sub-clades, while the sequences of Brazilian-3 and all the other sequences were in another sub-clade. The ArTuc (Argentina) sequence and two sequences from Mexico (MxTap, MxJIc) were isolated in basal terminal branches, similar to the UPGMA dendrogram, but the sequences MxQro and CoIba isolated as basal terminal branches in the UPGMA were included within the clades of the Mexican and Andean morphotypes, respectively, in the ML analysis. The second phylogenetic inference was made by the fast minimum evolution methodology, and the topology was very similar to that found by the ML inference. The main distinction is that the three Brazilian sub-clades were more distinctly separated from each other and from the other sequences, and two of the Mexican morphotypes (MxTap; MxJic) were still isolated basal terminal branches ( Figure S1).
A comparison of morphological distances among the three Brazilian groups exhibited high contrast, with the greatest similarity reported between Brazilian-1/Brazilian-3 (SMD = 7.17), whereas for Brazilian-1/Brazilian-2 (SMD = 45.9) and Brazilian-2/Brazilian-3 (SMD = 49.2), the values were up to six times lower. In other geographical regions in the Americas, a close similarity was also noted for the wing-shape between Peruvian/Mexican (SMD = 9.18), as compared with the low similarity Figure 5. Phylogenetic relationships of the ITS1 sequences of the AF complex. The evolutionary history was inferred by using the maximum likelihood method based on the general time reversible model. The tree with the highest log likelihood is shown. Initial tree(s) for the heuristic search were obtained by applying the BioNJ method to a matrix of pairwise distances estimated using the maximum composite likelihood (MCL) approach. A discrete gamma distribution was used to model evolutionary rate differences among sites. The numbers on the branches are bootstrap values >50% (1000 replications). All positions containing gaps and missing data were eliminated.
A comparison of morphological distances among the three Brazilian groups exhibited high contrast, with the greatest similarity reported between Brazilian-1/Brazilian-3 (SMD = 7.17), whereas for Brazilian-1/Brazilian-2 (SMD = 45.9) and Brazilian-2/Brazilian-3 (SMD = 49.2), the values were up to six times lower. In other geographical regions in the Americas, a close similarity was also noted for the wing-shape between Peruvian/Mexican (SMD = 9.18), as compared with the low similarity between Mexican/Andean morphotypes (SMD = 55). Meanwhile, minor similarities in wing shape were found between the Brazilian-3/Andean (SMD = 70.4) and Brazilian-3/Peruvian morphotypes (SMD = 74.2) ( Figure 6).  Partial-warp scores captured the shape variation to visualize the major trends associated with hypothetical changes found in the samples. These shifts in wing shape are depicted by an outline drawing and plot of the two first canonical variates by morphotype (Figure 7).
The centroid size used to measure geometrical dimensions was compared among morphotypes resulting in a significant differentiation of groups (Kruskal-Wallis: χ 2 = 50.93, p < 0.0001), showing that the Mexican and Peruvian groups possess smaller wing sizes than all the other morphotypes (Dunn test in all p < 0.010 comparisons) ( Figure 8). Furthermore, the model designed to distinguish among morphotypes was confirmed by the classification of the 239 individuals examined, as 99.2% were classified properly into expected clusters, proving the high reliability of the predictive model to distinguish them based on natural groups (Table 4).  Partial-warp scores captured the shape variation to visualize the major trends associated with hypothetical changes found in the samples. These shifts in wing shape are depicted by an outline drawing and plot of the two first canonical variates by morphotype (Figure 7).
between Mexican/Andean morphotypes (SMD = 55). Meanwhile, minor similarities in wing shape were found between the Brazilian-3/Andean (SMD = 70.4) and Brazilian-3/Peruvian morphotypes (SMD = 74.2) (Figure 6). Partial-warp scores captured the shape variation to visualize the major trends associated with hypothetical changes found in the samples. These shifts in wing shape are depicted by an outline drawing and plot of the two first canonical variates by morphotype (Figure 7).
The centroid size used to measure geometrical dimensions was compared among morphotypes resulting in a significant differentiation of groups (Kruskal-Wallis: χ 2 = 50.93, p < 0.0001), showing that the Mexican and Peruvian groups possess smaller wing sizes than all the other morphotypes (Dunn test in all p < 0.010 comparisons) (Figure 8). Furthermore, the model designed to distinguish among morphotypes was confirmed by the classification of the 239 individuals examined, as 99.2% were classified properly into expected clusters, proving the high reliability of the predictive model to distinguish them based on natural groups (Table 4).  The centroid size used to measure geometrical dimensions was compared among morphotypes resulting in a significant differentiation of groups (Kruskal-Wallis: χ 2 = 50.93, p < 0.0001), showing that the Mexican and Peruvian groups possess smaller wing sizes than all the other morphotypes (Dunn test in all p < 0.010 comparisons) (Figure 8). Furthermore, the model designed to distinguish among morphotypes was confirmed by the classification of the 239 individuals examined, as 99.2% were classified properly into expected clusters, proving the high reliability of the predictive model to distinguish them based on natural groups (Table 4).

Wings and ITS1 Data Relationships
A comparison between genetic and morphological data was conducted on the 19 population samples analyzed by the two methods ( Figure 9A, B).
Clustering analyses revealed an aggregation into two widely discrete sets of samples, the first one consisting of four Argentinian populations (ArCon, ArHmo, ArMis, ArTuc), along with two Brazilian ones (B1Bot, B1Ube)-all of them recognized as the Brazilian-1 morphotype [28] (sensu Hernández-Ortiz et al.)-and another one including samples B3Isa and B3Uba, which belong to the Brazilian-3 morphotype [19] (sensu Selivon et al.). Although somewhat divergent in topology, the clustering of ITS1 sequences produced a similar grouping, which separated the samples of Brazilian-1, Brazilian-3, and Argentina from all the other remaining samples. The samples from Argentina exhibited higher similarity with those of Brazilian-3 than those of Brazilian-1.
The second group consisted of all remaining populations clustered into four clusters, each of them identifying a distinct morphotype: 1) Andean, represented by a Colombian highland sample (CoIba) showing the largest divergence in wing shape, even among the nearest geographical samples, e.g., Peru and Ecuador; 2) Brazilian-2, depicted separately by only one Brazilian population (B2Isa); 3) Mexican, including all samples from Mexico and Guatemala (GuCit, MxQro, MxApa, MxJic, MxTeo, MxTap); and 4) Peruvian, composed of the Peruvian and Ecuadorian Pacific lowland samples (PeMol, PePiu, EcGu4). The clustering of ITS1 sequences from these samples yielded a highly similar dendrogram that separates the Brazilian-2 (Br2Isa) morphotype sequence from two other clades: one grouping the Mexican morphotype sequences and another one consisting of the three Peruvian

Wings and ITS1 Data Relationships
A comparison between genetic and morphological data was conducted on the 19 population samples analyzed by the two methods ( Figure 9).
Clustering analyses revealed an aggregation into two widely discrete sets of samples, the first one consisting of four Argentinian populations (ArCon, ArHmo, ArMis, ArTuc), along with two Brazilian ones (B1Bot, B1Ube)-all of them recognized as the Brazilian-1 morphotype [28] (sensu Hernández-Ortiz et al.)-and another one including samples B3Isa and B3Uba, which belong to the Brazilian-3 morphotype [19] (sensu Selivon et al.). Although somewhat divergent in topology, the clustering of ITS1 sequences produced a similar grouping, which separated the samples of Brazilian-1, Brazilian-3, and Argentina from all the other remaining samples. The samples from Argentina exhibited higher similarity with those of Brazilian-3 than those of Brazilian-1.
Insects 2019, 10, x FOR PEER REVIEW 16 of 23 morphotype samples, along with the Colombian (CoIba) sample from the Andean morphotype, which was not differentiated, as shown in the morphotype dendrogram.

Discussion
From a morphological point of view, the AF complex comprises at least eight morphotypes. The genetic divergence correlates with morphological differences for six of these morphotypes, as described in the present analysis, strongly confirming the existence of various taxonomic species within the AF complex. A previous analysis of ITS1 in samples of the AF complex revealed high distinctiveness throughout the Neotropical Region [31,32,38]. However, the lack of morphological characterization in these studies precluded the association of the ITS1 sequence types with the various morphotypes. This report presents results from the correlated analysis of morphometry and ITS1 data in specimens from the same populations.

ITS1 as a Molecular Marker
This molecular marker has been used in the genetic analysis of very closely related taxa for which a recent evolutionary history is assumed [49]. The multiple copies of the ITS1 present in a genome might complicate a genetic analysis if differences do exist among them, but the cloning of ITS1 fragments in the present analysis detected no intra-individual variations. As the intra-genome similarity of ITS1 copies has been reported for several organisms, molecular-drive mechanisms are assumed to explain the observed intra-individual homogeneities [50,51].
Differences in the fragment size were not detected in a previous ITS1 analysis of Argentinian samples [52]. If only the size of the amplified fragment had been considered herein, we would also have found no relevant differences between the samples, but the scenario is quite different when the sequences of such fragments are considered.
According to the data, the overall nucleotide variability of ITS1 sequences is generally low, varying from 0.006 (Brazilian-3) to 0.039 (Mexican). This is surprising, since non-coding regions have a tendency to accumulate mutations faster than coding sequences [53,54]. Regardless of the molecular marker employed, low genetic variability in Anastrepha samples has often been observed. Some hypotheses have been proposed to explain this phenomenon although the underlying cause(s) remain unknown [11,12,[55][56][57][58]. The low genetic distance estimated through ITS1 sequences between AF complex morphotypes could be attributed to a possible recent evolutionary divergence of the complex [11,20]. Nevertheless, this does not obscure distinction among the morphotypes, since peculiarities in ITS1 sequence structures in each of the six studied morphotypes really do exist. ). The clustering of ITS1 sequences from these samples yielded a highly similar dendrogram that separates the Brazilian-2 (Br2Isa) morphotype sequence from two other clades: one grouping the Mexican morphotype sequences and another one consisting of the three Peruvian morphotype samples, along with the Colombian (CoIba) sample from the Andean morphotype, which was not differentiated, as shown in the morphotype dendrogram.

Discussion
From a morphological point of view, the AF complex comprises at least eight morphotypes. The genetic divergence correlates with morphological differences for six of these morphotypes, as described in the present analysis, strongly confirming the existence of various taxonomic species within the AF complex. A previous analysis of ITS1 in samples of the AF complex revealed high distinctiveness throughout the Neotropical Region [31,32,38]. However, the lack of morphological characterization in these studies precluded the association of the ITS1 sequence types with the various morphotypes. This report presents results from the correlated analysis of morphometry and ITS1 data in specimens from the same populations.

ITS1 as a Molecular Marker
This molecular marker has been used in the genetic analysis of very closely related taxa for which a recent evolutionary history is assumed [49]. The multiple copies of the ITS1 present in a genome might complicate a genetic analysis if differences do exist among them, but the cloning of ITS1 fragments in the present analysis detected no intra-individual variations. As the intra-genome similarity of ITS1 copies has been reported for several organisms, molecular-drive mechanisms are assumed to explain the observed intra-individual homogeneities [50,51].
Differences in the fragment size were not detected in a previous ITS1 analysis of Argentinian samples [52]. If only the size of the amplified fragment had been considered herein, we would also have found no relevant differences between the samples, but the scenario is quite different when the sequences of such fragments are considered.
According to the data, the overall nucleotide variability of ITS1 sequences is generally low, varying from 0.006 (Brazilian-3) to 0.039 (Mexican). This is surprising, since non-coding regions have a tendency to accumulate mutations faster than coding sequences [53,54]. Regardless of the molecular marker employed, low genetic variability in Anastrepha samples has often been observed. Some hypotheses have been proposed to explain this phenomenon although the underlying cause(s) remain unknown [11,12,[55][56][57][58]. The low genetic distance estimated through ITS1 sequences between AF complex morphotypes could be attributed to a possible recent evolutionary divergence of the complex [11,20]. Nevertheless, this does not obscure distinction among the morphotypes, since peculiarities in ITS1 sequence structures in each of the six studied morphotypes really do exist.

ITS1 and Morphotype Relationships
Sequence type TI and the variant TIa, detected in an extensive area throughout south/south-eastern Brazil, Argentina, Bolivia, the highlands of Peru (Cusco region), and in Paraguay, were putatively identifiable as Anastrepha sp.1 aff. fraterculus (Brazilian-1 morphotype) [32]. However, the present data showed quite different results. Although all the samples from the highlands of Peru and Bolivia do have sequence types TI or TIa, those from Paraguay were found to be closely related to Brazilian-3, hence corroborating previous data [38]. Moreover, the Brazilian morphotypes [32] collected from the inland plateau of south/south-eastern Brazil presented sequence types that were similar, but not identical to, TI or TIa. Notably, no shared ITS1 sequences were detected in sympatric zones, e.g., the valley of the "Paraiba do Sul" river, where the three Brazilian morphotypes coexist (samples Br1Isa, Br2Isa, Br3Isa), nor in the coastal region of south-eastern Brazil, where A. sp.2 and A. sp.3 coexist (samples Br2Uba, Br3Uba) [19,20,34]. Distinct sequences characterize the other morphotypes: the Mexican morphotype is characterized by sequence type TII but with variants in three samples from Mexico; the Peruvian morphotype (lowlands of Peru and Ecuador) by sequence type TIIIA; and the Andean morphotype (highlands of Colombia) by sequence type TIV and a variant TIVa. Figure 10 summarizes information on the geographical occurrence of each morphotype associated with ITS1 sequences from Sutton et al. [32] and the present data. The specificity of the ITS1 sequences in each of these morphotypes is not surprising when considering the geographic distances between their occurrence zones. Although contact between the different morphotypes could have happened through anthropogenic actions, an important mechanism of fruit fly dispersion [15], the present results do not provide evidence of this.
In Figure 10, symbols of different colors indicate the ITS1 sequence types found for each morphotype. The circles indicate the ITS1 sequence types described by Sutton et al. [32] but with putative identification of morphotypes based on data from the present study. The inset shows the sympatric zones in south-eastern Brazil where the Brazilian morphotypes coexist.
Morphotype distinctiveness was in line with genetic differences revealed by the ITS-1 analysis, which, in turn, complied with the mutual reproductive incompatibility encountered in experimental crosses carried out with almost all of the morphotypes [20,30]. Variations observed in cytogenetic and genetic analyses of some populations from Argentina were interpreted as polymorphisms, and mating compatibility between some populations led to the assumption of the existence of a single regional taxon [59][60][61][62]. However, reproductive isolation, or conversely, compatibility, may not be the sole criterium to define species boundaries [63]. The present ITS1 and morphometric data, as well as previous linear morphometric results [27,28], revealed differences in Argentinian samples of such a magnitude that the presence of distinct entities may not be ruled out.
TIVa. Figure 10 summarizes information on the geographical occurrence of each morphotype associated with ITS1 sequences from Sutton et al. [32] and the present data. The specificity of the ITS1 sequences in each of these morphotypes is not surprising when considering the geographic distances between their occurrence zones. Although contact between the different morphotypes could have happened through anthropogenic actions, an important mechanism of fruit fly dispersion [15], the present results do not provide evidence of this. Figure 10. Approximate geographical location of Anastrepha fraterculus morphotypes and related ITS1 sequence types identified herein.
In Figure 10, symbols of different colors indicate the ITS1 sequence types found for each morphotype. The circles indicate the ITS1 sequence types described by Sutton et al. [32] but with putative identification of morphotypes based on data from the present study. The inset shows the sympatric zones in south-eastern Brazil where the Brazilian morphotypes coexist. Morphological variation among samples from Mexico has already been described [27][28][29]. The present results additionally confirm these morphological variations and reveal a partial congruence with genetic data found in the ITS1 sequences. The data suggest that further analysis is necessary to clarify the real meaning of such differences.
The existence of various biological entities within this species complex may be related to the wide continental distribution encompassing the diverse ecological conditions and geographic areas where they occur. Supposedly, the geographic distance could prevail in the maintenance of divergence. Nevertheless, reproductive isolation laboratory tests between distant geographic samples [64][65][66][67] suggested that the isolation between morphotypes has been established in the past with enough time for the appearance of reproductive barriers.
Since the Brazilian morphotypes, as a whole, have been well characterized through integrated analyses of several biological parameters, including reproductive isolation assays, they have been described as different species in previous studies (A. sp. 1, A. sp. 2 and A. sp. 3) [19,20,67]. A genetic distance of 0.128 was estimated by the isozymic analysis of 19 loci between Brazilian morphotypes. In contrast, a detailed analysis of 20 gene sequences, all subjected to strong purifying selection, registered the absence of genetic divergence among 16 samples of the nominal Anastrepha fraterculus from Brazil [68]. This is remarkable, since more than one species of the AF complex occurs in some of the sampled localities [19,20], and at least four isozyme diagnostic loci exist between A. sp.1 and A.sp.2 [20]. Thus, the absence of genetic divergence observed in [68] aligns with the low genetic divergence (0.004) between the ITS1 of the Brazilian morphotypes, which could contradict the recognition of morphotypes as units with distinct evolutionary pathways. However, as the authors themselves pointed out [68], their analysis involved limited parts of the genome, thereby limiting the detection of differentiation among the samples.

Scenarios of Divergences Within the AF Complex
No areas of sympatry have been described for the AF complex, except for the Brazilian morphotypes, and the knowledge of their host fruits is scarce. Crosses between specimens from south Brazil (A. sp.1 according to the authors) and flies from Argentina were found to be fertile [69], but it must be noted that in south Brazil, A. sp.3 is found in sympatry with A. sp.1 [19,70,71]. Thus, based on the existing information, it is reasonable to assume that divergence has been established in allopatric scenarios. It must be noted that in the present analysis, no evidence of hybridization was found among the Brazilian morphotypes, even among samples collected in zones of sympatry, a prerequisite for the occurrence of hybridization. This is contrary to the presumed evidence of hybridization inferred through an analysis of a sample from a laboratory colony of A. fraterculus from south-eastern Brazil, which led to the proposition that homoploid hybrid speciation would be a prevalent process [72].
Considering the simplest model of ecological speciation, the occupation of different ecological niches or different geographic environments could result in reproductive isolation [73]. It is recognized that allopatry and secondary gene flow can be determinants in the speciation processes, as was hypothesized for host plant shifts and host race formation in Rhagoletis. It was postulated that initial allopatry would furnish the genetic architecture conducive to a sympatric host shift for Rhagoletis pomonella divergence [74]. One may assume that such processes are involved in the evolutionary history of the Brazilian morphotypes, which are found in both allopatry and sympatry, presenting differences in the host fruit usage and reproductive incompatibilities ensured by pre-and post-mating barriers [19,20,67]. The lack of more detailed information on the different morphotypes of the AF complex, such as geographic distribution and host fruit usage, make the reconstruction of its evolutionary history premature. Still, it is reasonable to suppose a dynamic scenario in which allopatry and the exploration of different ecological niches act as driving forces in the divergence process. A dynamical analysis based in "species network" performed by generalized logical formalism methodology involving the three Brazilian morphotypes, their host fruits, and geographic distribution, revealed an unidirectional incompatibility between A. sp.1 and A. sp.2, in addition to an environmental factor related to altitude, that are involved in the distribution pattern currently observed [34]. The detailed knowledge of factors involved in the relationships between species as obtained by the dynamic analysis associated with results of integrative studies would be critical to reconstruct the scenario of divergence within the AF complex.

Conclusions
The present report describes the results of an integrated genetic and morphologic comprehensive analysis of individuals belonging to the same population samples of the AF complex morphotypes. The morphometrical analysis clearly identified six morphotypes in the samples out of the eight described for the AF complex: Brazilian-1, Brazilian-2, Brazilian-3, Peruvian, Andean, and Mexican. The analysis of ITS-1 detected specific sequences for each morphotype, thus showing a strong congruence between genetic and morphological data. Moreover, the morphological and genetic variations revealed in the Argentinian and Mexican samples indicate that further research should be done to better understand its true meaning. Overall, the data clearly confirm that morphotypes are distinct biological species bearing genetic differences, and potential tools for recognizing the taxa of the AF complex are described, which may be useful for academic and applied purposes.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-4450/10/11/408/s1, Figure S1: Fast Minimum Evolutionary inference of the ITS1 sequences of the AF complex; Table S1: Genetic distance among ITS1 sequences from the Brazilian morphotypes; Table S2: Genetic distance among ITS1 sequences from the Mexican morphotype; Table S3: Genetic distance among the ITS1 sequences from the Peruvian morphotype; Table S4: Genetic distance among ITS1 sequences from the Andean morphotype.