Pheromone‐binding proteins based phylogenetics and phylogeography of Maruca spp. from Asia, Africa, Oceania, and South America

Abstract Variations in the functional response of legume pod borer (Maruca vitrata) populations to sex pheromone blends were observed in Asia and Africa. Hence, this study was carried out to understand the differences in pheromone‐binding proteins (PBPs) among Maruca populations in Asia, Africa, Oceania, and South America. A de novo transcriptome assembly was adopted to sequence the entire transcribed mRNAs in M. vitrata from Taiwan. The raw‐sequence data were assembled using homologous genes from related organisms in GenBank to detect M. vitrata PBPs (MvitPBPs). Sections of the cDNA of MvitPBP of different length were used to design primers to amplify the full‐length cDNA of PBPs. All three PBP sequences comprised three exons interspersed by two introns. In total, 92 MvitPBP1 haplotypes, 77 MvitPBP2 haplotypes, and 64 MvitPBP3 haplotypes were identified in 105, 98, and 68 Maruca individuals, respectively. High pairwise F ST values (0.41–0.73) and phylogenetic analyses distinguished the putative Maruca species in South America from those occurring in rest of the world, and possibly two putative subspecies in Asia and Africa. The haplotype networks and Automatic Barcode Gap Discovery analyses also confirmed these results. The negative Tajima's D and Fu's F S values showed the recent demographic expansion of Maruca populations. Thus, this study confirmed the presence of different Maruca species and/or subspecies in different continents based on the diversity within PBP genes. Additional sampling and studies are suggested for Oceania and South America. The genetic differences among Maruca populations should be carefully considered while using sex pheromone lures and bio‐control agents.

