Introduction

The Mexican transition zone is the area in which the Neotropical and Nearctic biotas overlap (^{Halffter, 1976}). In a strict sense, this zone includes the highlands of Mexico and Guatemala, whereas northern Mexico and southern United States are clearly Nearctic, and the lowlands of southern Mexico and Central America are clearly Neotropical (^{Morrone, 2015}). Few assemblages of interactions among insects and their host plants have been documented in this Nearctic-Neotropical transition zone, even though their study offers a unique opportunity to link coevolving interactions with different evolutionary affinities.

The interaction between Trichobaris (Curculionidae) and their host plants (Solanum L. and Datura L.) occurs in this diverse and complex Nearctic-Neotropical transition zone. The weevils of the genus Trichobaris have a Nearctic distribution whereas South America is the ancestral area of origin of the Solanaceae family (^{Dupin et al., 2016}), and the region in central Mexico is the area with the greatest species richness of Datura and is thought to be the diversification center (^{Luna-Cavazos & Bye, 2011}).

From the 12 species of Trichobaris, only Trichobaris compacta (Casey, 1892) uses several Datura species as a food. Weevils oviposit into immature fruit, and then larvae eat developing seeds and associated tissues, whereas adults eat leaf tissue. Adults of Trichobaris compacta mate and rest on the host plant (^{Barber, 1935} and personal observations). Some studies show that this species is unable to develop in other solanaceous species, such as Nicotiana spp. (^{Diezel et al., 2011}). Thus, this weevil is restricted to the Datura plant species. Trichobaris compacta occurs from southeastern USA to northern and central Mexico (^{Barber, 1935}). The geographic distribution of T. compacta (^{Barber, 1935}) differs from that of Datura species that uses as host plants. Datura wrightii is distributed in the deserts of USA and Mexico (Sonoran Desert), whereas D. discolor inhabits the Mexican Pacific Coast and Balsas basin. Datura inoxia occurs in the highlands in central Mexico and USA (Chihuahuan Desert), whereas D. reburra and D. pruinosa are only present in small areas in Mexico (northwest Mexico: states of Sinaloa and Sonora; southern Mexico: states of Morelos, Guerrero and Oaxaca, respectively) (^{Luna-Cavazos & Bye, 2011}). ^{Barber (1935)} noted that T. compacta uses several host plant species distributed across a wide geographic area thus possibly promoting the formation of “races” or genetically distinctive populations associated with a particular plant species. However, we have not found evidence of lineages associated to a particular host plant species of Datura (^{De la Mora et al., 2018}).

Here, we aim to elucidate how the genetic variation of Trichobaris compacta is geographically distributed. To this end, we analyzed a fragment of COI gene in individuals from different populations. This genetic marker has been used to test the presence of major genetic discontinuities causing population structure, most likely due to extrinsic barriers (^{Avise, 2004}). Here, we show the (1) genetic diversity of 29 populations of T. compacta, population differentiation and structure estimates, (2) haplotype network, (3) the phylogeny of T. compacta haplotypes, and (4) historical demography.

Materials and methods

Trichobaris compacta occurs in southwestern USA (^{Barber, 1935}) where it feeds upon Datura species. We collected Trichobaris compacta in populations of Datura plants along the Pacific coast of Mexico. For this study, T. compacta was collected on D. wrightii, D. inoxia, D. pruinosa, D. reburra and D. discolor. We selected localities in order to cover the entire distribution of T. compacta (based on ^{Barber, 1935}). Collected adult weevils were preserved in absolute alcohol for posterior analysis. Several fruits per plant were collected. To avoid sampling relatives inside a given fruit, we performed genetic analysis on only 1 weevil from each fruit. A total of 29 locations were sampled along the host plant distribution (Table 1; Fig. 1).

Figure 1 Sampled locations of populations of Trichobaris compacta collected upon several Datura plant species (indicated on the figure), different color dots represent localities, pie charts depict the frequency of each haplotype at given locality, and the size of the circle is proportional to the sample size; each COI haplotype is indicated by a corresponding number (see Table 1 for locality codes and sample size). 

