Multi-gene Phylogenetic Analysis and Genetic Diversity of Discrete Elytral Color Phenotypes in Menochilus sexmaculatus (Coleoptera: Coccinellidae)

Weidong Huang1,2,3, Xinyue Liang1,2,3, Xiufeng Xie4, Xingmin Wang2, 3 and Xiaosheng Chen1,2,3* 1Department of Forest Protection, College of Forestry and Landscape Architecture, South China Agricultural University, Guangzhou 510640, China. 2Key Laboratory of Bio-Pesticide Innovaation and Application, Guangdong Province, Guangzhou 510642, China 3Engineering Research Center of Biocontrol, Ministry of Education and Guangdong Province, Guangzhou 510640, China. 4Guangdong Agriculture Industry Business Polytechnic College, Guangzhou 510507, China. Article Information Received 10 December 2019 Revised 23 March 2020 Accepted 25 March 2021 Available online 21 October 2021


INTRODUCTION
I nsect body coloration often shows genetically different forms even within a population, which may be important for intra-and interspecific communication and adaption to the local environment (Noriyuki and Osawa, 2015). Color polymorphisms provide some of the best characterized examples of functionally and ecologically important polymorphism. Insect melanism describes the occurrence of varied pigment patterns both within and between closely linked species making them polymorphic (True, 2003). During the past two decades, the genetics underlying variation in melanism have been unraveled for several O n l i n e

