Introduction
The independent emergence of agriculture took place in widely dispersed regions of the world between 10,000 and 6000 years before present (yBP) and allowed for the domestication of important crops, such as cereals and pulses in the Fertile Crescent (Lev-Yadun et al., Reference Lev-Yadun, Gopher and Abbo2000; Brown et al., Reference Brown, Jones, Powell and Allaby2009), rice in Southeast Asia (Kealhofer & Piperno, Reference Kealhofer and Piperno1994), cereals in Northern India (Fuller et al., Reference Fuller, Harvey and Qin2007), and maize in Mesoamerica (Iriarte et al., Reference Iriarte, Holst, Marozzi, Listopad, Alonso, Rinderknecht and Montaña2004; Piperno et al., Reference Piperno, Ranere, Holst, Iriarte and Dickau2009). Insects and insect remains, especially of stored grain insects, are important sources of information for human history in a variety of contexts. They are frequently present in archeological sites and have been important in tracing past urban environments, grain origin and likely routes of grain trade, and history of storage (Solomon, Reference Solomon1965; Buckland, Reference Buckland1981; Oliveira et al., Reference Oliveira, Corrêa, de Souza, Guedes and de Oliveira2013).
Historically, several species of pest insects have been associated with stored grains and grain products. Three species of Sitophilus Schönherr, 1838 (Coleoptera: Curculionidae) are amongst the most damaging pest insects of stored cereals, leading to severe losses worldwide: the granary weevil Sitophilus granarius (L., 1758), the rice weevil S. oryzae (L., 1763), and the maize weevil S. zeamais (Motschulsky, 1855) (Longstaff, Reference Longstaff1981; Throne & Cline, Reference Throne and Cline1989; Levinson & Levinson, Reference Levinson and Levinson1994; Plarre, Reference Plarre2010). The granary weevil is a flightless and fully synanthropic species that is adapted to artificial grain stores (Plarre, Reference Plarre2010); this species is particularly abundant in archeological sites, but is absent from natural reservoirs (Solomon, Reference Solomon1965; Buckland, Reference Buckland1981; Levinson & Levinson, Reference Levinson and Levinson1994; Panagiotakopulu, Reference Panagiotakopulu2001). The rice weevil is a scarce flier and recorded outdoors, it occasionally occurs in cereals in the field; the maize weevil is a more intensely flying species and frequently infests maize in the field before harvest (Longstaff, Reference Longstaff1981; Throne & Cline, Reference Throne and Cline1989; Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013). The grain weevils (i.e., the granary, rice, and maize weevils) seem to exhibit ecological preferences involving climatic gradients: the granary weevil in temperate climates, rice weevil in more subtropical areas, and the maize weevil in warmer climates (Longstaff, Reference Longstaff1981; Throne & Cline, Reference Throne and Cline1989; Levinson & Levinson, Reference Levinson and Levinson1994; Plarre, Reference Plarre2010; Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013). Although host specificity is not strict (Haines, Reference Haines1981; Throne & Cline, Reference Throne and Cline1991) host preferences seem to distinguish the grain weevils to a certain extent. The maize weevil shows preference toward larger grains (e.g., maize), while both the rice weevil and the granary weevil favor small grains (e.g., wheat and rye) (Haines, Reference Haines1981; Throne & Cline, Reference Throne and Cline1991; Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013). The rice weevil is almost morphologically indistinguishable from the maize weevil, except for differences on both male and female genitalia, which are difficult to dissect and examine – particularly for females (Kuschel, Reference Kuschel1961; Halstead, Reference Halstead1963; Hidayat et al., Reference Hidayat, Phillips and Ffrench-Constant1996). The rice weevil and the maize weevil are regarded as ‘sibling species’; they have been considered two races of the same species in the past and later recognized as distinct species (Richards, Reference Richards1944; Kiritani, Reference Kiritani1956; Halstead, Reference Halstead1963; Hidayat et al., Reference Hidayat, Phillips and Ffrench-Constant1996).
The likely geographic origin of the genus Sitophilus is not yet known. Several species of Sitophilus feed on seeds of wild trees from forest areas around the Himalayas and the Indian subcontinent (Plarre, Reference Plarre2010). The tamarind weevil S. linearis (Herbst, 1797), for example, is native to India and has dispersed to all areas where tamarind is grown (Cotton, Reference Cotton1920). The grain weevils most likely originated from acorn-feeding individuals that became associated with natural stores or reservoirs of acorns associated with bird and rodent nests (Levinson & Levinson, Reference Levinson and Levinson1994; Plarre, Reference Plarre2010). In Asia, the onset of agriculture during the Neolithic period (about 10,000 yBP) may have provided the grain weevils with the opportunity to expand their territories and change their behavior to feed on stored grains. Weevil fossil impressions (dubiously recognized as maize weevil) from Jomon pottery in Japan (ca. 10,500 yBP) (Obata et al., Reference Obata, Manabe, Nakamura, Onishi and Senba2011) and fossils of the granary weevil from Israel and Europe (ca. 7000 yBP) (Buckland, Reference Buckland1981; Panagiotakopulu, Reference Panagiotakopulu2001) and Egypt (ca. 4300 yBP) (Solomon, Reference Solomon1965; Buckland, Reference Buckland1981; Panagiotakopulu, Reference Panagiotakopulu2001) attest the early association of these pest insects with stored cereals. The subsequent expansion of cereal crops, initially on Asia and later on toward other continents, likely led to the dispersal of grain weevils across Asia at first and subsequently across Europe, Africa, and the Americas (Plarre, Reference Plarre2010; Obata et al., Reference Obata, Manabe, Nakamura, Onishi and Senba2011).
Despite the economic, archeological, and historical importance of Sitophilus, the origin and diffusion of the grain weevils remain controversial – particularly the maize weevil. The ecological differences among the grain weevils together with the fossil records suggest two plausible hypotheses to account for the origin of the maize weevil: (a) a common, ancient origin of the maize weevil on Asia alongside the granary and rice weevils, or alternatively (b) a derived origin from a strain of the rice weevil – its ‘sibling species’ – that recently adapted to maize as the cultivation of this cereal became widespread as a result of the Columbian Exchange (Crosby, Reference Crosby2003). Moreover, the unprecedented speed with which insect pests can disperse worldwide owing to modern agricultural settings and human-mediated transport of infested grains likely allows for adaptive gene combinations (e.g., insecticide resistance) to become widespread rapidly. Thus, it is of phytosanitary concern to uncover the extent of which the maize weevil retains population structure at both the regional and continental levels.
Phylogeography can shed light onto the evolutionary history of the maize weevil, therefore disentangling both its genealogical relationship to the rice weevil and the timing of its association with maize. Herein we explored the late evolutionary history of the maize weevil using mitochondrial and nuclear gene sequences and microsatellite markers from populations sampled from around the world. We addressed the following four questions: (a) Did the maize weevil emerge recently from the rice weevil? That is, the rice weevil is the actual ancestor of the maize weevil owing to an agriculturally driven speciation event. (b) Alternatively, did both the maize weevil and the rice weevil share a most recent common ancestor that predates the onset of agriculture? That is, the two weevil species emerged from a common ancestor and their speciation was not an agriculturally driven event. (c) Which are the current genetic structure and levels of gene flow amongst maize weevil populations, considering either the regional, continental, or intercontinental levels? (d) What are the implications of these findings to justify phytosanitary concerns?
Material and methods
Sampling
Sampling of the maize weevil took place in 53 localities spread throughout 15 countries; while the rice weevil was obtained from 20 localities spread throughout 13 countries (fig. 1, tables S1 and S2 in the supplementary material). Specimens (n ≥ 20) from each location were fixed in 95% ethanol and kept at −20°C until subsequent use. Species identification was carried out inspecting the insect genitalia in the laboratory (Halstead, Reference Halstead1963); additionally, identification was confirmed using species-specific molecular tools (Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013).
Geographic distribution and sampling localities of the 53 populations of maize weevil (Sitophilus zeamais) and 20 populations of rice weevil (Sitophilus oryzae) used in this study (for details, tables S1 and S2 in the supplementary material).
DNA extraction, polymerase chain reaction (PCR) amplification, and sequencing
Total genomic DNA was extracted following Clark et al. (Reference Clark, Meinke and Foster2001). For the maize weevil, the PCR amplified fragments of both the mitochondrial cytochrome oxidase subunit I (COI) gene and the cytochrome oxidase subunit II (COII) gene following protocols described previously (Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013, Reference Corrêa, Tomé, Braga, Martins, Oliveira and Guedes2014). For the rice weevil, amplification of the COI fragment was carried out with a primer pair described previously by Corrêa et al. (Reference Corrêa, Oliveira, Braga and Guedes2013), while the amplification of the COII fragment was carried using a pair of primers designed with information from public databases (GenBank accession no. AY014881) and the Primer 3 software (Whitehead Institute for Biomedical Research, Cambridge, MA, USA): forward primer (5′–TTTCTTCAAGATAGAGCCTCACC–3′) and reverse primer (5′–GCTCCGCAAATTTCAGAACA–3′). The nuclear internal transcribed spacer (ITS) region (ITS1–5.8S gene–ITS2) from the ribosomal DNA gene (nrDNA) was amplified using the universal ITS1 primer (White et al., Reference White, Bruns, Lee, Taylor, Innis, Gelfand, Sninsky and White1990) in combination with a species-specific ITS2 reverse primer (Peng et al., Reference Peng, Lin, Chen and Wang2003): either 5′–CGATTGTACGAGACGGGCA–3′ (for the maize weevil) or 5′–CCGTTTAAACGATTTCATCC–3′ (for the rice weevil). For the amplification of the ITS region, we used the PCR conditions mentioned previously by Corrêa et al. (Reference Corrêa, Tomé, Braga, Martins, Oliveira and Guedes2014), except that we increased the elongation time to 2.0 min. Amplicons were sequenced using the DNA sequencing services of Macrogen Inc., South Korea (www.macrogen.com).
Assembly of sequence datasets
Sequences were imported into the program Sequencer ver. 4.8 (Gene Codes Corp.) for alignment and editing. For each mitochondrial gene (COI and COII), we obtained 360 sequences (maize weevil, n = 261; rice weevil, n = 99). Sequence for the ITS region, we obtained 96 sequences of the maize weevil and nine sequences for the rice weevil. After all of the sequences were aligned, their ends were trimmed to eliminate fragments that could not be obtained for all sequences. The alignments had 806, 551, and 1230 bp, for COI, COII, and ITS, respectively. Some sequences of ITS1 region of the maize weevil harbored a 3-bp indel (insertion/deletion event), which consisted of a microsatellite (ACG4 versus ACG5). Sequences obtained for this study were deposited in GenBank (COI: KJ397732–KJ397813, KX190215–KX190493; COII: KJ397814-KJ397895, KX190494–KX190772; ITS: KX190064–KX190214).
The presence of numts (nuclear paralogs of mitochondrial origin (Lopez et al., Reference Lopez, Yuhki, Masuda, Modi and O'Brien1994); amongst COI and COII sequences was inspected with the help of MEGA ver. 5 (Tamura et al., Reference Tamura, Peterson, Peterson, Stecher, Nei and Kumar2011). The following signatures of numts were searched: (a) indels that introduce frameshifts, (b) out-of-place inframe stop codons that lead to premature termination of protein translation, and (c) lack of codon position substitution bias toward the third position that lead to a higher rate of non-synonymous mutations. The presence of signatures (a) and (b) was enough to consider a given sequence as a numt; signature (c) was used to confirm the numt status of the sequence. The sole presence of signature (c) was not used to declare a numt.
Subsequently, we assembled three datasets. Dataset A (n = 360; 1357 bp) contained the sequences of both COI and COII concatenated. Dataset B (n = 96; 1230 bp) contained following regions: ITS1 (763 bp, partial); 5.8S rRNA (63 bp, complete); and ITS2 (399 bp, partial). Dataset C (n = 4; 806 bp) contained the most prevalent haplotype of COI from each of the two weevil species (GenBank accessions KX190227 and KX190400, respectively), supplemented with sequences available from the GenBank for the granary weevil (AY131101) and the tamarind weevil (AY131102).
Demographic statistics
Measures of nucleotide diversity were carried out according to the geographical origin of the specimens: Americas, Europe–Africa, and Asia–Australia. We assembled species-specific sub-datasets from Dataset A. Haplotype diversity and nucleotide diversity parameters were calculated using DnaSP ver. 5 (Librado & Rozas, Reference Librado and Rozas2009). Both Tajima's D and Fu's Fs tests of selective neutrality were performed in Arlequin 3.1. (Excoffier et al., Reference Excoffier, Laval and Schneider2006). The neutrality tests were tested for significance by generating 1000 random samples using coalescent simulations. The ‘Infer from distance matrix’ option for ‘Haplotype definition’ was activated, following the recommendation in the Arlequin manual; Fu's Fs statistics were considered as significant at 5% if P < 0.02. Analysis of molecular variance (AMOVA) was performed using the Arlequin 3.1 with parametric bootstrapping (1000 replicates) and significance at the 5% level (Excoffier et al., Reference Excoffier, Smouse and Quattro1992).
Genealogical inferences
Haplotype definition and genealogical relationships among haplotypes in datasets A and B were obtained independently, using the median joining (MJ) method (Bandelt et al., Reference Bandelt, Forster and Rohl1999) as implemented in Network 4.5.0.2 (Fluxus Technology Ltd.). We allowed network analyses to infer genetic connections, without taken into consideration that the sequences were obtained from distinct weevil species. We obtained two networks, one for each dataset. Ambiguities in each haplotype network were resolved using the criteria of coalescent theory (haplotype frequency) and population geography (geographical proximity) (Crandall & Templeton, Reference Crandall and Templeton1993).
Divergence dating
The times of divergence among species of weevil were estimated using Beast v1.7.5. (Drummond et al., Reference Drummond, Suchard, Xie and Rambaut2012). Beauty version 1.7.5 was fed with sequence information from dataset C to generate an XML input file (available upon request). The XML file implemented the following steps: the General Time Reversible model plus invariant sites (GTR + I), which the Akaike Information Criterion (Akaike, Reference Akaike1974) had selected as the best fit model amongst 24 models of molecular evolution using MrModeltest v. 2.3. (Nylander, Reference Nylander2004); the strict molecular clock assumption; the Yule speciation process model selected for tree prior; and a mean substitution rate of 0.0177 substitutions per site per million year (My) for the COI gene (Papadopoulou et al., Reference Papadopoulou, Anastasiou and Vogler2010), which correspond to 3.54% pairwise divergence per My. The XML file assumed that the mutation rate of the COI gene followed a normal distribution (mean = 0.0177; SD = 0.001), matching the 95% interval [0.0157, 0.196]. This analysis allowed us to estimate the age of the most recent common ancestor (tMRCA) for the maize weevil and the rice weevil using the highest known mutation rate for COI gene in beetles (Papadopoulou et al., Reference Papadopoulou, Anastasiou and Vogler2010). The analysis was run for 500 million generations, samples taken every 1 million generation, with three independent replications. Results of each replication were checked for sampling, mixing, and convergence of the Markov chains to a stationary distribution using Tracer v.1.5 (Drummond et al., Reference Drummond, Suchard, Xie and Rambaut2012). Subsequently, we pooled the results of the three replications using Tracer. These settings ensured that both model parameters and time estimates were sampled adequately, as the effective sample size (ESS) values were above 200 for all statistics in Tracer. Means and 95% highest posterior density (HPD) intervals were determined in Tracer.
Microsatellite-based analyses in maize weevil
Ten microsatellite markers (Barat et al., Reference Barat, Bravo, Chandra, Corrêa, Giombini, Guedes, Huailei, Lal, Liang, Matura, Mohindra, Oliveira, Patangia, Qiyong, Sah, Singh, Singh, Singh, Tosto, Tripathi and Vinson2012) were applied to 309 specimens, which had been sampled from 15 populations: California (USA), Mexico (MEX), Panama (PAN), Colombia (COL), Brazil (seven locations, BR1 to BR7), Mozambique (MOZ), China (CHI), India (IND), and Thailand (THA). Amplifications were performed in multiplex sets using two triplex sets (loci 3H6-1D10-3G7 and 1A1-1B1-1G7) and two duplex sets (3A11-3G1 and 1E1-1B10). Reactions were performed in a final volume of 13 µl containing template DNA (50 ng), 1× PCR buffer (10 mM Tris–HCl pH 8.4, 50 mM KCl, 1% Triton X-100, 1.5 mM MgCl2), MgCl2 (3.0 mM), dNTPs (50 µΜ each), forward and reverse primers (0.3 mM each), and 1 U Taq DNA polymerase (Phoneutria). The amplification conditions were 4 min at 94°C for the initial denaturation followed by 35 cycles of 30 s denaturation at 94°C, 45 s annealing at 55°C and 1 min elongation at 72°C, with a final elongation at 72°C for 20 min. The amplicons were sized on an ABI PRISM 3100xl DNA Analyser using the GeneScan 500 ROX Size Standard (Applied Biosystems, Foster City, CA, USA)
Hardy–Weinberg equilibrium and linkage disequilibrium were estimated on each of the 15 populations using the F STAT software (Goudet, Reference Goudet2002) and the frequency of null alleles was estimated using the applicative MICRO-CHECKER ver. 2.2.3. (Van Oosterhout et al., Reference Van Oosterhout, Hutchinson, Wills and Shipley2004). Allelic frequencies, private alleles (A PRIV), expected/observed heterozygozity and coefficient of fixation (F) were calculated using the software GDA 1.1 (Genetic Data Analysis) (Lewis & Zaykin, Reference Lewis and Zaykin2001). AMOVA was performed with Arlequin 3.1 (Excoffier et al., Reference Excoffier, Laval and Schneider2006) for two hierarchical levels (among and within populations) using option R ST; the parameters of molecular variance were tested using 1000 permutations. The genetic differentiation among populations was estimated using the parameters F IS, F IT, and F ST (Weir & Cockerham, Reference Weir and Cockerham1984) G ST (Nei, Reference Nei1973) and R ST (Slatkin, Reference Slatkin1995) with the software FSTAT ver. 2.9.3.2 (Goudet, Reference Goudet1995), and the significance was tested following 1000 permutations with 95% confidence interval.
The Bayesian clustering approach of Structure ver. 2.3.4 (Pritchard et al., Reference Pritchard, Stephens and Donnelly2000) inferred the number of genetic groups of the maize weevil using the Monte Carlo Markov Chain (MCMC) approach. The following parameters were enforced: allele frequencies independent; use no admixture model; no location priors were used. We set runs with a burn-in period of 250,000 steps followed by 750,000 steps, with 20 independent replications for every K. As suggested by Evanno et al. (Reference Evanno, Regnaut and Goudet2005) the K was set from 1 to 18, which is from one to the number of sampled populations (15) plus three. The best K was determined according to the ΔK method of Evanno et al. (Reference Evanno, Regnaut and Goudet2005) calculated with the program Structure Harvester (Earl & vonHoldt, Reference Earl and vonHoldt2012). We converged the data of the 20 replications in the best K with the software Clumpp (Jakobsson & Rosenberg, Reference Jakobsson and Rosenberg2007).
Results
Demographic statistics
Signatures of numts were absent both from COI and COII sequences; therefore, we kept the sequences for subsequent analyses. Measures of nucleotide diversity and neutrality test statistics were carried out on each species separately given that these are intraspecific measures.
For the maize weevil, there were 18 haplotypes among the 261 sequences of the subdataset of mitochondrial genes (COI + COII) and six haplotypes among the 96 sequences of the subdataset of ITS sequences (table 1). There were instances in which we found intraindividual polymorphism for ITS, which were restricted to a single site across the dataset. We created two alternative ITS sequences for each specimen that exhibited the ambiguity. Haplotype diversity (Hd) for sequences of mitochondrial origin reached values that were similar to those found for sequences of nuclear origin (about 0.72). The sample size was much smaller in ‘Asia-Australia’ (COI + COII, n = 31; ITS, n = 11) than in ‘Americas’ (COI + COII, n = 200; ITS, n = 73); however, both of these groups exhibited similar values of nucleotide diversity. These results indicated that the small sample size of ‘Asia–Australia’ was able to capture a relatively large amount of polymorphism. For the rice weevil, there were ten haplotypes among the 99 sequences of the mitochondrial genes (table 2); however, no polymorphisms were found among the nine sequences of the ITS region. Once again, measures of diversity in ‘Asia–Australia’ reached high values, even though this group had a small sample size.
Measures of genetic diversity and neutrality test statistics for Sitophilus zeamais, based on two concatenated mitochondrial genes (COI and COII) and the nuclear internal transcribed spacer (ITS) region.
Measures of genetic diversity and neutrality test statistics for Sitophilus oryzae, based on two concatenated mitochondrial genes (COI and COII).
Tests of selective neutrality in both weevil species showed non-significant values for Tajima's D (P > 0.05) and Fu's Fs (P > 0.05) for most geographic groups (tables 1 and 2). Such non-significant values are consistent with COI, COII, and ITS sequences harboring randomly evolving mutations (neutrally evolving DNA). The only exception was the significant negative values of Fu's Fs in ‘Asia–Australia’ of the maize weevil (table 1). Such significant negative values of Fu's Fs indicated that the populations of the maize weevil in ‘Asia–Australia’ harbored an excess of low-frequency polymorphism and suggested either population expansion or purifying selection.
The results of AMOVA suggested that both weevil species lacked spatial genetic structure at regional, continental, and intercontinental spatial scales (tables S3–S5 in the supplementary material). Differences within groups accounted for most of the genetic variances for the maize weevil (mtDNA 86.6%, P < 0.001, table S3 in the supplementary material; ITS 99.6%, P = 0.33, table S4 in the supplementary material) and the rice weevil (mtDNA 86.2%, P < 0.001, table S5 in the supplementary material).
Haplotype networks
Genealogical relationships among haplotypes of the maize weevil and haplotypes of the rice weevil were assessed using MJ networks, which were constructed for the concatenated mitochondrial genes (fig. 2) and ITS region (fig. 3), separately.
Median-joining network for the mitochondrial haplotypes of maize weevil (Sitophilus zeamais) and rice weevil (Sitophilus oryzae). The network was based on a dataset containing two gene regions concatenated: the cytochrome oxidase subunit I (COI) gene and the cytochrome oxidase subunit II (COII) gene. A circle represents a given haplotype (coded with numbers); circle size is proportional to the relative frequencies. Numbers of mutational steps are indicated with small (empty) circles when more than one (unless indicated otherwise). Color codes according to the geographic origin of the sequence.
Median-joining network for the ITS haplotypes of maize weevil (Sitophilus zeamais) and rice weevil (Sitophilus oryzae). A circle represents a given haplotype (coded with letters); circle size is proportional to the relative frequencies. Numbers of mutational steps are indicated with small (empty) circles when more than one (unless indicated otherwise). Color codes according to the geographic origin of the sequence.
There were 28 haplotypes for the concatenated set of the mitochondrial genes (fig. 2). Amongst the 28 haplotypes, 18 were recovered from specimens of the maize weevil, while the remaining ten haplotypes were from specimens of the rice weevil. For the maize weevil, haplotype 1 was recovered from 47.5% of the specimens and was the most widely distributed haplotype, with worldwide occurrence. For the rice weevil, the most frequent haplotype (haplotype 26; 48%) was also the most widely distributed haplotype. From both weevil species, low-frequency haplotypes were geographically restricted. The topology of the mitochondrial network showed a clear differentiation between the two weevil species. Most notably, there was no haplotype sharing between the two species. Each species exhibited its own distinct mitogroup (i.e., a group of closely related mitochondrial haplotypes). There were 167 mutational steps (from haplotype 9 to either haplotype 20 or 23) between the two mitogroups.
When compared to the mitochondrial network, the ITS network exhibited much less haplotype diversity. There were only seven haplotypes: sequences of the maize weevil displayed six haplotypes, while sequences of the rice weevil exhibited a single haplotype (fig. 3). There was no haplotype sharing between the two species of weevils; thus, each weevil species exhibited its own ribogroup (i.e., a group of closely related ITS haplotypes). The small sample size (n = 9) may account for the presence of a single ITS haplotype in the rice weevil ribogroup. There were ambiguities that we could not resolve using the criteria of Crandall & Templeton, (Reference Crandall and Templeton1993), therefore, we maintained these ambiguities in the network – showed as tip haplotypes having more than one connection to central haplotypes (fig. 3). Haplotypes A, B, C, and D (on the ribogroup of the maize weevil) and haplotype G (on the ribogroup of the rice weevil) displayed the highest frequencies and were distributed over a vast geographical area. Haplotypes E and F were singletons (i.e., they occurred only once in our dataset); these two haplotypes were obtained from a single specimen of the maize weevil from Thailand, which exhibited intragenomic polymorphism for ITS. Overall, there were no clear relationship between network architecture and geographic origin of the specimen (figs 2 and 3). In both networks, there was a tendency of finding haplotypes of medium to high frequency with origins within two to three ranges, while most of the haplotypes of low frequency were found within a single range.
Estimated date of divergence
The Beast analysis estimated the time of divergence between the maize weevil and the rice weevil. The results suggested that the two weevil species split from the tMRCA about 8.7 My ago (95% HPD: 4.0–15.0).
Microsatellite-based structure in maize weevil
Our analyses used eight out of the ten microsatellite loci we tested. Locus 1B10 was monomorphic, while locus 1B1 produced bands from a very low number of specimens; both loci were removed from subsequent analyses. After Bonferroni`s correction, no linkage disequilibrium was observed. The mean number of alleles was 3.7, ranging from 3.0 (population BR2) to 4.9 (population BR7); average heterozigosity (expected/observed) was 0.58/0.44 (table 3, table S6 in the supplementary material). The observed heterozygosity in all maize weevil populations was low, with inbreeding coefficients (F) varying from 0.06 (population BR3) to 0.46 (population THA). Overall mean was 0.23. A total of 19 private alleles were identified. These private alleles were found in ten populations: BR1, COL, MEX, and THA (a single private allele per population); BR4, BR6, and USA (two private alleles per population); and BR7, CHI, and MOZ (three private alleles per population) (table 3). Results of the AMOVA revealed that differences within populations accounted for 81.5% of the total variance. The genetic differentiation was moderate among all populations (F ST = 0.115, G ST = 0.110, R ST = 0.105).
Estimates of genetic diversity of Sitophilus zeamais based on eight microsatellite loci (see table S6 in the supplementary material for raw data), with mean number of individuals with locus successfully amplified (n), average number of alleles (A/locus), number of private alleles (A priv), allelic richness (AR), expected heterozygosity (He), observed heterozygosity (Ho), and inbreeding coefficient (F).
The analysis carried out on structure revealed that the 309 specimens of the maize weevil contained ancestry in three Bayesian groups (best K = 3). Overall, specimens from any given population displayed varying proportions of membership coefficients in the three Bayesian groups – depicted in gray, orange, and blue; respectively (fig. 4). Most of the 15 populations contained specimens that displayed proportion of membership in at least two Bayesian groups: more noticeable the joint presence of the gray and orange Bayesian groups, which were widespread in the Americas and Asia. There were populations in which most of the specimens exhibited a high proportion of membership in a single Bayesian group. For example, four populations from the Americas (EUA, MEX, PAN, and BR7) together with two populations of Asia (CHI and THA) showed the highest proportion of membership assigned to the gray Bayesian group, which ranged from 83% (EUA) to 40% (BR7). Six populations from Brazil (BR1 to BR6) and the population from India (IND) showed the highest proportion of assignment in the orange group, ranging from 76% (BR4) to 52% (IND). The presence of the third Bayesian group (depicted in blue) was most obvious in populations COL (90%) and MOZ (81%), which were each collected from sites located very far apart: Colombia and Mozambique. These results provided evidence for the lack of spatial genetic structure at the regional, continental, and intercontinental levels and suggested high gene flow amongst maize weevil populations.
Clustering analyses based on eight microsatellite loci for 309 specimens of maize weevil (Sitophilus zeamais) and geographic distribution of the three Bayesian groups (coded gray, orange, or blue) across 15 study populations. Along the x-axis of each plot, each vertical bar represents a weevil specimen; along the y-axis, membership coefficient of a specimen for a Bayesian group represents the fraction of its genome that has ancestry in that Bayesian group. In the insert, the best K (K = 3) calculated according to the ΔK method (Evanno et al., Reference Evanno, Regnaut and Goudet2005).
Discussion
The maize weevil and the rice weevil are each a distinct species
Our molecular evidence suggests that speciation is complete; therefore, the maize weevil should not be considered as a derived strain of its ‘sibling species’ – the rice weevil. Although they share a great deal of morphological similarities (Kuschel, Reference Kuschel1961; Halstead, Reference Halstead1963), the maize weevil and the rice weevil are reproductively isolated and therefore distinct species (Hidayat et al., Reference Hidayat, Phillips and Ffrench-Constant1996). Molecular evidence for complete speciation comes from the fact that there was no haplotypes shared between species; moreover, the networks displayed discontinuous distributions of haplotypes, with large gaps between groups within each network. Each species exhibited its own mitogroup (for the mitochondrial genome) and ribogroup (for the nuclear genome). The large mutational distances – 167 steps between mitogroups; 39 steps between ribogroups – indicate that gene flow between species has ceased completely. Since the time of their split from the tMRCA, the two species evolved without further genetic admixture. Over time, non-shared mutations accumulated within each mitogroups (or ribogroup) and gave rise to the large gaps uncovered on the networks. The large pairwise differences in the COI genes allowed for the use of molecular tools to discriminate between weevil species (Hidayat et al., Reference Hidayat, Phillips and Ffrench-Constant1996; Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013).
Ancient origin of the maize weevil
Molecular data provide support for an ancient origin of the maize weevil – about 8.7 My ago, with the 95% HPD lower bound at 4.0 My ago, which is significantly older than the establishment of agriculture between 10,000 and 6000 yBP. Our estimated time of divergence is very coarse and should be taken with caution, given that it was obtained with sequence data from a single gene set and used a mean substitution rate (0.0177 substitutions per site per My) estimated for the COI gene of beetles (Papadopoulou et al., Reference Papadopoulou, Anastasiou and Vogler2010). Estimating precisely when the maize weevil and the rice weevil split from the tMRCA is not a simple task. Although with caveats, our estimated time of divergence is valuable because it shed light on the temporal placement of the origin of the maize weevil relative to the onset of agriculture. This is compelling evidence that both the maize weevil and the rice weevil emerged as independent species prior to the onset of agriculture.
Additionally, in order to constrain the origin of the maize weevil within the time frame maize cultivation left Mesoamerica (about 400 years ago), the divergence dating analysis would require the mean substitution rate for the COI gene of Sitophilus to run 20,000 faster (354 substitutions per site per My) than the fastest known rate for the COI gene of beetles (0.0177 substitutions per site per My) (Papadopoulou et al., Reference Papadopoulou, Anastasiou and Vogler2010). This is compelling evidence that the maize weevil emerged prior to the human-mediated dissemination of maize as part of the Columbian Exchange (Crosby, Reference Crosby2003).
Our results do not allow us to deduce where the ancestral ranges of the maize weevil resided. A body of evidence suggests that both Southeast Asia and the Indian subcontinent are likely regions where the maize weevil may have emerged. These regions provide home for the remaining, non-pest species of Sitophilus (Plarre, Reference Plarre2010). Archeological findings attest the early association of Sitophilus with human activities in Southeast Asia. Deposits of ca. 3000 yBP in China contained Sitophilus specimens, which very likely are rice weevils (Chu & Wang, Reference Chu and Wang1975; Buckland, Reference Buckland1981). Pottery dated from ca. 10,500 yBP in Japan showed impressions of Sitophilus, which the authors claimed as belonging to maize weevils (Obata et al., Reference Obata, Manabe, Nakamura, Onishi and Senba2011). We hypothesize that the early cultivation and storage systems of rice and other cereal on Asia may have led to a habitat shift in the grain, from the prior acorn-feeding on forest environments to the later stored grains from early agricultural settlements.
Recent association between maize and maize weevil
Most of today's leading cereals were domesticated on the Old World; maize is an exception. Maize has been domesticated in the Mexican highlands about 9000 yBP (Piperno et al., Reference Piperno, Ranere, Holst, Iriarte and Dickau2009; Matsuoka et al., Reference Matsuoka, Vigouroux, Goodman, Sanchez, Buckler and Doebley2002). During the next few millennia, its cultivation quickly spread northwards into North America and southwards into the Andes and the Atlantic coast of South America, including Brazil (Iriarte et al., Reference Iriarte, Holst, Marozzi, Listopad, Alonso, Rinderknecht and Montaña2004; Pohl et al., Reference Pohl, Piperno, Pope and Jones2007). Currently, maize has the broadest cultivation range of all cereals (Van Heerwaarden et al., Reference van Heerwaarden, Doebley, Briggs, Glaubitz, Goodman, Gonzalez and Ross-Ibarra2011; Mir et al., Reference Mir, Zerjal, Combes, Dumas, Madur, Bedoya, Dreisigacker, Franco, Grudloyma, Hao, Hearne, Jampatong, Laloë, Muthamia, Nguyen, Presanna, Taba, Xie, Yunus, Zhang, Warburton and Hearne2013). The timing of which maize weevil became associated with maize is intriguing, given that the origin of the pest insect took place likely on Asia, whereas maize originated on Mesoamerica.
Archeological records that could attest the presence of Sitophilus on Mesoamerica prior to the Columbian Exchange are silent. Without human intervention, it is highly unlikely that the maize weevil could leave its ancestral ranges on Asia, expand its geographical range to Mesoamerica through a transoceanic journey, and finally find suitable habitats at early maize cultivation systems. Had the maize weevil left Asia and arrived on Mesoamerica prior to – or immediately after – maize domestication, our microsatellite analyses would had uncovered the signatures of this geographic range expansion on the population structure of the maize weevil. In contrast, we found no molecular evidence to support the Americas as a source of genetic variation for other geographic regions; there was no phylogeographic signals preserved on a region-dependent manner. The fact that genetic variation of the American populations reached values similar to those found elsewhere is consistent with the Americas as part of recent colonization events associated with high gene flow (see next section for details). Thus, we hypothesize that the maize weevil became associated with maize only very recently; this association took place after maize became widespread as a global crop. It is plausible that the relatively large grain size of maize and the warmer conditions in which maize is cultivated favored the maize–maize weevil association (Van Heerwaarden et al., Reference van Heerwaarden, Doebley, Briggs, Glaubitz, Goodman, Gonzalez and Ross-Ibarra2011; Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013; Mir et al., Reference Mir, Zerjal, Combes, Dumas, Madur, Bedoya, Dreisigacker, Franco, Grudloyma, Hao, Hearne, Jampatong, Laloë, Muthamia, Nguyen, Presanna, Taba, Xie, Yunus, Zhang, Warburton and Hearne2013).
Weak spatial structure at diverse geographic scales
Microsatellite marker analyses yielded concordant results with DNA sequence data. However, microsatellite markers analyses provided additional details about the distribution of genetic diversity of the maize weevil across cultivation areas at the regional, continental, and intercontinental scales that were not noticeable from DNA sequence data alone.
Sampling maize weevil across Brazil for microsatellite analyses – 142 specimens from seven populations that cover most of the maize cultivation area in the country – allowed us to investigate extant levels of gene flow among populations at the regional level. The lack of spatial genetic structure in Brazil is consistent with a pest species that: (a) has been introduced recently – more likely during the last few hundred years; (b) is capable of dispersing to short distances owing to its flying capabilities (Throne & Cline, Reference Throne and Cline1989; Rees, Reference Rees, Subramanyam and Hagsrum1996); (c) is capable of reaching new, distant territories owing to human-mediated transportation of infested grains (Longstaff, Reference Longstaff1981; Throne & Cline, Reference Throne and Cline1989; Plarre, Reference Plarre2010). High levels of gene flow and low genetic differentiation among populations constitute a pattern that very likely is not restricted to the maize weevil populations from Brazil; his pattern is probably widespread at a regional scale across other maize cultivation areas around the world (Semeao et al., Reference Semeao, Campbell, Beeman, Lorenzen, Whitworth and Sloderbeck2012; Coelho-Bortolo et al., Reference Coelho-Bortolo, Mangolin and Lapenta2016; Thagaraj et al., Reference Thagaraj, McCulloch, Subbarayalu, Subramaniam and Walter2016).
To explore the spatial genetic structure of the maize weevil further, we analyzed data from additional populations from the Americas (USA, MEX, PAN, and COL), Africa (MOZ), and Asia (IND, CHI, and THA). The lack of spatial genetic structure that has been detected at the regional level (Brazilian populations) was also found at the continental and intercontinental levels, with two Bayesian groups dispersed worldwide. The restricted distribution of one of the three Bayesian groups (coded blue in fig. 4) indicated a recent gene flow between Colombia and Mozambique. Direction of this gene flow could not be determined using our current dataset, but the Colombia–Mozambique pattern of gene flow is consistent with the major role that Colombia played in the Portuguese slave trade possibly leading to maize introduction from Colombia to Africa, including Mozambique (Madeira Santos & Ferraz Torrão, Reference Madeira Santos, Ferraz Torrão and Chastanet1998; Mir et al., Reference Mir, Zerjal, Combes, Dumas, Madur, Bedoya, Dreisigacker, Franco, Grudloyma, Hao, Hearne, Jampatong, Laloë, Muthamia, Nguyen, Presanna, Taba, Xie, Yunus, Zhang, Warburton and Hearne2013). Altogether, these findings suggest that the maize weevil underwent recent range expansions at a global scale, most likely mediated through human-related activities.
The lack of spatial genetic structure at varying geographic levels has implications for pest management. Physiological and behavioral differences among populations of the maize weevil have been reported and such differences incur in likely differences in dispersion and damage potential to the host grains (Morales et al., Reference Morales, Cardoso, Della Lucia and Guedes2013; Carvalho et al., Reference Carvalho, Vieira, Haro, Corrêa, Ribon, Oliveira and Guedes2014; Malia et al., Reference Malia, Rosi-Denadai, Cardoso and Guedes2016). However, such differences seem more prevalent within populations rather than between populations, suggesting a high multidirectional gene flow, which is consistent with the findings we reported herein and with recent studies using sets of integrated behavioral tendencies in populations of the maize weevil (Morales et al., Reference Morales, Cardoso, Della Lucia and Guedes2013; Malia et al., Reference Malia, Rosi-Denadai, Cardoso and Guedes2016).
Phytosanitary concerns for the maize weevil
Grain weevils are very damaging to stored cereals. While the maize weevil is particularly important as a pest of stored maize on both Neotropical America and Africa, the rice weevil is important for stored maize in temperate climate regions such as North America, Europe, and Asia (Longstaff, Reference Longstaff1981; Rees, Reference Rees, Subramanyam and Hagsrum1996; Plarre, Reference Plarre2010; Corrêa et al., Reference Corrêa, Oliveira, Braga and Guedes2013). The maize weevil exhibits a range of trait variations to support its status as insect pest, including flight proficiency, walking ability, body mass, grain consumption, and genes for insecticide resistance (Grenier et al., Reference Grenier, Nardon and Nardon1994; Daglish, Reference Daglish2004; Guedes et al., Reference Guedes, Oliveira, Guedes, Ribeiro and Serrão2006, Reference Guedes, Tolledo, Corrêa and Guedes2010; Corrêa et al., Reference Corrêa, Pereira, Cordeiro, Braga and Guedes2011). The unintended introduction of novel genotypes may increase the fitness of local populations of the maize weevil because the newly introduced specimens may carry a vast array of favorable alleles. Upon mating, reshuffling of the newly acquired genes together with the genes of the local specimens may give rise to novel genomic combinations, which otherwise would not take place on either local or regional populations of the maize weevil. Thus, phytosanitary agencies should pay attention to prevent, at the maximum extent possible, the introduction of novel genetic variation, especially genetic variation with origin from populations that might have experienced high selection pressure toward adaptive traits.
Enhanced dispersal ability and insecticide resistance are some of the traits that may potentially aggravate grain losses when specimens that harbor them are re-introduced into occupied territories or introduced into newly occupied areas. Insecticide use and insecticide resistance of maize weevil for instance are frequent problems in Neotropical America (Champ & Dyte, Reference Champ and Dyte1977; Subramanyam & Hagstrum, Reference Subramanyam and Hagstrum1996; Ribeiro et al., Reference Ribeiro, Guedes, Oliveira and Santos2003; Pimentel et al., Reference Pimentel, Faroni, Guedes, Sousa and Tótola2009; Haddi et al., Reference Haddi, Mendonça, Santos, Guedes and Oliveira2015). Unintended human-mediated transport may disperse adaptive traits, eventually leading to the establishment of insecticide resistant genotypes of the weevils elsewhere. Therefore, phytosanitary measures to prevent further transfer of both the maize weevil and the rice weevil, particularly from warm maize producing areas, are worthwhile and should be implemented.
In conclusion, we uncovered evidence for an ancient origin of both the maize weevil and its morphologically indistinguishable congener, the rice weevil; both species likely emerged prior to the beginning of the agriculture on Asia. The association between the maize weevil with maize, its preferred host, was established recently, after maize cultivation became widespread. Lack of spatial structure characterized the extant populations of the maize weevil, which is suggestive of high gene flow.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485316000687.
Acknowledgements
We would like to express our deepest gratitude to a multitude of colleagues, farmers and technicians who assisted us during field sampling (most listed in supplementary tables S1 and S2). We also thank Prof. Anete de Souza and her group at the CBEMEG – University of Campinas (Campinas, SP, Brazil) for assisting us with the development of microsatellite markers and EMBRAPA Café at the Center of Biotechnology Applied to Agriculture of the Federal University of Viçosa (Viçosa, MG, Brazil) for providing us with facilities for microsatellite genotyping. The financial support was provided by the Minas Gerais State Foundation of Research Aid – FAPEMIG (grants PPM 00561-15 to L.O.O. and APQ-01977-09 to R.N.C.G.), National Council of Scientific and Technological Development – CNPq (fellowships 305827/2015-4 to L.O.O. and PQ 304174/2010-6 to R.N.C.G.), and the CAPES Foundation (grant AUXPE 191/2009 to R.N.C.G.). A.S.C. received a fellowship from CNPq.