Each insect was frozen in liquid nitrogen (-196 °C), then macerated in an eppendorf tube using a micropestle. Genomic DNA was extracted using the DNeasy Tissue Kit (QiagenTM) according to the manufacturer’s protocol for animal tissue. We then amplified a 663-bp nucleotide protein coding region of the mitochondrial gene cytochrome c oxidase I (COI) by polymerase chain reaction (PCR), using the following primers: COA3107 (5´-TCT ATT ARD GGD GAD GCD CTA TCT TG-3´) and COS2183N (5´-CAR CAY YTA TTY TGR TTY TGR TTY TTY GG-3´) (^{Sota et al., 2004}).

Each PCR (total volume of 25.5 ml) contained 1 ml of DNA (20 nM), 1 ml of each primer (10 mM), 0.2 ml of Taq polymerase (1U of Go Taq (Promega), 5 ml of Buffer 5x, 0.5 ml of each nucleotide (10 mM), 3 ml of MgCl2 (25 mM) and 12.3 ml of ddH₂O. The thermal cycling conditions were as follows: an initial period of 5 min at 95 ºC, followed by 35 cycles of 60 s at 95 °C, 1.2 min at 55 °C, 60 s at 72 °C, and as a final extension for 7 min at 72 °C. PCR products were sequenced at Washington University using an ABI 3730xl sequencer (Applied BiosystemsTM). All nucleotide sequences obtained were compared, edited manually with SequencherTM 4.7, and aligned with MAFFT (^{Katoh & Standley, 2013}).

To describe the neutral variation of T. compacta, we estimated the total number of haplotypes, mutations, segregating sites (S), nucleotide diversity π, haplotype diversity (h) and the theta index of genetic diversity ϴ using the ^{Nei’s (1987)} equations implemented in DNAsp 5.1 (^{Rozas et al., 2003}).

In order to compare the levels of genetic diversity among populations of T. compacta, we estimated the fixation index (F _ST ) between populations using Arlequin v. 3.11 (^{Excoffier et al., 2005}). To assess potential genetic clustering due to geographic distribution, we analyzed population structure using 2 methods. First, we used a Bayesian clustering analysis in BAPS for linked loci (^{Corander & Tang, 2007}; ^{Corander et al., 2008}); we determined the number of genetic groups with the maximum likelihood value. We performed a mixture analysis with K values from 2 to 12, with 3 iterations each. Then, we used the output file to perform the admixture analysis. The minimum size of populations to be taken into account was set to 1, while 100 iterations were applied to estimate the admixture coefficients of individuals, 10 reference individuals from each population were used, and finally, 100 iterations were performed to estimate the admixture coefficients of the reference individuals. Second, we used a SAMOVA approach (spatial analysis of molecular variance; ^{Dupanloup et al., 2002}), which defines populations that are geographically homogeneous and maximally differentiated from each other. We repeated the analyses with 2 (the minimum for SAMOVA) to 20 numbers of groups (K) until the F _CT value (the proportion of genetic variation among groups) reached a plateau. Each run consisted of 100 replications.

The method is based on a simulated annealing procedure that aims to maximize the proportion of total genetic variance due to differences among groups of populations (F _CT ). The most significant statistical output indicates the number of groups. Simulations have shown that the SAMOVA algorithm does indeed find maximally differentiated groups, especially when data are derived from a single locus.

To visualize the frequency and distribution of haplotypes, we constructed an un-rooted haplotype network using the Median-joining algorithm. This analysis was performed in Network v. 4.6.1.1. (^{Bandelt et al., 1999}). A haplotype phylogeny was estimated with a Bayesian inference analysis in BEAST v.1.4.7 (^{Drummond & Rambaut, 2007}). The parameters of MCMC runs were 30 million generations, sampling every 1,000 generations and discarding the first 10,000 trees as burn-in. We used a GTR+I model of nucleotide substitution and a fixed substitution rate of 2% per million years for COI, as previously reported for coleopterans (^{Nakamine & Takeda, 2008}). To identify ancestral and derived haplotypes in the phylogeny, the trees were rooted using sequences of Trichobaris soror. We analyzed the sets of trees using DensiTree, showing only consensus trees (BEAST v.1.4.7 package software; ^{Drummond & Rambaut, 2007}).

To determine whether demographic events have shaped the genetic variation of T. compacta, we estimated ^{Fu and Li’s F (1993)}) and Tajima’s D statistics (^{Tajima, 1989}). These estimators assume neutrality to infer population historical demography. The Tajima’s D statistic is expected to be negative when the genetic structure has been influenced by a rapid range expansion, positive when the population has passed through a bottleneck, and zero when there is an equilibrium between mutation and drift (Tajima, 1989). To compare the observed frequencies of pairwise differences with those expected under a model of demographic expansion, mismatch distributions were generated using DNAsp v.4.10 (^{Rozas et al., 2003}). A multimodal distribution is expected when there are no changes affecting population size, but unimodal distributions are expected when sudden demographic expansions have occurred (^{Rogers & Harpending, 1992}).

