The Bemisia tabaci (Gennadius) (Hemiptera: Aleyroridae) sibling species group is comprised of genetic variants defined by biological differences or a fragment of mitochondrial cytochrome oxidase I gene sequence (mitotype) that allows for phylogeographic affiliation. Some mitotypes may cause damage to crop plants by feeding and transmission of plant viruses. In Pakistan, cotton-vegetable agroecosystems are vulnerable to whitefly-transmitted virus (genus Begomovirus; family Geminiviridae) infection. The identity and distribution of the whitefly B. tabaci mitotypes associated with the cotton crop were studied in 8 districts in the Punjab Province from 2014 to 2016. Phylogenetic analysis of the 3-fragment of the mitochondrial cytochrome oxidase I gene indicated the predominant haplotypes belonged to the Asia II-1 mitotype, with pairwise distances ranging from 0.15 to 3.2%. Pairwise distances showed that B. tabaci haplotype diversity varied by district, with the Khanewal harboring the highest divergence at 1.37%, compared to the lowest at 0.50% in the Dera Ghazi Khan district. The median-joining network analysis showed genetic expansion, or a ‘recovery’ trend, following the declining genetic diversity that occurred during the late 1990s to the early 2000s. The Asia II-1 mitotype group was the predominant whitefly vector species in Punjab Province. The haplotype network provides documentation of continued genetic expansion among the B. tabaci populations in the Punjab, which is consistent with previously reported trends among whiteflies sampled in the same or nearby districts from 2012 to 2014. Genetic expansion varied among districts and could be explained by factors unique to each district, i.e., management practices that influence B. tabaci mitotype composition, whitefly susceptibility to cotton leaf curl disease complex, and cotton genotype
The Bemisia tabaci (Gennadius) (Hemiptera, Aleyrodidae) sibling (or cryptic) species group is comprised of an unknown number of morphologically indistinguishable phenotypic variants, for which a small number have been characterized biologically (Bedford et al. 1994; Brown et al. 1995; Sseruwagi et al. 2005, 2006; Caballero 2007; Gill & Brown 2010; Hadjistylli et al. 2016). As a species group, they are endemic to tropical, subtropical, and temperate ecosystems. As a group, they also have an extensive host range among eudicots (Mound & Halsey 1978; Cock 1993), with more than 200 plant species recorded as hosts in Pakistan (Attique et al. 2003).
Recently, B. tabaci mitotypes and their genetic groups have been differentiated based on phylogenetic analysis of an approximately 580 to 780 base pair fragment of the 3′-end of the mitochondrial cytochrome oxidase I gene (mtCOI-??3′) sequence (Frohlich et al. 1999), with several theoretical estimates proposed for species delineation at 3.5 to 4.0% nucleotide (nt) divergence, irrespective of corresponding biological species criteria (Dinsdale et al. 2010; Lee et al. 2013). To date, no consensus has been reached for a globally applicable species demarcation system within the B. tabaci taxon. Even so, understanding the mitotype composition and population structure of B. tabaci in agroecosystems is important, because management practices may cause shifts in the prevalence of the different B. tabaci mitotypes, which in turn may govern begomovirus transmission competency, potentially influencing the dynamics of the different mitotypes and cotton leaf curl disease outbreaks, possibly resulting in shifting of 1 or more components of the pathosystem. Different mitotypes of the B. tabaci sibling species group have been shown to respond differently to various pesticide regimes, the range of plants suitable as hosts, and environmental factors. These and other factors strongly influence the composition of B. tabaci mitotypes and their dynamics in agroecosystems (Bedford et al. 1994; Brown 2007, 2010; Gill & Brown 2010).
Plant viruses transmitted by whiteflies belong to 1 of 5 genera (Jones 2003; Navas-Castillo et al. 2011). Among the economically important plant viruses transmitted by whiteflies, the genus Begomovirus (Geminiviridae) is the most widespread, comprising more than 400 species ( https://talk.ictvonline.org/taxonomy/), including many known to cause economically important plant diseases (Brown et al. 2012). In Pakistan and India, over 60 begomoviral species or strains have been described from cultivated and wild plant species (Ho et al. 2017). Cotton leaf curl disease is caused by one of several monopartite begomoviruses or strains, and associated betasatellites that are endemic to much of southern Asia, including Pakistan (Zhou et al. 1998; Saleem et al. 2016). During the past several decades, the emergence of cotton leaf curl disease has resulted in major economic losses in the cotton producing areas of Pakistan (Briddon et al. 2001; Zubair et al. 2017). Understanding the distribution of the different whitefly mitotypes and their population dynamics, which can be predicted by spatio-temporal patterns of mtCOI diversity, have become essential for devising effective management approaches. Such strategies are further enhanced by knowledge of the host range, pesticide resistance, and other biological characteristics associated with the range of different mitotypes that transmit the equally diverse begomoviral species and strains that cause cotton leaf curl disease, as well as mitotype-virus specific factors that directly or indirectly influence the efficiency of whitefly-mediated transmission.
The objective of this study was to determine the mitotype composition and spatial patterns of genetic divergence of the B. tabaci whitefly vector collected from commercial cotton fields in the southern Punjab districts of Pakistan from 2014 to 2016. Mitotypes of the whitefly B. tabaci were identified phylogenetically based on the mtCOI-3′ gene sequence, and the structure of haplotypes was studied using a median-joining network analysis.
Materials and Methods
WHITEFLY FIELD COLLECTIONS
Whiteflies were collected from infested cotton plants in 8 districts of the southern Punjab Province, including Bahawalpur, Dera Ghazi Khan, Khanewal, Lodhran, Multan, Muzaffargarh, Rahim Yar Khan, and Vehari, Pakistan, from 2014 to 2016 (Fig. 1). Adult whiteflies were collected using a hand-held aspirator, transferred to a microfuge tube containing 95% ethanol, and stored at –20 °C. Standard management practices and the same cotton varieties were implemented at all of the whitefly collection sites. Twenty-four whiteflies per district (192 total) were sequenced, and haplotypes were assigned to a sister or major clade based on the mtCOI sequence fragment.
WHITEFLY DNA ISOLATION, POLYMERASE CHAIN REACTION AMPLIFICATION, AND DNA SEQUENCING
Total DNA was purified from individual whiteflies using a previously published protocol (Zhang et al. 1998) with modifications (Paredes-Montero et al. 2019). Briefly, whiteflies were dipped in double-distilled water to dilute the ethanol, followed by blotting against the torn edge of a piece of No. 1 Whatman filter paper (Sigma-Aldrich, Darmstadt, Germany) to remove the distilled water. Each whitefly was ground with a micro-pestle in the bottom of a 1.0 mL microfuge tube, containing 250 µL lysis buffer. The lysis buffer contained 100 mM Tris-HCl, 1.4 M NaCl, 20 mM EDTA, 2% hexadecyltrimethylammonium bromide (CTAB) (Sigma-Aldrich, Darmstadt, Germany), and 20 µg per reaction of proteinase K. The lysate was incubated at 55 °C for 1 h followed by 10 min at 65 °C. One volume of chloroform:isoamyl-alcohol (24:1) was added, the contents were mixed by inversion, and the emulsion was separated by centrifugation, at 14,000 rpm for 3 min. The DNA was precipitated by the addition of 1 volume isopropanol and 40 µg glycogen, followed by incubation at 4 °C for 10 min. The pellet was collected by micro-centrifugation at 14,000 rpm for 10 min, washed with 70% ethanol, air dried, dissolved in 30 µL ddH2O, and stored at –20 °C.
To detect the possible presence of the exotic B mitotype [North Africa-Mediterranean-Middle East major clade [NAF-MED-ME (III)] among the endemic Asia I and II mitotypes (major clades) (Brown 2010), samples were subjected to polymerase chain reaction amplification using the previously reported B-mitotype-specific primers 5′-BtBF1F-3′ [TATTTCACTTCAGCCACTATAA] and 5′-WfBr2R-3′ [GCTTAAATCTTACTAACCGCAG], yielding an expected size fragment of 550 base pairs (Andreason et al. 2017). The negative samples were subjected to polymerase chain reaction amplification using the universal mtCOI primers C1-J-2195F 5′-TTGATTTTTTGGTCATCCAGAAGT-3′ and L2-N-3014R 5′-TCCAATGCACTAATCTGCCATATTA-3′ that yield an expected size amplicon of about 850 base pairs for B. tabaci (Frohlich et al. 1999). The polymerase chain reaction cycling parameters were 95 °C for 2 min, followed by 30 cycles at 95 °C for 60 s, 52 °C for 60 s, 72 °C for 60 s, with a final extension at 72 °C for 5 min. Reactions contained 1X Jumpstart REDTaq Ready-Mix (Sigma-Aldrich, St. Louis, Missouri, USA), 0.4 µM primers, 20 ng DNA, and ddH2O, to a final volume of 25 µL.
The amplicons were ligated into the pGEM®-T Easy plasmid vector (Promega Corp., Madison, Wisconsin, USA). The presence of an expected size insert of 875 to 880 base pairs (non-B mitotype) and 550 base pairs (B mitotype), respectively, was confirmed by colony polymerase chain reaction amplification (Gussow & Clackson 1989). Amplicons were sequenced bi-directionally by capillary Sanger DNA sequencing (University of Arizona Genetics Core, University of Arizona, Tucson, Arizona, USA). The sequences were assembled, trimmed, manually edited, and aligned using Lasergene software (DNASTAR, Madison, Wisconsin, USA).
SEQUENCE EVALUATION FOR PRESENCE OF NUCLEAR INSERTIONS OF MITOCHONDRIAL DNA SEGMENTS
To identify amplified nuclear insertions of mitochondrial DNA segments (Song et al. 2008), mtCOI sequences were examined for nuclear insertions of mitochondrial DNA signatures, including indels, stop co-dons, and multiple base calls on chromatograms (Song et al. 2008). Whitefly sequences from the same sample having inconsistent BLASTn results (GenBank) were removed from the alignment. The redundant haplotypes sharing 100% identity were identified using FABOX v1.5 (Aarhus University, Aarhus, Denmark) (Villesen 2007) ( http://usersbirc.au.dk/palle/php/fabox/) and removed (collapsed) to obtain the final set of mtCOI sequences consisting of 1 representative sequence for each unique haplotype.
MITOTYPE IDENTIFICATION BY PHYLOGENETIC ANALYSIS
Mitotypes were identified by phylogenetic analysis of about 780 base pair fragments of mtCOI-3′, with selected reference COI sequences from the J.K. Brown Laboratory database. Reference sequences were representative of the 2 major clades of B. tabaci known to occur in Pakistan, the Asia [ASIA (IV)], and the Asia-Australia-Pacific [AS-AUS-PAC (V)] (Ashfaq et al. 2014; Masood et al. 2017; Islam et al. 2018).
The mtCOI-3′ sequences were aligned using MUSCLE (Edgar 2004), implemented in MEGA v6 (Tamura et al. 2013). The alignment was trimmed to result in a sequence length of about 760 base pairs. The phylogenetic tree was reconstructed using Maximum Likelihood, implemented in MEGA v6 software, with 1,000 bootstrap iterations (Tamura et al. 2013). The best-fitting model of evolution was determined using jmodeltest v2.1.7 (Free Software Foundation, Boston, Massachusetts, USA) (Darriba et al. 2012, and the decision theory selection method, available in CIPRES ( www.phylo.org). The best-fitting model was identified as Hasegawa-Kishino-Yano (Hasegawa et al. 1985), with gamma distributed-rate variation among sites (HYK+G). The tree was rooted using the mtCOI-3′ sequence for Bemisia afer (Priesner & Hosny) (Hemiptera: Aleyrodidae) (J.K. Brown Laboratory, COI database), and bootstrap values were placed at the major nodes on branches having a bootstrap value of ≥ 70%. The tree was drawn using FigTree v1.4.2 (University of Edinburgh, Edinburgh, United Kingdom; http://tree.bio.ed.ac.uk/software/figtree/). The map showing mitotype distributions was drawn using the CorelDraw X7 software (Corel Corporation, Ottawa, Ontario, Canada; https://www.coreldraw.com/en/pages/coreldraw-x7/).
GENETIC DIVERSITY
The extent of genetic diversity among the mtCOI-3′ sequences was estimated using a corrected distance analysis that calculates evolutionary distance based on the HYK+G model that accounts for multiple substitutions at a single site (Hasegawa et al. 1985). The genetic diversity was calculated using the ‘compute overall mean distance' option in the ‘distances’ menu, available in MEGA v6.0 (Tamura et al. 2013). The barplots of the average pairwise distance per districts were built using barplot function in the R version 3.5.2 (R Core Team 2019).
NETWORK ANALYSIS
The aligned matrix, which contained 168 mtCOI-3′sequences, was used to reconstruct a median-joining haplotype network (Bandelt et al. 1999). The sequences were input into the “nexus sequential” format to the Network v 5.0.0.3 algorithm (Fluxus Technology Ltd, Suffolk, United Kingdom) available at http://www.fluxus-engineering.com/sharenet.htm. The analysis was carried out by weighting each site uniformly, with the epsilon (η) value set to zero.
Results
MITOTYPE IDENTIFICATION
We obtained 168 sequences, which were collapsed to 107 unique haplotypes. Based on the phylogenetic analysis, all 107 unique mtCOI-3′ haplotypes (about 760 base pairs) grouped with reference sequences of the Asia II-1 mitotype (Fig. 1, Supplementary Table S1). The Asia II-1 mitotype and 11 other mitotype groups, arbitrarily designated as subclades 1 to 12 (Dinsdale et al. 2010; Firdaus et al. 2013; Masood et al. 2017), cluster within the ‘ASIA’ phylogeographic clade that spans all of the Asian continent, Australia, and the Pacific islands (Brown 2010) (Fig. 2).
GENETIC DIVERSITY AMONG ASIA II-1 HAPLOTYPES
The corrected pairwise distance among Asia II-1 haplotypes was 0.15 to 3.2%, with extent of divergence varying primarily by district (Fig. 3). Whiteflies collected in the districts of Muzaffargarh and Khanewal exhibited the highest pairwise divergence at 1.06 and 1.37%, respectively, whereas the range of pairwise distances for the other 6 districts was 0.50 to 0.77%, with the Dera Ghazi Khan district showing the lowest divergence at 0.50% (Fig. 3). The differences in genetic diversity are reflective of the dynamic patterns of diversification, which has been shown to influence mitotype distribution. Accordingly, when certain mitotypes show differences in transmission competency over others, the distribution and spread of cotton leaf curl disease-begomoviruses in cotton and vegetable crops can be altered dramatically, or result in the emergence of previously unrecognized species or strains, often due to recombination (Bisaro 1994; Brown 2001; Fauquet et al. 2005; Saleem et al. 2016).
NETWORK ANALYSIS FOR ASIA II-1 HAPLOTYPES
Based on the median-joining network, 4 predominant haplotypes were resolved within the Punjab region Asia II-1 population, and the haplotypes in cotton crops were found to be ‘geographically unrestricted' or independent of district. Also, several rare haplotypes were identified, based on the accumulation of unique mutations, relative to the ‘common’ haplotype. The rare haplotypes in the network were identified based on their affiliation with a star-like pattern surrounding the common haplotypes (Fig. 4). This kind of pattern can be indicative of population recovery, which is characteristically reflected as increased genetic diversity following a decline in diversity within a population (Fig. 4).
Discussion
The results indicated that Asia II-1 group sub-mitotype 1 was the only B. tabaci variant among the samples collected from the cotton-growing areas of the Punjab Province. This observation is consistent with recent studies that reported it to be the predominant mitotype in agroecosystems in the Punjab and Sindh provinces (Ashfaq et al. 2014; Masood et al. 2017; Islam et al. 2018; Paredes-Montero et al. 2019). Five other endemic mitotypes or sub-mitotypes have been documented previously in Pakistan agroecosystems, namely Asia I and Asia II-1, Asia II-5, Asia II-7, and Asia II-8 (Firdaus et al. 2013; Masood et al. 2017; Islam et al. 2018; Paredes-Montero et al. 2019). Additionally, there is one report of the B mitotype (Ashfaq et al. 2014) which is phylogeographically affiliated with the NAF-MED-ME (III) major clade, and so is thought to have been introduced only recently (Ashfaq et al. 2014).
The diversity of haplotypes given by the corrected pairwise distances varied by district. The most highly divergent haplotypes at 1.37% were found in Khanewal, whereas in Dera Ghazi Khan, haplotypes differed only by 0.5%. The range of genetic divergence among whiteflies from the different districts studied was found to be similar to divergences reported in several previous studies (Ashfaq et al. 2014; Masood et al. 2017; Paredes-Montero et al. 2019), suggesting that differences in whitefly management practices may have favored or limited establishment or displacement of certain haplotypes compared to others (Chu et al. 2010; Horowitz & Ishaaya 2014). Certain classes of insecticides over others have been shown to favor outbreaks of different B. tabaci mitotypes/sub-mitotypes, which appears to be due to differences in their ability to develop insecticide resistance (Bedford et al. 1994; Coats et al. 1994; Anthony et al. 1995; Denholm et al. 1998; Ahmad et al. 2002; Dennehy et al. 2010; Luo et al. 2010; Horowitz & Ishaaya 2014).
The network of haplotypes exhibited a high frequency of unique mutations that gave rise to several low-frequency haplotypes that diverged from the predominant haplotype in 1 or 2 mutations, forming a star-like structure (Mirol et al. 2008). This pattern is suggestive of genetic expansion or recovery from a recent bottleneck within the Asia II-1 populations. A previous study conducted in these same locations from 2012 to 2014 showed a similar scenario, with haplotype diversity levels starting to recover after a dramatic decline in genetic diversity among Asia II-1 haplotypes starting in the early 2000s (Paredes-Montero et al. 2019). It was hypothesized that this pattern of reduced genetic variation could have been due to the reasonably well-documented overuse of pesticides in cotton agroecosystems over the last 2 decades (Ahmad et al. 2002). Two yr after initial signs of recovery occurred, the rate of genetic expansion was about the same for B. tabaci collected from 2014 to 2016; however, the number of predominant haplotypes increased from 1 to 4 in relation to that reported by Paredes-Montero et al. (2019). In this context, and despite the smaller sample size analyzed in this study, at 107 compared to 538 samples collected in many of the same or additional districts analyzed by Paredes-Montero et al. (2019), the genetic signatures indicated the continued predominance of haplotypes classified in the subclade Asia II-1 in cotton grown in the Punjab Province from 2014 to 2016. Therefore, based on these observations, sustained genetic expansion of the populations is predicted to reach a plateau in subsequent yr in the absence of changes in cotton varieties or management practices.
The exclusive presence of Asia II-1 in cotton-agroecosystems in the southern Punjab Province may explain in part the recent spread of several predominant strains and species of cotton leaf curl disease there (Paredes-Montero et al. 2019). This hypothesis is supported by the high incidence of cotton leaf curl disease observed in cotton plantings where the Asia II-1 mitotype predominates. In contrast, lower disease incidence has been reported in the southern province of Sindh (Simón et al. 2003; Ahmed et al. 2011; Pan et al. 2018), where the B biotype occurs (Ashfaq et al. 2014) together with mitotypes endemic to the Asian continent (Ashfaq et al. 2014; Masood et al. 2017; Islam et al. 2018; Paredes-Montero et al. 2019). Based on one study, the (exotic) B mitotype showed low efficiency to no transmission of an isolate of the cotton leaf curl multan virus introduced recently in China (Pan et al. 2018), compared to variable rates of transmission by the introduced Q mitotype, or high rates of transmission by mitotypes endemic to Asia. This evidence for differential transmission specificity, together with a low transmission frequency of cotton leaf curl multan virus-transmission by the B mitotype suggests that, although the introduced B mitoype has established in at least 1 locale in Pakistan, it is unlikely to have a prominent role in the transmission of the majority of begomoviruses occurring in Pakistan. However, the B miotype has been shown to be a competent vector of other New and Old World begomoviruses, and so it could pose a threat as a begomovirus vector if any of those viruses are introduced into Pakistan in the future.
Conclusion
The Asia II-1 was identified as the predominant mitotype among B. tabaci samples collected from cotton in southern Punjab, Pakistan. Here, 2 lines of evidence underscore differences in diversification among and between different B. tabaci mitotypes associated with the different districts. First, the intra-mitotype diversity was found to vary by district, which most likely is indicative of differential mitotype responsiveness to locally imposed agricultural practices. Second, a continued increase in haplotype diversity, or genetic expansion, was observed following a notable decline in diversity beginning in about 2010 with the initiation of recovery during 2014 to 2016 that could in time lead to restored diversity, perhaps approaching the previously high levels, and the geographic redistribution of the Asia II-1 and other sub-mitotypes in the southern Punjab Province of Pakistan. Several interacting factors are thought to be involved, including altered management practices implemented since 2010, including use of different pesticide regimes for whitefly control, particularly those that cause minimal agroecosystem and natural enemy disturbances. Additionally, with the release of tolerant or resistant cotton leaf curl disease cotton varieties, transmission frequencies may have been reduced, which if so encourages the routinely judicious use of insecticides. Routine monitoring of the population structure and distribution of B. tabaci mitotypes and cotton leaf curl disease-begomoviruses in cotton or cotton-vegetable agroecosystems in the Punjab Province is expected to contribute to implementing best practices that contribute to avoiding upsurges in whitefly that lead to cotton leaf curl disease outbreaks. Consequently, the new information presented here as well as that from previous studies of historical fluxes in B. tabaci mitotype and begomovirus distribution are extremely valuable for predicting these outbreaks.
Acknowledgments
The authors wish to thank Muhammad Naeem Aslam of CAB International-South Asia for his assistance with whitefly sample collections. The first author expresses sincere gratitude to the Higher Education Commission of Pakistan for the International Research Support Initiative Program fellowship that provided funding to cover the cost of travel and lodging, and partial support of the research reported here, carried out at the School of Plant Sciences, University of Arizona, Tucson, Arizona, USA.