Molecular Analysis of Targeted Insecticide Resistance Gene Mutations in eld Aedes Aegypti and Culex Quinquefasciatus Mosquito Vectors of Arboviruses from Saudi Arabia

Background In absence of effective and safe vaccines or drugs against mosquito-borne arboviral diseases (such as dengue, Zika and Rift Valley and West Nile), there have been an increasing pattern of insecticide use against the mosquito vectors of these diseases in integrated vector management (IVM) programs in Saudi Arabia. However, the ecacy of IVM programmes is threatened by the development of insecticide resistance in vectors. Gene mutations on target sites can be a valuable reference to the status of insecticide resistance. Jeddah, a global commercial and major port-of-entry city, is bearing the most (> 70%) dengue disease burden and the population of the dengue vector Aedes aegypti in Saudi Arabia. Culex quinquefasciatus is a second example as one of the major arboviral vectors in the region and a proven vector of Rift Valley fever and wide-spread in Jeddah and the rest of the western region of Saudi Arabia. However, the status of insecticide resistance and targeted site mutations on the responsible genes are not fully characterized. Methods We randomly sampled both mosquito species, Ae. aegypti and Cx. quinquefasciatus across Jeddah by daily mosquito surveillance in 2016 to detect the resistance-associated target site mutations on the voltage gated sodium channel (VGSC) and acetylcholinesterase 1 (ace-1) genes. Results Our ndings showed that Ae. aegypti resistance-associated VGSC gene mutations revealed polymorphic mutations on the 989 (allele frequency = 65.91%), 1016 (allele frequency = 65.91%), and 1534 (allele frequency = 52.27%) sites. Additionally, we documented two types of introns between exons 20 and 21, however, the I1011M point mutation was undetected. Linkage disequilibrium associations were shown between V1016G with S989P, V1016G with F1534, and V1016G with type A intron. Furthermore, no mutation on ace-1 was identied in Ae. aegypti. In Cx. quinquefasciatus, homozygous L1014F/L1014F (95.23%) on the VGSC and heterozygous G119/G119S (100%) on ace-1 were widely distributed in the samples studied. Conclusions To the best of our knowledge, this is the rst record of the intron types between exons 20 and 21 on VGSC of Ae. aegypti populations, and the rst report of insecticide resistance gene mutation being present in eld caught Cx. quinquefasciatus in Jeddah, Saudi Arabia. High prevalence of insecticide resistance gene mutations in local primary mosquito vector species, implies that the local Ae. aegypti mosquitoes are resistant to both pyrethroid and DDT, while the resistance to all four insecticide groups may possessed in Cx. quinquefasciatus populations. Our ndings alert the urgent need to carry out a comprehensive insecticide resistance gene mutation surveillance and monitoring to guide sustained and effective IVM planning and innovative guidelines in local predominate mosquito populations and mosquito-borne disease control and elimination in Saudi Arabia. The target resistance-associated genes successfully amplied from a total of 20 Ae. aegypti and 21 Cx. quinquefasciatus. The sequences originally obtained in this study have been deposited in the GenBank under the accession numbers MN997310–MN997399. Three non-synonymous VGSC mutation alleles, S989P, V1016G and F1534C were identied from specimens of wild Ae. aegypti. The presence of a silent mutation (TTG to TTA) at position (Leu, L) 982 (14/20) was found. No I1011M and F1552C mutations were identied. In addition, two types of introns between exons 20 and 21, type A (250 bp) and type B (234 bp) were identied. Six genotypes of the combination 989, 1016, and 1534 positions were documented including SVF/PGF, SVF/SVC, PGF/PGF, PGF/PGC, PGC/PGC Among them, PGF/PGC was the dominant genotype, detected in 45% of the examined samples, by the SVC/SVC with 25%. Genotypes of SVF/PGF, SVF/SVC, PGC/PGC observed less frequently. Allele frequencies The V1016G mutation be more widely distributed with high frequency (67.50%) in Jeddah, that the F1534C mutation Both TTG982TTA (χ 2 = 22.000, df = P S989P (χ 2 df epidemiological IVM and pesticide choices implementation amongst all stakeholders. Furthermore, strengthening laboratory and epidemiologic capacity building programs for health workers, health professional, and community engagement including and private-public sectors partnerships on mosquito vector, dengue and emerging threat and epidemics sustained control towards elimination.


Introduction
Mosquitoes-borne diseases, including dengue, malaria, West Nile, Japanese encephalitis, Rift Valley fever, Chikungunya, Zika, and other emerging arboviral diseases are prevalent and are enormous threat to public health worldwide. Intensive international trade and travel, in uenced by global economic integration, accelerated the spread of mosquito-borne viruses, resulting in occasional disease outbreaks even in previously endemic or disease-free areas. Since no protective vaccines and effective drugs are available for majority of these infectious diseases, the strategies on reducing the risk of transmission is mainly focused on mosquito vector control. Vector control has been relying heavily on insecticides application to reduce mosquito larval and adult populations in their natural habitats as well as personal protection using anti-mosquito repellents [1]. Selection pressure on insecticide target sites of mosquito vectors, resulting in the occurrence and developing of single or multiple-site mutations in the responsible genes have been reported and subsequently detected globally in major mosquito vector species belonging to the genera Anopheles, Aedes and Culex in the last few decades [2,3]. Increased insecticide target site insensitivity due to gene point mutations, and the increased or elevated activity of insecticide detoxifying enzymes are two of the major mechanisms evolved and developed by mosquito vectors for surviving the lethal effect of insecticides [4,5].
The most prevalent insecticide resistance gene mutation is the knockdown resistance (kdr) mutation on the para voltage sodium channel gene (VGSC). The mutation sites on the amino acid 1014 position on VGSC IIS for both Anopheles and Culex, and 1016 position for Aedes, were well con rmed to reduce the mosquito sensitivity to pyrethroid, and the organochlorine DDT due to cross-resistance effect [21][22][23]. In speci c, several point mutation sites on VGSC are speculated to get involved in pyrethroid and DDT resistance in Ae. aegypti, including G923V, L982W, S989P, M1011I/G, V1016G, F1269C in domain II, and T1520I, F1552C, F1534C, D1763Y in domain III of VGSC [1,[24][25][26]. Among them, two point mutations, V1016G and F1534C that confer resistance to pyrethroids and DDT, are the most widespread [27][28][29][30][31][32]. Four non-synonymous (A99S, L1014F/S/C, V1016G, W1594R) and six synonymous, (L884L, G923G, A1262A, D1266D, P1270P, G1754G) mutations were found to be linked with the levels of permethrin resistance in Cx. quinquefasciatus populations [33][34][35]. In addition, the replacement of the amino glycine (Gly, G) with serine (Ser, S) on position 119 of acetylcholinesterase 1 (ace-1) was found to be associated with mosquito resistance to organophosphate and carbamate insecticide groups [36]. The G119S mutation was commonly detected in Anopheles and Culex, while seldomly recorded in Aedes species [1]. In Saudi Arabia, S989P, V1016G and F1534C mutations on VGSC were reported in Ae. aegypti samples collected from Jeddah and Makkah in 2015 [18]. It is the rst on insecticide resistance gene detection in Ae. aegypti populations from Saudi Arabia.
However, little is documented on the trend and impact of insecticides application in integrated vector management (IVM) programs implementation nearly for 20 years (1995-2016) on vector-borne disease threats and epidemics prevention and control in the Kingdom of Saudi Arabia [16,17].
Additionally, point mutations in the VGSC gene of Ae. aegypti populations in Jeddah has not been fully characterized, especially the types of introns between exons 20 and 21, as well as whether the ace-1 mutations in Jeddah Ae. aegypti populations. In Cx. quinquefasciatus, the V1016G mutation was detected in the JPP-B strain, originated from Jeddah and selected by permethrin for 20 consecutive generations [20,34]. This mutation has not been reported in other populations of Cx. pipiens complex [37]. To the best of our knowledge since the rst studies performed in the mid of 1980s, no study on insecticide resistance gene(s) in Cx. quinquefasciatus populations have been reported in the so far in Saudi Arabia. Thus, currently, there are limited information on insecticide resistance gene mutations in Cx. quinquefasciatus against chemical insecticides in Jeddah. For the bioassay test requires numerous live mosquitoes and time-consuming [38]. Early detection and characterization of target mutations of insecticide resistance genes, as potentially valuable molecular markers, can provide a reference to estimate the local insecticide resistance level for monitoring, planning and evaluation of vector control strategy [1,39,40].
This study aimed to determine the status of insecticide resistance in Ae. aegypti and Cx. quinquefasciatus mosquito vectors, by mapping the point mutations associated with kdr mutations and the phylogeny of the introns between the exons 20 and 21 on VGSC, and ace-1. The results of our study will be very important for the early-detection of resistance development and in guiding IVM and insecticide application policies for the prevention and control of mosquito-borne diseases in Saudi Arabia.

Mosquito Collection
A total of 25 Ae. aegypti and 24 Cx. quinquefasciatus female mosquitoes were analyzed, representing different districts and municipalities of Jeddah Governorate. These mosquitoes were collected during routine mosquito surveillance in 2016. Collection sites for detection of target site mutations on insecticide resistance genes are shown in Fig. 1, and a map was generated using ArcGIS 10.1 ArcMap software (ESRI, Redlands, CA, USA).
Taxonomic identi cation and ecological information on the analyzed samples are mentioned elsewhere [41].
DNA extraction and target gene mutation sequencing For DNA extraction, individual mosquito samples, each was homogenized in a mixer mill (Jingxin, Shanghai, China) in 250 µl of ATL (Qiagen, Hilden, Germany) and 20 µl of proteinase K (Qiagen, Hilden, Germany) with two 3-mm steel balls were added to each tube. The homogenates were incubated at 56 °C overnight in an oscillating thermo-block. The homogenate tubes were then centrifuged at 1000 × rpm for 3 min at room temperature. A For Ae. aegypti, the primer set IIS5-6(3)F and IIS5-6(3)R were used to amplify the IIS6 domain of the VGSC, spanning the codons 989, 1011, 1016, and the introns between exons 20 and 21; while the codon 1534 domain IIIS6 was ampli ed by the primer set AaNa31F and AaNa31R [42]. The ace-1 gene region spanning the position 119 was ampli ed by the primer set Ace-1F and Ace-1R [43]. For Cx. quinquefasciatus, the VGSC region spanning the 1014 codon was ampli ed by the primer set Cq1 and Cq2 [44]. The fragment containing the ace-1 mutation was ampli ed by the primer set Moustdir and Moustrev [45].
The PCR ampli cation products were separated on 1.5% agarose gel electrophoresis, puri ed, and sequenced in both directions by Sangon Biotech (Shanghai, China). The peak spectrum (chromatogram) of the sequencing data was determined using Bio-edit v.7.2.6 [46]. The allelic information for certain ambiguous sequences, cloning and plasmid sequencing were performed for con rmation.

Phylogenetic analysis
The obtained VGSC IIS6 fragments spanning the introns anked by exons 20 and 21 were aligned with available sequences retrieved from GenBank using ClustalW2 [47] with default settings, which were manually adjusted when was necessary. Both intra-and inter-genotypes (type A and type B) of introns spanning exons 20 and 21 divergences were examined using the Kimura's 2-parameter (K2P) distance model [48] in MEGA v.7.0 [49]. Based on the Akaike information criterion, the best-t model for the alignment was determined using Model test 3.7 [50], in cooperation with PAUP* v.4.0b10 [51]. Consequently, calculation of the maximum likelihood (ML) and Bayesian likelihood trees were completed under the TVM + G model for introns between exons 20 and 21 on the VGSC. Both the ML and neighbour-joining trees were constructed using MEGA v.7.0 software, with 1,000 bootstrap replicates. The Bayesian tree was built with MrBayes v3.2.6 [52] on the CIPRES portal (www.phylo.org/) [53], run for 10 million generations, with the rst 25% generations discarded as burn-in. The trees were visualised using Figtree v1.4.2 (http://tree.bio.ed.ac.uk/software/ gtree/). Both intra-and interspeci c divergences of introns spanning between exons 20 and 21 were calculated using K2P distance model [48] in MEGA v.7.0.

Statistical analysis
Pearson's Chi-squared test was applied to analyze associations among mutation codons and type of introns between exons 20 and 21 on VGSC. The tests were performed in SPSS v.20 (IBM Corp., Armonk, NY, USA).

Results
Mapping of the target resistance-associated genes in Ae. aegypti and Cx. quinquefasciatus  (Table 1). Among them, PGF/PGC was the dominant genotype, detected in 45% of the examined samples, followed by the SVC/SVC genotype with 25%. Genotypes of SVF/PGF, SVF/SVC, and PGC/PGC were observed less frequently. Allele frequencies are shown in Table 2. The V1016G mutation seems to be more widely distributed with high frequency (67.50%) in Jeddah, than that of the F1534C mutation (55.00%). Both TTG982TTA (χ 2 = 22.000, df = 2, P < 0.001), and S989P (χ 2 = 44.000, df = 4, P < 0.001) co-occur with the presence of V1016G. All TTG982TTA, and S989P coexist with the haplotype harboring the V1016G mutation. In individuals carrying homozygous F1534C mutation in the IIIS6 region, no coexistence of other mutations on 982, 989, and 1016 codons in IIS6 domain was observed. The type A intron between exons 20 and 21 was strongly linked to the V1016G mutation.  Phylogenetic analysis of intron types identi ed in the voltagegated sodium channel gene ( VSGC ) of Aedes aegypti The intron located between exons 20 and 21 on VGSC identi ed in our study and those retrieved from the GenBank were compared and used to calculate the evolution of the kdr mutations. The studied 158 sequences can be clearly divided into two clades with high bootstrap values, one clade contained the sequences possessing type A intron, while the second clade contained those sequences having type B intron (Fig. 2) position (χ 2 = 49.910, df = 4, P < 0.001) than to the 1534 position (χ 2 = 16.297, df = 4, P = 0.003). While the intron polymorphism was more correlated to the 1016 position (χ 2 = 34.333, df = 2, P < 0.001) than the 1534 position (χ 2 = 18.005, df = 2, P < 0.001) ( Table 3).  [3,18,21,33,36,37].
No samples that harbored the wild-type at the three mutation loci simultaneously were identi ed. The I1011M/V mutation was absent in Jeddah, even though this substitution of isoleucine on 1011 position is prevalent in Latin America countries [39,42,[54][55][56], and spread in Thailand and Vietnam but at low frequencies [57,58]. Yet, the situation has persisted with increasing urbanization, climate changes and sustained Aedes-Culex endemicity estimated at 31-56% prevalence rate. The mitigation and slowing down the spread and impact of insecticide resistance development in disease vector populations, needs integrated laboratory assessment coupled with surveillance, monitoring and evaluation to better understand and obtain the necessary scienti c knowledgebased guidance and informing pesticides policies on selection, purchase and the long-term and sustainable vector control efforts.
So far, to our knowledge, Saudi Arabia is the place with the most variation in kdr mutations associated with high resistance to pyrethroids in the Middle East countries, probably due to long-term usage of high doses of pyrethroids in vector control programs in and all major cities, especially in Jeddah since the rst dengue epidemics in 1994. This is different from the reported kdr polymorphisms of Ae. aegypti populations from India and Pakistan. High prevalence of the F1534C mutation, and a low frequency of I1011M/V and V1016I mutations, with the absence of the 989P allele, were observed in India [25]. In Pakistan, all the studied mosquito samples harbored the homozygous F1534C mutation, whereas no mutation on the 1016 position has been identi ed yet [59].
The percentage of homozygous V1016G and S989P mutations co-occurrence was high (65.0%) in Jeddah. The V1016G mutation plays a critical role in kdr resistance against pyrethroid compounds [22,23,59]. Linkage association has been observed between the 989 and 1016 loci in Ae. aegypti populations of Saudi Arabia [18], Myanmar [60], Thailand [61]. The occurrence of S989P was highly linked to the V1016G mutation as a compensatory mutation to reduce the tness cost on the mosquito [22] and strengthen the response of V1016G to permethrin and deltamethrin [62].
Additionally, double homozygous mutations at positions 989 and 1016 combined with either wild-type or heterozygous mutation at 1534 position were common in Jeddah's population, in agreement with the speculation that the 1016 and 1534 loci were not independent for linkage disequilibrium analysis [63]. Further, the genotype with the triple mutation haplotype predominated in the tested mosquito samples, with an occupation of 45% (PGC/PGF). Additionally, one individual possessing triple homozygous mutation on 989, 1016, and 1534 loci was found in this study. In eld-collected mosquitoes, the triple homozygous mutation was in low frequency in the population [32,59,60,64,65], probably due to tness cost exerted on the mosquito by this triple mutation [62,66]. Whereas, the arti cially-introduced triple mutation, using the Xenopus oocyte expression systems, appeared to confer a higher level of resistance to both permethrin and deltamethrin than those carrying the individual mutations, probably for synergistic effect of the combination of mutant alleles [22,23,62]. Li et al., [30] suggested that the F1534C mutation could possibly serve as a compensatory mutation for reducing the tness cost on the mosquito induced by the V1016G mutation. Thus, whether the homozygous triple mutations on 989, 1016, and 1534 loci showing insecticide resistance in eld mosquitoes is yet to be ascertained.
The intron polymorphism of VGSC could serve as a marker to study the evolution of kdr mutations [38,54]. Based on sequence size, the introns located between exons 20 and 21 were classi ed into type A (250 pb) and type B (234 pb) [54]. Analysis of intron sequences obtained in our study and homologous sequences retrieved from GenBank showed that the intron type was signi cantly associated with 1016 (χ 2 = 34.333, df = 2, P = 0.000), found to be strictly associated with the type B intron, and the V1016G mutation coexisted with the type A intron [38]. By contrast, in Africa, the point mutation at 1534 site was found to be strongly linked with the type A intron in Ghana, but only one heterozygote point mutation V1016I was recorded, except the F1534C mutations [67]. Additionally, V1016I was reported to locate at the haplotype possessing the type A intron, with no 1534 variation observed in Brazil [54,68]. The mutation at position 1016 is associated with resistance to both type I and II pyrethroids, while the F1534C allele is primarily associated with resistance to type I pyrethroids, and the V1016I mutation did not alter channel sensitivity to pyrethroids in Xenopus oocytes [22,23]. Based on the evidences above, we speculate that kdr mutations incline to coexist with the type A intron. When a mutation occurred in either 1016 or 1534 positions, it would be linked to the type A intron. While in case of co-mutation, the type A intron distributes more frequently with the V1016G/I mutation rather than the mutation F1534C. It is probably due to reducing the tness cost on the mosquito by increasing compensatory advantage or linkage disequilibrium [62,69].
It is speculated that the type B intron is the ancestral clade of clade 1, because that the Ae. aegypti was originally from Africa and the majority of African individuals belong to type B [67] [70]. In view of the phylogenetic tree obtained with available global data of the intron sequences, 100%, 83.33% and 70.40% individuals from North America, Latin America and Asian countries respectively, were type A intron. By contrast, 73.91% of sequences from Africa belonged to the type B intron. According to our analyses above, the type A intron coexistence with kdr mutations is common, especially the V1016G mutation; while the type B couples with the wild-type or F1534C mutation when the V1016G mutation is absent. These variable pattern of association, may re ect the history of insecticide treatment in different continents. The major disease burden of malaria is in Africa and South-East Asia [71]. The mutation on 1534 locus may have arisen in Aedes species due to earlier or current use of DDT for larvicides spraying and IRS to control malaria since the middle of the last century [42,72].
Whereas both in Latin and North America countries, vector control is mainly focused on Aedes species (mainly Ae. aegypti and Ae. albopictus) vectors of dengue, Chikungunya, yellow fever, and Zika [73]. For the domestic and endophagic habitats of Ae. aegypti, pyrethroid insecticides are widely and randomly used in household aerosol insecticides for public and personal protection, since that pyrethroids have the advantages of high e cacy against mosquito vectors, low mammalian toxicity and short residual action [74]. In South Asia countries, the higher frequency of type A intron than that of type B may indicate different selective pressures resulting in different histories of intense and unrestricted insecticide usage. It is probably that rst, the F1534C mutations occurred due to the extensive use of DDT or type I pyrethroids, as both insecticides have the same target site [42], then V1016G has been selected by the subsequent and wide use of type II pyrethroids, as cypermethrin and deltamethrin [60].
Reports on target-site mutations on VGSC and ace-1 genes of Cx. quinquefasciatus are rare compared of that on Ae. aegypti. However, the insecticide resistance level of the Cx. quinquefasciatus in Jeddah may even worse as inferred from the xation of mutation allele both on VGSC and ace-1 genes. In the studied mosquito samples in Jeddah, the kdr homozygous genotype TTT/TTT occupied most (90.48%) of the individuals. While the heterozygous genotype GGC/AGC on ace-1 gene was detected in all samples of Cx. quinquefasciatus. This is probably related to constant exposure due the extensive use of insecticides for routine spraying campaigns of public health, or for domestic use against mosquito vectors and other insect pests in Jeddah region. Similarity, the homozygous susceptible 1014 genotype was absent in samples collected in 2011 from Zanzibar, East Africa [75]. In India, most of the studied mosquito samples harbored the heterozygous genotype TTT/TTA, with the frequency of 75% in Cx. quinquefasciatus samples [44].
Both the non-synonymous mutation G119S (GGC/AGC) and silent mutation T506T on ace-1 [43] were not detected in our sequences in Ae. aegypti samples from Jeddah. While it does not indicate that the Ae. aegypti population was completely susceptible to organophosphate and carbamate compounds, as both types of insecticides have the same target site. In bioassay tests, the Ae. aegypti population from Jeddah was tolerant to the fenitrothion [18], indicating that other than kdr mutations, or more likely, enzymatic mechanisms are involved in conferring resistance in mosquitoes. In Central Java, Indonesia, although the G119S mutation was absent, a degree of resistance to malathion was observed in Ae. aegypti [1].

Limitations of this study
One major limitation of our study is the small sample sizes of mosquitoes used for detection on target site mutations of insecticide resistance genes. Whereas, the primary survey reveals the widespread of kdr in predominant mosquitoes in Jeddah. Further, bioassay tests on both larvae and adult mosquitoes are required to evaluate the status of insecticide resistance in practice. And the underlying metabolic mechanism responsible for insecticide resistance is worth exploring.

Conclusions
The current study involved a wider geographic region and covering a larger sample of Ae. aegypti vector populations in Jeddah. The results showed that six genotypes of the three codon combinations, 989, 1016, and 1534 were spreading in Jeddah's Ae. aegypti population, with the PGF/PGC as the dominate one. Two types of introns between exon 20 and 21 on VGSC for the rst time are identi ed in Saudi Arabia. Linkage associations have been shown between V1016G with S989P, V1016G with F1534, and V1016G with type A intron. High prevalence of kdr mutations implies that the local Ae. aegypti mosquitoes are pyrethroid-resistant, and the level of resistance may increase with the extensive and intensive use of pyrethroid, as the active ingredients in household insecticide and spatial repellents. To the best of our knowledge, this is the rst report of resistantassociated mutations in eld caught Cx. quinquefasciatus in Saudi Arabia. The xation of L1014F allele on VGSC and G119S on ace-1 gene was detected in local Cx. quinquefasciatus populations, with the frequency of 95.24% and 100%, demonstrating that traditional insecticides are likely to be losing e cacy through fogging activities or IRS.
The primary analyses showed the widespread of kdr mutation both in Ae. aegypti and Cx. quinquefasciatus, thus threating the e ciency of vector control against mosquito-borne diseases. It raises the attention that widerange of surveillance the polymorphism and distribution of insecticide target site mutations in mosquitoes in this region. The mapping of the nation-wide distribution and polymorphism of target site mutation against insecticide resistance will provide valuable consensus for forecasting the status of insecticide resistance in mosquito populations of Jeddah and other regions of Saudi Arabia. Additionally, regular monitoring and documenting the dynamics of target mutation gene and their associations with the insecticide resistance phenotype is important in the surveillance system of preventing mosquito-borne diseases through effective vector control strategies.
Meanwhile, the above ndings suggest that alternative insecticides in areas, and proper use of insect growth regulator or biocontrol approaches for integrated vector control programs should be considered in Jeddah and neighborhoods. Albeit, substantial Saudi government and Jeddah municipality commit and invest in innovative operational research, research and development on more early-detection of insecticide resistance and effective control measures and tools for discovery and utilization, including vaccines for local and pilgrims' immunization programs implementation. With the establishment of the Public Health Pest Laboratory of Jeddah Governorate, sustainable funding is imperative to provide evidence for coherent and coordinated temporal and spatial epidemiological IVM and pesticide choices implementation amongst all stakeholders. Furthermore, strengthening laboratory and epidemiologic capacity building programs for health workers, health professional, and community engagement including and private-public sectors partnerships on mosquito vector, dengue and emerging threat and epidemics sustained control towards elimination. No speci c permits were required for this study. The study did not involve endangered or protected species.

Consent for publication
Not applicable.

Availability of data and materials
All data generated or analyzed during this study are included in this published article. The newly generated sequences were submitted to the GenBank database under the accession numbers MN997310-MN997399.

Competing interests
The authors declare that they have no competing interests.

Funding
The Special Foundation of Basic Science and Technology Resources Survey of Ministry of Science and Technology of China (grant no. 2017FY101200).

Author contributions
The paper was conceived Tambo