Results

The southwest Mexico has not previously been reported as a distribution area of T. compacta (^{Barber, 1935}), and most sampling was performed in the northern part of its distribution on D. wrightii and D. discolor (Table 1, from localities 1-22). However, in this study, we report new data for T. compacta distribution and host plant species: D. pruinosa, D. reburra and D. inoxia (Table 1, from localities 23-29). These populations were sampled in southwest Mexico (Fig. 1).

We found 49 haplotypes, 33 segregation sites (S), and 56 singletons (Table 1). Sequences without missing bases were submitted to GenBank (KX359683 to KX359723; Table 2). Trichobaris compacta possess high genetic diversity (h = 0.704, π = 0.002, ϴ = 0.014) (Table 1). The most common haplotype, Co1, is present in the entire distribution range of this weevil (Fig. 1).

A significant genetic differentiation among populations was detected. F _ST values range from 0.076 to 1.0 among populations. The most differentiated populations were 14SD, 19SD and 3SD (Table 3, Fig. 1), indicating that the haplotypes of this region are more exclusive along the distribution of T. compacta. The Bayesian clustering analysis performed with BAPS showed that the maximum likelihood value is reached at K = 3 (Fig. 2): cluster 1 in green, cluster 2 in yellow and cluster 3 in blue. Much of the genetic variation is specific to each group 1 (97%), group 2 (100%), and group 3 (87%) with very little gene flow between groups (for example, between group 3 and group 1 is only 1.1%). With SAMOVA, the highest value of F _CT was found at K = 3, p-value = 0.05083 ± 0.00800, the criterion for choosing this number of groups was the peak of F _CT before the plateau phase (after K = 9 and up to K = 20 values F _CT for which the analysis was run) (Table 4).

Figure 2 A) Population structure of Trichobaris compacta and B) genetic flow among clusters. Each color represents a genetic group determined by Bayesian clustering from COI gene (663 bp). The localities are described in Table 1. Gene flow is indicated as weighted arrows, such that the weights equal relative average among of ancestry in the source cluster, among the individuals assigned to the target cluster. The localities are described in Table 1. 

The most frequent and shared haplotypes among populations were Co1 (52% of total sampling), Co26 (13%), Co33 (3%), and Co10 (3%) (Figs. 1, 3). From these, Co1 and Co10 were sampled in more than 1 host-plant species (D. wrightii, D. discolor, D. inoxia and D. pruinosa), while 7 haplotypes were found exclusively in D. discolor, 5 in D. reburra and 35 in D. wrightii.

Phylogenetic analysis of COI haplotypes did not support the formation of main clades into T. compacta (Fig. 4). Posterior probability values lower than 0.50 are not show in the phylogeny, only 3 clades with 2 or 3 haplotypes have support higher to 0.90.