F i r s t A r t i c l e
ladybird H. axyridis is an emblematic species of elytral color pattern polymorphism, with more than 200 distinct color patterns described around the world (Ando et al., 2018;Gautier et al., 2018), and have brought difficulties in their identification or nomenclature. Intraspecific variation in color patterns of ladybird beetles has been studied for many decades (Millar et al., 1999). Recently, the phylogenetic relationships of different color patterns within H. axyridis was investigated based on 12S rRNA and 16S rRNA genes (Yao et al., 2011). The color polymorphism of this species was controlled by the transcription factor pannier (Gautier et al., 2018). Menochilus sexmaculatus is a polymorphic aphidophagous ladybird with wide prey range and distribution (Kawakami et al., 2013). However, the broad distribution and multiple color forms in M. sexmaculatus were usually confused for their identification (Fig. 1). This species is known to have 20 phenotypes of elytral polymorphs according to the ratios of elytral red and black areas (Kawakami et al., 2013(Kawakami et al., , 2018. Meanwhile, the phylogenetic relationships and genetic diversity of different color patterns within M. sexmaculatus is still poorly known. Recently, polymorphism has become the focus of attention in population ecology as well as evolutionary biology as it can contribute to population productivity, stability and persistence (Wennersten and Forsman, 2012;Forsman, 2013;Takahashi et al., 2014). Here we investigate the molecular phylogenetic analyses and genetic diversity of different elytral color forms in M. sexmaculatus based on two mitochondrial genes (cytochrome oxidase subunit I (COI) and II (COII)) and two nuclear genes (carbamoyl phosphate synthetase (CAD) and histone subunit 3 (H3)). The aim of this study is testing the different color patterns of M. sexmaculatus whether possess phylogenetically specific relationships.

Sampling
Multiple elytral color forms in M. sexmaculatus ( Fig.  1) were collected from Guangdong Province, China during September 2017. Each phenotype was sampled by two individuals. After collection, the specimens were preserved in absolute ethanol. Detailed sample information is shown in Table I.

DNA extraction and PCR
Total genomic DNA was extracted using the TIANGEN DNA extracting kit (TianGen Biochemistry, Beijing, China) following the manufacturer's instructions. Four gene regions (two mitochondrial genes COI, COII, and two nuclear genes CAD, H3) were amplified and sequenced. The primers information was listed in Table II. Polymerase chain reactions (PCR) were performed in 25 μL volumes containing 12 μL 2 × EasyTaq PCR SuperMix (TransGen Biotech, Beijing, China), 10 μL ultrapure water, 1 μL of each primer and 1 μL DNA template. PCR cycling conditions consisted of initial denaturation at 94℃ for 3 min, 35 cycles of 94℃ for 30 s, 50℃ for 30 s, 72℃ for 1 min, and ending with a final extension at 72℃ for 5 min. For COII and CAD, we used a hemi-nested approach. We performed an initial PCR reaction using following primer pairs: CD806F3/CD1098R2 for CAD, COIIF-leu/ COIIR-lys for COII. One microliter of product from initial reaction was used as template for the hemi-nested reaction using primers CD821F/CD1098R2 for CAD, COIIF-leu/ COIIR-9b for COII. Successful amplification was assessed using gel electrophoresis on 1% agarose gels by adding 5 μL PCR product. All obtained sequences were compared using BLAST against GenBank to ensure that the target sequences were amplified.

Alignment and sequence analyses
Target sequences were manually cleared, trimmed and O n l i n e H3R ATATCCTTRGGCATRATRGTGAC aligned using Geneious 9.1.5 (Kearse et al., 2012), the four gene alignments were then concatenated to obtain a supermatrix using the program Sequence Matrix 1.7.8 (Vaidya et al., 2011). MEGA 7.0 (Kumar et al., 2016) was used to analyze nucleotides composition and pairwise genetic distance based on Kimura-2-parameter (K2P) (Kimura, 1980) models using combined dataset. DnaSP 5.1 (Librado and Rozas, 2009) was performed to calculate the number of polymorphic sites, haplotype diversity (H d ) and nucleotide diversity (P i ) used single gene and combined dataset.

Phylogenetic analyses
Both methods of maximus likelihood (ML) and Bayesian inference (BI) were employed to explore the phylogenetic relationships within different elytra color forms in M. sexmaculatus. Three datasets were assembled for phylogenetic analyses: (1) the P1 matrix, including only mtDNA; (2) the P2 matrix, including only nuclear genes; (3) the P3 matrix, including both mtDNA and nuclear genes. Psyllobora vigintiduopunctata (Linnaeus, 1758) and Halyzia straminea (Hope, 1831) were chosen as outgroups. Together with M. sexmaculatus they belong Phylogenetic Analysis and Genetic Diversity of Color in M. sexmaculatus O n l i n e

F i r s t A r t i c l e
to the tribe Coccinellini according to previous studies (Escalona et al., 2017). Partition Finder 1.1.1 (Lanfear et al., 2012) was used to infer the optimal partition schemes and models of molecular evolution for the concatenated data sets, applying an all search approach with branch lengths unlinked across partitions and the Bayesian information criterion (BIC). The ML analyses were conducted with the program RAxML 8.0 (Stamatakis, 2006). Since it is not currently possible to specify different models of substitution for different partitions in RAxML, we used GAMMA model for each ML analysis. Branch support was estimated with 500 replicates using a rapid bootstrapping algorithm (Stamatakis et al., 2008). The BI analyses were calculated in MrBayes 3.2 (Ronquist et al., 2012). Two Markov Chain Monte Carlo (MCMC) runs were performed with one cold and three heated chains for 30 million generations and sampled every 1000 generations. The consensus tree was estimated after a burn-in of 25% of the sampled trees. The chain stationarity was visualized by plotting likelihoods against the generation number using the program Tracer 1.6 (Rambaut et al., 2014).

Genetic diversity and structure
To gain insight into the genetic diversity of the different elytral forms of M. sexmaculatus, we analyzed the genetic variation of 34 individuals standing for 17 phenotypes of elytral forms, based on 759 bp for CAD, 829 bp for COI, 722 bp for COII and 328 bp for H3, 2638 bp in total. These sequences were deposited in GenBank under the accession number MH589128-MH589261 (Supplementary Table  S1). The average nucleotide contents of A, T, G and C were 35.1%, 25.9%, 22.4% and 16.6% for CAD, 33.7%, 37.9%, 13.8% and 14.6% for COI, 35.3%, 38.9%, 11.2% and 14.6% for COII, 27.7%, 20.8%, 23.8% and 27.7% for H3, and 33.8%, 32.6%, 16.8% and 16.8% for combined dataset, respectively. The nucleotide compositions of CAD, COI and COII were similar, having high A + T content, but H3 gene G + C content were slightly higher than A + T content. In four gene regions, there were 15 variable sites for CAD gene sequences, 9 of which were parsimony informative, 17 variable sites for COI gene sequences, 11 of which were parsimony informative, 11 variable sites for COII gene sequences, 9 of which were parsimony informative, 2 variable sites for H3 gene sequences, 2 of which were parsimony informative, whilst 45 variable sites for combined dataset sequences, 31 of which were parsimony informative.
Based on single gene, the genetic structure of 17 different phenotypes of color forms in M. sexmaculatus was also analyzed. The number of CAD haplotypes was higher, with 13 haplotypes, whilst the number of haplotypes for H3 was least, only 3 haplotypes were found. The haplotype diversity (h) was remarkably high (0.731) for CAD, and the nucleotide diversity (π) of COI was highest (0.007) (Table  III). Besides, the combined dataset was used to analyse the genetic structure of the different phenotypes of elytral forms in M. sexmaculatus. Results show that the number of haplotypes was 17, the haplotype diversity was 0.902 and the nucleotide diversity was 0.004. For the haplotypes of the combined dataset, eight individuals shared the H4 haplotype and seven individuals shared the H7 haplotype, H1-H3, H5-H6, H8-H10, H12, H14-H17 were represented by only one individual (Table I). However, the same elytral color patterns in M. sexmaculatus did not shared the same haplotype.
The estimated intra-specific genetic distance based on combined dataset among the different elytral forms in M. sexmaculatus ranged from 0 to 0.011 (Supplementary Table S2). The average genetic distance was 0.005, and the maximum genetic distance between M. sexmaculatus-C2 and M. sexmaculatus-H1 was 0.011.

Phylogenetic analyses
The best-fit partition scheme with corresponding substitution models for each dataset was shown in Table  IV. Based on P1 and P3 datasets, the ML and BI topologies respectively resulting from RAxML and MrBayes analyses O n l i n e

F i r s t A r t i c l e
Phylogenetic Analysis and Genetic Diversity of Color in M. sexmaculatus 5 were largely congruent except for several specimens showing the different placement (Figs. 2 and 4). We recover the thirty-four specimens of M. sexmaculatus forming monophyletic clade with strong support in both ML and BI analyses (P1: BS/PP = 100/1; P3: BS/PP = 88/1). Two major clade groups were recognized one clade includes 16 individuals, the second clade includes 18 individuals. However, when P2 dataset was used the resulted topology was different that the tree topologies inferred from P1 and P3 datasets (Fig. 3). This difference is mainly manifested in incongruent basal branching.  Table I. Bootstrap support and posterior probabilities shown near the nodes.

DISCUSSION
Morphological polymorphisms are well studied and provide evidence of natural variation and micro evolutionary processes occurring in nature (Ford, 1964). Ladybird beetles are one of the classical groups in studying the mechanisms that determine local and temporal trends in color polymorphism. Hence among the frequently studied polymorphic taxa ladybird beetles occupy the main position, especially the results of studies on the polymorphism in the pattern and color of their head (Rogers et al., 1971), pronotum (Blehman, 2007) and elytra (Yao et al., 2011).  Table I. Bootstrap support and posterior probabilities are shown near the nodes.
It is important to derive information from molecular research, because they provide direct evidence of interaction between environmental factors and animal characteristics at the organismic level (Lowe et al., 2004). Based on two O n l i n e

F i r s t A r t i c l e
W. Huang et al.
mitochondrial (COI and COII) and two nuclear (CAD and H3) genes, we analyzed the phylogenetic relationships and genetic diversity of seventeen different phenotypes of elytral forms in M. sexmaculatus. Base composition analyses indicated the M. sexmaculatus of mitochondria genes COI and COII with A + T content apparently higher than G + C content, according with the feature of A + T high content in insect mitochondrial genes (Simon et al., 1994) and the average genetic distance based on combined dataset was 0.005 and shows their close relationships. Either single gene or combined data sets consistently revealed the lower levels of nucleotide diversity in different phenotypes in M. sexmaculatus. In phylogenetic analyses of molecular sequence data, choosing an appropriate partitioning scheme is an important step in most analyses due to it can affect the accuracy of phylogenetic reconstruction (Lanfear et al., 2012). Based different datasets resulted in molecular phylogenies analyses showed thirty-four individuals of 17 different elytral forms in M. sexmaculatus formed a clade with high support value, and the two individuals of the same elytral forms do not cluster together indicated the different color patterns of M. sexmculatus do not possess specific relationship phylogenetically. Similar research has been reported by Yao et al. (2011) who conducted phylogenetic analyses based on mitochondrial genes among different elytral forms of H. axyridis. In most cases the occurrence of morphs appears to be associated with climatic factors such as temperature, visual predation (Brakefield, 1985) and industrial pollution (Zakharov, 2003). Dubey et al. (2016) used different temperature regimes to assess the mate choice, reproductive success and offspring coloration of typical and melanic morphs of the M. sexmaculatus. Their findings on offspring phenotype variation indicated that the degree of melanism in morphs is a result of environmentally regulated expression of the parental genotype. In addition to these factors have mentioned above, understanding the genetic mechanisms generating and maintaining such phenotypic variation within species is essential to comprehending morphological diversity. Indeed the genetic and ecological mechanism for the maintenance of elytral color polymorphism in ladybirds is not fully understood (Noriyuki and Osawa, 2015). However, Gautier et al. (2018) combined whole-genome sequencing, population genomic, gene expression and functional analyses, showed that the gene pannier controls melanic pattern polymorphism in H. axyridis. They also pointed out that pannier, which encodes an evolutionary conserved transcription factor, is necessary for the formation of melanic elements on the elytra. Allelic variation in pannier leads to protein expression in distinct domains on the elytra, and thus determines the distinct color patterns in H. axyridis. Meanwhile, Ando et al. (2018) through loss-of-function analyses, genetic association studies, de novo genome assemblies, and gene expression data revealed that repeated inversions within a pannier intron drive diversification of intraspecific color patterns of H. axyridis. These findings provide a reference for our understanding the genetic mechanisms of the elytral color polymorphism in M. sexmaculatus. Therefore, in order to definitely unravel specific causes result in phenotype diversity for M. sexmaculatus, combine genetic data and environmental data is needed in further research.  Table I. Bootstrap support and posterior probabilities are shown near the nodes.

ACKNOWLEGGEMENTS
We are grateful to Xiaoshuang Wang for her guidance on molecular experiments. We are deeply indebted to Dr.

O n l i n e F i r s t A r t i c l e
Phylogenetic Analysis and Genetic Diversity of Color in M. sexmaculatus Shaukat Ali (SCAU), who helped to check the English text. We also would like to express our great appreciation to the anonymous reviewers for their valuable suggestions and comments on our manuscript, which improved this article surely. The present study was supported by the Natural Science Foundation of Guangdong Province

Supplementary material
There is supplementary material associated with this article. Access the material online at: https://dx.doi. org/10.17582/journal.pjz/20210128120158

Statement of conflict of interest
The authors declare no conflict of interest.