M. vitrata caused more than 65% grain yield reduction in pigeon pea in Australia (Sharma, Saxena, & Bhagwat, 1999). Hence, farmers rely more on chemical pesticides to combat this pest. For instance, more than 80% of the yard-long bean growers in Cambodia, Lao PDR (Laos), Thailand and Vietnam predominantly relied on synthetic pesticides (Schreinemachers et al., 2017(Schreinemachers et al., , 2014. On an average, Thai yard-long bean growers used 16.3 kg/ha of pesticide formulations per cropping cycle (Schreinemachers et al., 2014), and Cambodian farmers mixed four pesticides together in a single spray (Schreinemachers et al., 2017). Such an intensive pesticide use has serious consequences on human and environmental health. Hence, alternative pest management strategies are warranted for legume growers.
Insect pheromones are an important component in pest management programs, especially as a monitoring, mating-disruption, and/or mass-trapping tool. M. vitrata sex pheromone consists of one major and two minor compounds (Adati & Tatsuki, 1999;Downham et al., 2003). A synthetic sex pheromone consisting of major [(E,E)-10,12-hexadecadienal] and minor [(E,E)-10,12-hexadecadienol and (E)-10-hexadecenal] compounds developed in a ratio of 100:5:5 attracted male moths in Benin and Ghana, whereas the major compound alone was most effective in Burkina Faso (Downham et al., 2003(Downham et al., , 2004. However, none of these blends attracted any males in Taiwan (Schläger et al., 2012), Thailand, and Vietnam (Srinivasan et al., 2015), although a variant blend was attractive in India (Hassan, 2007). These differential responses suggest the presence of genetically different M. vitrata populations.
An earlier study showed evidence for the presence of multiple Maruca species or subspecies (Margam et al., 2011). Herbison-Evans, Hacobian, and Crossley (2017) also reported two forms of M. vitrata in Australia. However, we undertook a detailed study investigating the mitochondrial cytochrome c oxidase I (coxI) diversity in populations from Southeast Asia (the probable center of origin for Maruca), South Asia, sub-Saharan Africa, and in reference populations from Oceania and Latin America. This study confirmed the presence of three putative Maruca species, including one in Latin America, one in Oceania (including Indonesia) and M. vitrata in Asia, Africa and Oceania (Malini, Schafleitner, Muthukalingan, & Ramasamy, 2015).
The results also showed the presence of two putative M. vitrata subspecies in Asia and Africa.
Since different species or subspecies seem to exist in the genus Maruca, the pheromone composition and their reception may not be uniform in different geographical locations. A recent study found only two pheromone compounds in M. vitrata populations from Taiwan, Thailand, Vietnam, and Benin (Schläger et al., 2015).
Similarly, different M. vitrata populations also produce pheromone compounds in different ratio. M. vitrata females from Wuhan and Huazhou provinces in China produced different ratio of the three compound pheromones (Lu, Qiao, & Luo, 2013). Thus, the pheromone composition in M. vitrata seems to vary across locations.
Hence, it has been hypothesized that variations in the M. vitrata male pheromone reception may be attributed to the presence of different pheromone strains in M. vitrata females.
Insect sex pheromones facilitate the mate-finding among the members of an insect species. In male moths, a specialized subset of chemosensilla contains pheromone-sensitive neurons, which are highly sensitive and specific to sex pheromone compounds produced by conspecific females (LaForest, Prestwich, & Löfstedt, 1999). At the molecular level, the reception of pheromones in male moths is mediated by pheromone-binding proteins (PBPs), a subfamily of odorantbinding proteins (OBPs). PBPs which are localized in the lymph of the sensilla surrounding the olfactory neuron cells on the moth antennae (Vogt, Rogers, Franco, & Sun, 2002) bind to the lipophilic pheromonal compounds (Bette, Breer, & Krieger, 2002;Lautenschlager, Leal, & Clardy, 2007;Maida, Ziegelberger, & Kaissling, 2003;Steinbrecht, Laue, & Ziegelberger, 1995;Vogt & Riddiford, 1981) and carry them to the receptor cells (Van den Berg & Ziegelberger, 1991). It has been demonstrated that the change in male pheromone response behavior is caused by differences in a sex-linked locus or set of linked loci (Willett & Harrison, 1999). The gene loci that are instrumental in conferring specificity in pheromone communication systems should show fixed amino acid differences between strains or species (Willett & Harrison, 1999). Thus, understanding the patterns of variation in the gene encoding PBP could provide insights into the population structure of Maruca spp., which differed in their responses to the same pheromone blend(s) in different geographical locations. Hence, this study was carried out to assess whether there are fixed nucleotide differences at the PBP locus between the pheromone strains of Maruca from different host plants and geographical origin.

| Insects
A Maruca vitrata colony was established at the Insectary of World Vegetable Center from a field population. The larvae were reared on Spodoptera exigua meridic diet (Bio-Serv, Frenchtown, NJ, USA) modified with cowpea powder, at 27 ± 1°C and 70 ± 10% relative humidity, photoperiod 14:10 hr (Light:Dark) until pupation. On pupation, they were sexed and placed in acrylic cylinders (30-cm long and 15-cm diameter), whose ends were covered with nylon-nets. Emerged adults were fed with 10% (w/v) sugar solution. Besides from Taiwan, M. vitrata larval populations from nine countries (Bangladesh, Benin, Indonesia, India, Kenya, Laos, Malaysia, Thailand, and Vietnam) from different host plants were collected (Malini et al., 2015).

| RNA extraction, complementary DNA (cDNA) synthesis and reverse transcription polymerase chain reaction (RT-PCR) amplification
About 100 antennae were used to obtain about 25 mg of the tissues.

| Molecular divergence and population genetic analyses
The MvitPBP1, MvitPBP2, and MvitPBP3 sequences were aligned and edited using BioEdit v7.0 (Hall, 1999 quences to identify if a population has undergone a recent population expansion event and is determined by the difference between average number of nucleotide differences and the number of segregating sites estimated from pairwise comparisons (Tajima, 1989).
Fu's F S test uses information from the haplotype distribution in a sample. The test estimates the probability of observing a random sample with equal or less singletons than the observed given a level of diversity. The test is based on the infinite site mutation model and assumes that all of the alleles are selectively neutral.
The genetic structure of M. vitrata populations based on various PBP sequences was examined by analysis of molecular variance (AMOVA) using Arlequin 2.001 (Schneider, Roessli, & Excoffier, 2000). This method was used to partition the genetic variance within and among populations as well as within and among groups. The populations were grouped by geographical locations (continents). Levels of significance were determined through 1,000 random permutation replicates. Pairwise F ST values used to appraise the genetic structure among populations were obtained with 1,000 permutations and at the significance level of 0.05 using the K2P model (Kimura, 1980).

| Phylogenetic, species delineation, and haplotype network analyses
The FASTA formatted coding regions of MvitPBP sequences were imported into the MEGA-X software package sequence alignment application, and a multiple sequence alignment was performed with the ClustalW algorithm using default parameters (Tamura et al., 2011).
The aligned sequences were used for phylogenetic analysis. The evolutionary history among the haplotypes of MvitPBP sequences was inferred by using the maximum likelihood method in MEGA-X (Kumar, Stecher, Li, Knyaz, & Tamura, 2018). The appropriate model of sequence evolution, including model parameters, was calculated using corrected Akaike Information Criterion and resulted in T92 + G+I (Tamura 3-parameter using a discrete Gamma distribution plus assuming that a certain fraction of sites is evolutionarily invariable) (Tamura, 1992) as the best model for MvitPBP1. The best model for MvitPBP2 was K2 (Kimura 2-parameter)+G + I, whereas K2 + G was selected for MvitPBP3. The models were also selected based on partitioning by codon position. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Joining and BioNJ algorithms to a matrix of pairwise distances estimated using the maximum composite likelihood approach, and then selecting the topology with superior log likelihood value. The bootstrap consensus tree inferred from 1,000 replicates (Felsenstein, 1985) was taken to represent the evolutionary history. Branches corresponding to partitions reproduced in less than 50% of the bootstrap replicates were collapsed. The percentage of replicate trees in which the samples clustered together in the bootstrap test is shown next to the branches (Felsenstein, 1985). The phylogenetic trees were rooted by the outgroup Conogethes punctiferalis.
The primary species hypothesis was evaluated using Automatic Barcode Gap Discovery (ABGD), a molecular species delineation method. ABGD is an automated procedure that clusters sequences into candidate species based on pairwise distances by detecting differences between intra-and interspecific variation without a priori species hypothesis (Puillandre, Lambert, Brouillet, & Achaz, 2012).
The program requires a prior limit to intraspecific diversity (P) and a proxy for minimum gap width (X  (Clement, Posada, & Crandall, 2000).

| Structure of M. vitrata PBP genes
The assembly of the candidate homologs from the transcriptome intron 2 was longer than intron 1, whereas intron 1 was longer than intron 2 in MvitPBP1. and Lys residues, while MvitPBP3 contains more Glu, Ala, and Val residues than other amino acids. The amino acid sequence analysis of MvitPBPs revealed that they consisted of seven α-helices and a conserved motif of six cysteine residues. The location of the α-helices has been predicted following Sandler, Nikonova, Leal, and Clardy (2000) to be located between residues 1-13 (α1a) MvitPBP1 exhibited the highest similarity with CpunPBP1 of C. punctiferalis (87%), followed by CmedPBP1 of C. medinalis (86%).
MvitPBP1 amino acid sequence also did not vary much with the PBPs

| PBPs haplotype variation in M. vitrata population and neutrality tests
The haplotypes identified in Maruca individuals were deposited in the NCBI GenBank (MvitPBP1: MK548942-MK549033, MvitPBP2: South America, which was the lowest (Table 4). However, it was the highest based on MvitPBP3 for South America and it was almost similar for all other continents (Table 4). On continental basis as well, the haplotype diversity was one or close to one for most of the sampled continents for all the MvitPBP genes (

| F-statistics (F ST ) and analysis of molecular variance
The

| Phylogenetic pattern based on MvitPBPs
The   There were only two groups for MvitPBP3 (Table 16). As showed in the phylogenetic tree, ABGD analysis for MvitPBP3 also suggests a single group for Maruca populations from Asia, Africa, and Oceania, whereas South America populations formed a separate group.

| Haplotype network
The haplotype network analysis involving the active MvitPBP1 haplotypes in this study revealed two distinct groups (Figure 9a).
Although the phylogenetic tree and the ABGD grouping clearly differentiated the South America Maruca populations from rest of the populations, they were placed at the periphery of the radial expansion of the major cluster that contained the Asian and Oceania populations in the haplotype network. Surprisingly, few African populations also aligned within this cluster. However, majority of the also obtained for the network based on active MvitPBP2 haplotypes ( Figure 9b). However, the results from the network based on active MvitPBP3 haplotypes clearly differentiated the populations in this study into three clusters-Asia and Oceania as the major cluster, Africa and South America as the two other minor, but distinct clusters ( Figure 9c).

| D ISCUSS I ON
The Another PBP was identified from M. vitrata female adults and named as MvitPBP3. It is common to have more than one PBP in moth species. Earlier studies also reported the occurrence of multiple PBPs in moths that produce multi-component sex pheromones, and each PBP may be encoded by a distinct locus (Newcomb, Sirey, Rassam, & Greenwood, 2002). Hence, it is possible that each PBP recognizes a specific compound in the multi-component pheromone blend. For instance, two PBPs were described in Lymantria dispar (Vogt, Köhne, Dubnau, & Prestwich, 1989), which selectively bound the two pheromone enantiomers (Bette et al., 2002;Du & Prestwich, 1995;Plettner, Lazar, Prestwich, & Prestwich, 2000). Although one of the three PBPs from Antheraea polyphemus (ApolPBP1) was shown to bind to all three pheromone compounds with high affinity at high pH, competitive assays showed considerable differences in affinity only for the major compound (Leal, Chen, & Erickson, 2005 (Breer, Krieger, & Raming, 1990). In addition, the amino acid sequences of the third exon in PBPs should possess three conserved cysteine amino acids (Willett & Harrison, 1999). Hence, the identified MvitPBPs are of the expected size with the presence of six highly conserved cysteine residues.
Pheromone-binding proteins have six α-helices with the pheromone ligand bound in an internal hydrophobic pocket (Sandler et al., 2000). Subsequent studies revealed a seventh α-helix, formed from the C-terminal tail (Horst et al., 2001). We ascertained the location of seven α-helices in MvitPBPs by aligning their amino acid sequences with PBPs and OBPs from Bombycidae, Saturniidae, Sphingidae, and Noctuidae (Malini, 2017). Interestingly, these locations were almost similar to the seven α-helices identified for B. mori and L. dispar (Yu et al., 2012). The three disulfide bonds in MvitPBPs are the same as the two that attach α3 to helices α1 and α6 (Cys19-Cys54 and Cys50-Cys108, but Cys50-Cys109 for MvitPBP3), and the third disulfide bond (Cys97-Cys117 but Cys98-Cys118 for MvitPBP3) connecting helices α5 and α6 reported in B. mori (Sandler et al., 2000). Met74 in α4 and Ile91 in α5 of B. mori PBP were substituted by Gln74 and Val91, respectively, in MvitPBP1. Although Met74 was not substituted by another hydrophobic amino acid, Ile91 was substituted by the hydrophobic amino acid. The amino acids of helices α5 and α6 used to form a hydrophobic assembly in B. mori (Sandler et al., 2000) are the same in MvitPBP1 except Ile93 in α5, which was replaced by hydrophobic Leu93. In other small interhelix contacts, especially between helices α2 and α3, three substitutions (Val48Thr, Leu52Ile and Met55Leu) were found in MvitPBP1. A loop formed by amino acid residues 60-69 is the most flexible region of the protein, and it serves as the lid into the pheromone-binding pocket (Nemoto, Uebayasi, & Komeiji, 2002). Thus, the identified MvitPBPs are similar to the structure of already reported lepidopteran PBPs or GOBPs.
Although structural modeling was used to predict the "presumed" structures of MvitPBPs (Mao et al., 2016), future studies should confirm their three-dimensional structures by X-ray diffraction and/or NMR spectroscopy.   (Sandler et al., 2000). It should also be noted that most of the residues lining the binding pockets were hydrophobic. However, hydrophilic residues, such as threonine present in the binding sites, are probably responsible for hydrogen bonding with the functional group of the ligand (Yu et al., 2012).  (Newcomb et al., 2002;Prestwich & Du, 1997;Rogers, Sun, Lerner, & Vogt, 1997). Thus, this warrants further detailed studies to understand whether these polymorphisms contribute toward the reported differential pheromone recognition patterns of male M. vitrata moths in different geographical regions (Downham et al., 2004;Schläger et al., 2012;Srinivasan et al., 2015).