Figure 3 Haplotype network of Trichobaris compacta calculated from COI gene (663 bp) by the “median joining” algorithm. Black dots represent vectors, and red numbers represent the number of mutational changes. 

Figure 4 Phylogeny of COI haplotypes (663 bp) of Trichobaris compacta inferred by Bayesian inference analysis. Blue numbers show the posterior probability for clades (only values above 0.90 pp are showed). The length of the branches represents the divergence among haplotypes. 

The haplotype network showed a star-like pattern indicative of a recent population expansion and congruent with both the mismatch distribution that displayed a unimodal pattern (Fig. 5) and with the negative values of Tajima’s D = -2.69330 (p < 0.001) and Fu and Li’s F = -5.10607 (p < 0.02).

Figure 5 Historical demography of Trichobaris compacta. Mismatch distribution calculated from 663 bp of the COI gene sequence. 

Discussion

In general, the COI sequence of weevils shows high levels of genetic variation (from 0.973 to 0.75) allowing to discern different evolutionary lineages and population structure of this species (^{Anducho-Reyes et al., 2008}; ^{Aoki et al., 2008}, ²⁰⁰⁹; ^{Kuester et al., 2012}; ^{Ruiz et al., 2010}). Here, we show that T. compacta have also higher level of genetic variation (h = 0.704) distributed mainly in the northern part of its distribution, USA, than in the south (Table 1; Fig. 1). Trichobaris soror, a weevil of the same genus but species-specific to a host plant and distributed on the highlands of Mexico, showed a similar level of genetic variation, h = 0.731 (^{De la Mora et al., 2015}). Trichobaris compacta diverged from T. soror ca. 1.75 (± 1) million years ago, with a recent diversification ca. 0.5 (± 0.25) million years ago (De la Mora et al., 2018).

Neutral theory postulates that genetic diversity will increase with larger effective population sizes due to a decreasing effect of genetic drift (^{Kimura, 1983}). According to the specialist-generalist variation hypothesis (SGVH), specialists are expected to have lower effective population sizes than generalists and consequently lower levels of genetic diversity due to the effect of drift (^{Li et al., 2014}). This pattern may result from their specialization in resource use; for instance, the habitat of specialist species may be more heterogeneous and patchier than that of generalists and, as a consequence, populations of specialist species may be less connected and more subdivided into smaller populations than those of generalist species. Therefore, we expected that the generalist T. compacta, which uses several host-plant species to feed and reproduce, would have higher levels of genetic variation than a specialist weevil. However, although Trichobaris soror feeds only on D. stramonium, our estimate of genetic diversity does not differ between species of Trichobaris. Therefore, no support for this hypothesis was found. It is likely that since the host plant of T. soror (D. stramonium) is highly abundant but with a patchy distribution along the Trans-Mexican Volcanic Belt, this has promoted high genetic variation of similar magnitude to that of T. compacta. However, the role of trophic host specialization and the distribution pattern of both species of Trichobaris on their genetic variation remains to be tested.

Phylogenetic structure. The population genetic structure of T. compacta investigated with COI does not reflect a significant genetic clustering, since the same haplotype can be found from Arizona in the USA, to Oaxaca in central Mexico. As mentioned below, a slight influence of the Colorado River basin may be disrupting the genetic flow in T. compacta. It is a bit difficult to compare the compact phylogeographic pattern of T. compacta with other insects since phylogeographic studies are scarce, in spite that other insect species distribute in the same geographical range (^{Anderson & O’Brien, 1996}). Unlike T. compacta, T. soror is distributed only along the Trans-Mexican Volcanic Belt and in the Sierra Madre del Sur (^{Barber, 1935}). The phylogeographical pattern of T. soror follows a partition due to the mountain systems, similar to other beetles distributed on the same region (i.e., Dendroctonus mexicanus; ^{Anducho-Reyes et al., 2008}).

Other taxa distributed in North American warm deserts have revealed vicariant events among their populations distributed in this area (^{Bryson et al., 2012}; ^{Leache & Mulcahy, 2007}; ^{Mantooth et al., 2013}), and the barriers to gene flow for this biota have been identified as the Central Valley, the Colorado River and floristic provinces. At the equatorial climate, along the western coast of Mexico, Neotropical biota have been greatly influenced by past events that have produced strong genetic differentiation within species (^{García-Trejo et al., 2004}; ^{Zarza et al., 2008}). We expected that at least the genetic variation of T. compacta would reflect these discontinuities. However, the genetic groups do not correspond to either the desert typical barriers or to a disjunction between the desert and equatorial climates. The Colorado River system seems to have slightly influenced the population structure of T. compacta because of the different haplotypic compositions found at each side (Fig. 2). However, the results of BAPS and SAMOVA do not support this barrier and estimations of gene flow are needed to corroborate this. Both analyzes suggest the formation of 3 groups but these groups are not equivalent. In BAPS there is a widely distributed group (cluster2-yellow) that includes several populations from the north and mainly from the south of the distribution (Fig. 2). While the other 2 form some geographical structure to the north of the T. compact distribution, SD populations in San Diego (cluster 1-green) and AZ populations in Arizona (cluster 3-blue). In SAMOVA, the K = 3 populations were grouped as: group 1 (0Az, 10SD, 11SD, 12SD, 13SD, 15SD, 1Az, 20SD, 2Az, 2SD, 314, 3Az, 405, 4Az, 516, 5SD, 6Az, 8SD, AP, OG, 129) group 2 (19SD, 3SD) and group 3 (14SD). Although the SAMOVA results, in general, were not statistically significant (Table 4), K = 3 was the peak of F _CT and was marginally significant. It is likely that T. compacta is broadly distributed over the Datura species that we report in Mexico and gene flow is maintained along their distribution; yet extensive sampling to delimit the distribution of T. compacta in all 5 Datura species is badly needed.

The distribution of haplotypes is not uniform in all Datura species (Fig. 1). The most frequent haplotype of T. compacta (Co1) was found on 4 species of Datura distributed in USA and Mexico. All other haplotypes were found in D. wrightii and D. discolor, but some low-frequency haplotypes were restricted to D. reburra, D. pruinosa and D. inoxia. The haplotype network shows few mutational steps among haplotypes, which suggests a recent divergence of the haplotypes, but none of these as previously reported is associated to particular Datura species (^{De la Mora et al., 2018}; Fig. 3).

Because plant-herbivore interactions have largely been maintained since the last glacial period, variation in mtDNA results useful to infer clades association to particular host plants (^{Aoki et al., 2011}). Like results reported in other weevil species with larval development inside the host plant tissues (^{Aoki et al., 2009}; ^{Barr et al., 2013}), we expected to find the COI gene haplotypes of T. compacta linked to different Datura species and clustered together with high support by a phylogeny, but we did not find such a pattern (Figs. 4, 5). The neutral nature of COI variation, in relation to host plant use made by weevils, does not mean that T. compacta is not adapted to Datura species. We are unable to assess adaptation to the host plant directly (Barr et al., 2013; ^{Morse & Farrell, 2005}). If weevils were adapted to host plants, analysis of other loci that code for adaptive traits to host plants and behavioral experiments of oviposition would be necessary. Rather, we determined that the relationship between T. compacta and Datura species is not reflected by clustering of mtDNA lineages linked to particular host plant species. It is possible that T. compacta had recently dispersed to Datura host species until attaining its current distribution range, as suggested by both the star-like pattern of the haplotype network and the historical-demographic analysis.

In the insect-plant interaction study, the phylogenetic study of variation in the COI gene revealed that interacting species share patterns of neutral genetic variation because they have been exposed to the same environmental and geological conditions (^{Aoki et al., 2009}, ²⁰¹¹). In this sense, this study constitutes one of the first steps to understand the coevolution of Trichobaris with Datura species.

Phylogeographic and phylogenetic analyses showed that T. compacta populations are slightly structured over a wide geographic distribution range. This genetic population structure is the result of the admixture of tree genetic groups, forming a single one genealogical lineage.