Skip to main content
  • Research article
  • Open access
  • Published:

Genomic insight into diet adaptation in the biological control agent Cryptolaemus montrouzieri



The ladybird beetle Cryptolaemus montrouzieri Mulsant, 1853 (Coleoptera, Coccinellidae) is used worldwide as a biological control agent. It is a predator of various mealybug pests, but it also feeds on alternative prey and can be reared on artificial diets. Relatively little is known about the underlying genetic adaptations of its feeding habits.


We report the first high-quality genome sequence for C. montrouzieri. We found that the gene families encoding chemosensors and digestive and detoxifying enzymes among others were significantly expanded or contracted in C. montrouzieri in comparison to published genomes of other beetles. Comparisons of diet-specific larval development, survival and transcriptome profiling demonstrated that differentially expressed genes on unnatural diets as compared to natural prey were enriched in pathways of nutrient metabolism, indicating that the lower performance on the tested diets was caused by nutritional deficiencies. Remarkably, the C. montrouzieri genome also showed a significant expansion in an immune effector gene family. Some of the immune effector genes were dramatically downregulated when larvae were fed unnatural diets.


We suggest that the evolution of genes related to chemosensing, digestion, and detoxification but also immunity might be associated with diet adaptation of an insect predator. These findings help explain why this predatory ladybird has become a successful biological control agent and will enable the optimization of its mass rearing and use in biological control programs.

Peer Review reports


The remarkable evolutionary success of insects is associated with adaptations to a vast diversity of food sources and access to multiple trophic niches. For example, the emergence of gene families encoding odorant binding proteins and odorant receptors allowed insects to locate new diet sources [1]. The expansion of diet range is associated with the expansion of gene families related to detoxification and digestion [2, 3]. In beetles, several studies have demonstrated that the adaptation to plant feeding includes the evolution of genes encoding chemosensors for finding appropriate food sources [4], digestive enzymes for breaking down plant cell walls [5,6,7], and detoxifying enzymes for eliminating harmful plant toxins [7, 8]. In addition, diet also affects insect immunity [9]. The evolution of insect immunity allows insects to change their phenotype in response to changes in the environment, including diet and microbiota [10]. For example, insect antimicrobial peptides (AMPs) can maintain a core microbiota while protecting against microbes [11]. In ladybird beetles, the invasive species Harmonia axyridis (Pallas, 1773) has more genes encoding AMPs than non-invasive species, which might reflect its invasive biology [12]. However, it is not clear whether feeding related traits, e.g. predatory efficiency, prey specialization and adaptability to feeding on unnatural or artificial foods in some beetle species is associated with similar patterns of genomic evolution, while these traits are of particular importance in biological control use.

The use of predatory insects in classical and augmentative biological control programs has yielded in some cases huge economic and ecological benefits. After the successful control of cottony cushion scales using the vedalia ladybird beetle Novius cardinalis (Mulsant, 1850) [13] in California in 1888–1889, hundreds of predatory insects were introduced from abroad for biological control purposes all around the world, but most of them failed to establish or provide pest control [14]. Some species used in classical or augmentative biological control programs even became invasive and harmed local biodiversity [15]. In contrast, the mealybug destroyer Cryptolaemus montrouzieri Mulsant, 1853 is a successful predator and is still being used worldwide [16]. This predatory ladybird beetle is native to Australia and has been introduced to at least 64 countries or regions for classical or augmentative biological control purposes since 1891 [16]. The success of C. montrouzieri can be attributed to its efficient predation of mealybug pests and easy mass rearing [16,17,18].

Mealybugs (Hemiptera, Sternorrhyncha, Pseudococcidae) are the predominant prey of C. montrouzieri. Whereas mealybugs produce wax secretions to protect themselves from a range of natural enemies, these wax secretions act as an attractant and oviposition stimulant for C. montrouzieri [16], indicating ladybird-mealybug specialization. Under laboratory and mass rearing conditions, C. montrouzieri can also feed on other Sternorrhyncha species (e.g., whiteflies, aphids and other coccids), lepidopteran eggs and even artificial diets [18,19,20,21]. Some of these alternative diets can support the complete life cycle of the ladybird (provided that an artificial oviposition substrate is supplied) but will to some extent decrease fitness of the predators. Previously, we detected a large number of differentially expressed genes (DEGs) in C. montrouzieri in response to a diet shift from mealybugs to aphids [21]. This suggests that C. montrouzieri can adapt to a variety of nutritional conditions via phenotypic and transcriptional plasticity.

In this study, we hypothesize that diet adaptation of C. montrouzieri is associated with evolution and regulation of genes related to chemosensing, digestion, detoxification and immunity. We used genomic and transcriptomic approaches to examine the extent of dietary adaptation in C. montrouzieri (Fig. 1). We assembled a high-quality genome of C. montrouzieri and compared its content to eight other Coleoptera genomes. We further tested for gene expression differences between C. montrouzieri larvae that were experimentally fed different diets.

Fig. 1
figure 1

The study design for exploring the genomic basis of diet adaptation of Cryptolaemus montrouzieri


General genomic features of C. montrouzieri

A total of 115.55 Gb of raw data and 106.63 Gb of high-quality clean reads were generated with PromethION DNA sequencing (Oxford Nanopore, UK). These Nanopore data together with additional 151.03 Gb Illumina data were assembled using Wtdbg, followed by Racon and Pilon polishing, which produced a 988.11 Mb genome assembly with 398 contigs and a contig N50 of 9.22 Mb (shortest: 39,165 bp; longest: 32,637,267 bp). This genome size of C. montrouzieri was larger than that of published ladybird and other Coleoptera genomes (largest among the ladybirds, Propylea japonica (Thunberg, 1781), 850.90 Mb; largest among Coleoptera, Anoplophora glabripennis (Motschulsky, 1853), 981.42 Mb) [5, 22]. Application of the Benchmarking Universal Single-Copy Orthologs (BUSCO, Insecta set) pipeline [23] showed that this C. montrouzieri genome compared well with the other insect genomes in the OrthoDB v10.1 database in terms of completeness, with 97.1% complete genes (96.0% single copy and 1.1% duplicated), 0.7% fragmented and 2.2% missing at the genome level.

Annotation of the C. montrouzieri genome using the Braker pipeline [24] yielded a final set of 27,858 genes and 32,187 protein sequences. Application of the BUSCO pipeline showed that this C. montrouzieri gene set has 93.1% complete genes (91.9% single copy and 1.2% duplicated), 4.2% duplicated and 2.7% missing at the protein level in the Insecta of OrthoDB database. In the functional annotation of this protein set, 31,632 were found in the National Center for Biotechnology Information (NCBI) nonredundant (NR) Hexapoda subset, 30,884 in Swiss-Prot, 16,042 in at least one protein domain in Pfam, 7613 in Gene Ontology (GO) and 8290 in Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (Additional file 2: Table S1).

Comparative genomics

A genome-wide scan of gene family evolution was performed among the genomes of C. montrouzieri and eight other Coleoptera species with different feeding habits (Table 1) [4, 5, 25,26,27,28,29]. As revealed by the clustering algorithm implemented in CAFE software [30], we found that 2426 and 2577 gene families of the C. montrouzieri genome underwent expansion and contraction, respectively (Additional file 1: Fig. S1). Among these, only 28 gene families underwent significant contraction (P < 0.05), among which two encode chemosensors and eight encode digestive or detoxifying enzymes (Table 2 and details in Additional file 2: Table S2). Of the 598 significantly expanded gene families in C. montrouzieri (P < 0.05), one encodes chemosensors, nine encode digestive or detoxifying enzymes, and one encodes immune effectors (Table 2 and details in Additional file 2: Table S2). Further identification of immunity-related genes showed that C. montrouzieri had a large number of genes encoding immune effectors including 18 antimicrobial peptides (AMPs) and 15 lysozymes. This number of immune effector genes is equal to those of the beetle Onthophagus taurus (Schreber, 1759) [25] and larger than the other beetles (Fig. 2). In contrast, only 33 genes of C. montrouzieri are involved in recognition, while some of the other beetles have around 60 of these genes (Fig. 2).

Table 1 Genomes of Coleoptera species with different feeding habits used for comparative genomic analyses. Species IDs were ordered based on the species tree topology (Additional file 1: Fig. S1)
Table 2 Significant expansion (E) or contraction (C) of Cryptolaemus montrouzieri in gene families encoding chemosensor, digestive and detoxifying enzymes and immunity as detected with CAFE in comparison to the other beetle genomes. Orthologous genes were identified by Orthofinder. The number of genes in each family are shown for each species. These gene families contain protein domain of: odorant receptors (OR), odorant binding protein (OBP), maltase, glycosyl hydrolase (GH), trypsin, cathepsin, cytochrome P450 (P450), UDP-glucuronosyltransferases (UGT), carboxylesterase (CE) and attacin. Abbreviations of the tested species are defined in Table 1
Fig. 2
figure 2

Evolution of immunity-related genes in Cryptolaemus montrouzieri. Number of genes related to immune recognition and response identified from nine beetle genomes. The species’ ultrametric tree was adapted from Mckenna et al. [6]. Each term contains genes that produce the same protein. Abbreviations of the tested species can be found in Table 1

Feeding experiment and transcriptome profiling

The responses of both life history traits and gene expression to different diets were experimentally studied. Second-instar C. montrouzieri larvae were raised on 13 diets, including one natural prey diet and 12 factitious prey or artificial diets (Table 3). The natural prey Planococcus citri (Risso, 1813) (MEALYBUG) has been used to maintain the tested laboratory population for more than 10 years. Thus, the use of this prey for C. montrouzieri as a control allows comparison of the responses to different diets.

Table 3 Diet design for Cryptolaemus montrouzieri larvae

In the comparison of life history traits, we found large differences in the performance of the larvae among these 13 diets (Fig. 3 and details in Additional file 1: Table S3). The natural prey diet (MEALYBUG) was clearly favorable, with the highest adult weight, second shortest development time and lowest mortality rate (Fig. 3a). The two factitious prey diets PEAAPHID and FLOURMOTH were second only to the natural prey diet in terms of adult weight. Individuals in the remaining ten diet treatments performed much worse, with > 70% mortality or failure to develop to the adult stage on six of those diets. The 12 unnatural diets overall led to significantly longer larval survival than no food (Fig. 3b), especially the POLLEN diet, which sustained larvae for up to 50 days.

Fig. 3
figure 3

Comparison of life history traits of C. montrouzieri fed different diets. a Effect of different diets on the development, adult weight and mortality of C. montrouzieri larvae. Error bars show the standard deviation. b Survival time of larvae with different diets that did not allow development to the adult stage

Gene expression was then profiled in 4th-instar larvae (< 24 h after molting) fed different diets. An overview of relative gene expression levels in the 13 treatments is presented by a heatmap of r2 values in Fig. 4. All of the treatments had r2 values between the two replicates exceeding 0.88, indicating repeatability within treatments. The top three favorable treatments in life history trait comparisons (MEALYBUG, PEAAPHID and FLOURMOTH) shared also high r2 values among each other (0.82–0.93, mean = 0.88). Gene expression patterns in these three treatments were usually more different from the inferior diet treatments (r2: 0.56–0.88, mean = 0.73). When comparing the 12 unnatural (i.e. factitious prey or artificial diet) treatments with the MEALYBUG treatment, DEGs were enriched in 32 KEGG pathways (Q < 0.05), among which 29 were related to nutrient or toxin metabolism (Fig. 5). Similarly, the DEGs were enriched in GO terms mainly related to nutrient metabolism processes and catalytic/oxidoreductase activities (Q < 0.05, Additional file 1: Fig. S2).

Fig. 4
figure 4

Relationship (r2) of gene expression between the studied transcriptomes with different diet treatments. Abbreviations of diet treatments can be found in Table 3

Fig. 5
figure 5

Heatmap of adjusted P values (Q) in KEGG pathway enrichment analysis for the transcriptome comparisons of alternative diets versus the natural prey of C. montrouzieri larvae. Enrichment with Q < 0.05 is marked with an asterisk. Twenty-nine out of 32 enriched pathways were related to metabolism

As we found a significant expansion in a gene family involved in the immune response in C. montrouzieri (Table 2), the pattern of expression of immune effector genes including those encoding antibacterial peptides (AMPs) and lysozymes was specifically analyzed. We found a general down regulation of the immune effector genes (log2-fold change mean ± SE: − 1.73 ± 0.15) as compared to the natural prey MEALYBUG treatment. Among them, 6 of 28 were dramatically downregulated when larvae shifted their diet from mealybugs to unnatural diets, with most of the log2-fold change values lower than − 3 (i.e. nine times lower than those in the natural prey control, Fig. 6). These genes include two attacin genes, two defensin genes, one coleoptericin gene and one cwh gene. In contrast, only slight regulation of expression (log2-fold change mean ± SE: − 0.69 ± 0.08) was detected in genes related to immune recognition including those encoding c-type lectin, peptidoglycan recognition protein (PGRP) and gram-negative binding protein (GNBP) (Fig. 6).

Fig. 6
figure 6

Transcriptome profiling of C. montrouzieri larvae fed different diets. Heatmap shows the log2(fold change) for the transcriptome comparisons of alternative diets versus the natural prey of C. montrouzieri. Hierarchical clustering of the heatmap was carried out using pheatmap package in R. Only genes related to immune recognition and response are shown. Abbreviations of diet treatments can be found in Table 3. PGRP: peptidoglycan recognition protein. GNBP: gram-negative binding protein


Gene expansion/contraction related to feeding habits

The order Coleoptera is the most speciose group of animals with highly diverse feeding habits. Most of the species in the suborder Adephaga are predaceous while Polyphaga (e.g. weevils, longhorn beetles and leaf beetles) are predominantly phytophagous species. The high diversity of phytophagous beetles can be explained by their complex interactions with flowering plants [8, 31, 32]. However, Polyphaga also includes the ladybird beetles (Coccinellidae), most of which are predaceous [33]. Evolutionary studies have suggested that the ancestral ladybirds have switched from mycophagy to a predatory life style [33,34,35]. This is associated with the diversification of ladybirds into more than 6000 reported species [34].

In this study, we explored the evolutionary patterns of gene families involved in the functions of chemosensing, digestion, detoxification and immunity in the predatory ladybird C. montrouzieri as compared with other beetles with different feeding habits. We found that the C. montrouzieri genome has undergone significant expansion or contraction of several gene families encoding chemosensors, digestive and detoxification enzymes. It seems that these gene families are usually involved in diet adaptation of not only phytophagous but also predatory beetles. The evolution of these gene families of C. montrouzieri might be associated with adaptation to mealybug feeding. However, we also need to be aware of the other potential factors (e.g. the distinct phylogenetic relationship of the tested species, the different qualities of genome assemblies and annotations) that affect these patterns of genome evolution.

Deficiencies of novel diets and diet-induced gene regulation

The availability of a cost-effective factitious prey or artificial diet is key to the successful mass production of arthropod natural enemies for use in biological control [36]. The availability of molecular markers that could be used as early indicators of insect responses to general and specific nutritional levels may greatly advance the practice of insect mass rearing [37,38,39,40]. However, very few biological control agent species have completely sequenced genomes. In this study, with the help of our new completely sequenced C. montrouzieri genome, diet-induced gene regulation was explored, which might further benefit diet development. We found that the performance of C. montrouzieri larvae fed unnatural diets decreased to different degrees. Both life history trait performance and the pattern of gene expression congruently revealed that the feeding treatments with another Sternorrhyncha species (PEAAPHID) or lepidopteran eggs (FLOURMOTH) were relatively close to those with natural prey. In comparison, individuals in the rest of the treatments performed much worse and had more divergent patterns of gene expression. DEGs in the unnatural diet treatments compared to the natural prey treatment were mainly enriched in nutrient metabolism. This finding explores a relationship between diet and regulation of gene expression, and suggests that the decrease in performance might be caused by a dietary nutrient imbalance. Several studies have tried to explore the relationship between gene expression patterns and diet limitation of predatory biological control agents, and have tried to further optimize the diet formulations [37,38,39,40,41]. However, the mechanism of how diet components affect gene expression of these predators is still not clear. Similarly, our current data are insufficient to pin down the exact deficiencies of the tested diets. More detailed studies are needed to explore the relationship between diet components and gene expression.

Potential roles of immune enhancement in prey adaptation

Previous studies have demonstrated that many genes encoding AMPs and lysozymes of ladybirds are dramatically upregulated under the challenge of bacterial infection [12, 42]. In our study, the expression of some genes encoding AMPs and cell wall hydrolases was dramatically downregulated when using unnatural diets. Change of nutritional condition is known to cause significant changes in the physiology of ladybirds [43] and other predatory insects [44]. Thus, it is possible that using an unnatural diet has a negative effect on the physiology of ladybirds, including immune defense. Several studies have demonstrated the impact of diet on insect immune response [45,46,47]. However, their findings show that diet changes can positively impact on certain immune traits but negatively affect others. In addition, little evidence in insects supports a relationship between diet changes and downregulation of immune effector genes [48]. Alternatively, it seems that using natural prey, i.e. mealybugs, is one of the factors inducing immune responses in C. montrouzieri. Furthermore, the gene family encoding attacin has expanded in comparison to other Coleoptera genomes. These two pieces of evidence together suggest a potential role of immune enhancement in prey adaptation of C. montrouzieri.

The main prey of predatory ladybirds include aphids, coccids (e.g., mealybugs and scale insects), psyllids and whiteflies, all of which are in the suborder Sternorrhyncha (Hemiptera). These sap-feeding insects have a specific bacteriome that harbors a large number of symbiont bacteria [49]. It is possible that the symbiont bacteria of mealybugs cause an immune response in ladybird predators. A growing body of evidence shows that bacterial symbionts can protect their hosts from parasites and predators [50,51,52]. For example, the symbiont of Paederus beetles synthesizes a chemical toxin that beetles can use as a defense against predators [53]. Also, the symbionts of the ladybird H. axyridis produce pyrazines that have a function in defense behavior [54]. It would be interesting to investigate whether and which bacteria in prey cause immune responses in their predators.


The high-quality whole genome assembly of the predatory ladybird C. montrouzieri provides insights into its diet adaptation, including the expansion or contraction of the gene families encoding chemosensors and digestive and detoxifying enzymes, and the enrichment of differentially expressed genes in pathways of nutrient metabolism when using unnatural diets. We highlight the potential role of immune enhancement and evolution of genes encoding immune effectors in diet adaptation of this species. This genomic study of a biological control agent is valuable for improving our basic understanding of its feeding habits, and may assist in improving its utilization in biological control.


DNA extraction, genome sequencing and assembly

DNA was extracted from the whole body of ten female adults of C. montrouzieri. These individuals were derived from a population that has been reared under laboratory conditions (27 ± 1 °C, 80 ± 1% relative humidity (RH) and a 14:10 (L:D) h photoperiod) at Sun Yat-sen University since 2005 [55]. Genomic DNA was extracted using the CTAB method [56]. The quality and concentration of the extracted genomic DNA were checked using 1% agarose gel electrophoresis and a Qubit fluorimeter (Invitrogen, Carlsbad, CA, USA). High-quality DNA was used for subsequent Nanopore and Illumina sequencing.

Approximately 15 μg of genomic DNA was used to generate Oxford Nanopore long reads, and the sequencing reaction was performed in a PromethION DNA sequencer (Oxford Nanopore, Oxford, UK). The raw data were then filtered to remove short sequence reads (< 5 kb) and reads with low-quality bases (Q30 < 90%) using Nanofilt v2.3.0 [57]. For assembly of Nanopore sequencing data, Canu v1.5 [58] was implemented to generate more accurate self-corrected reads with a corrected error rate of 0.05. Assembly was then performed by Wtdbg ( with default settings. Racon v1.32 [59] was implemented to correct the assembly with Nanopore reads through two rounds with default settings. For further error correction, genomic DNA was also sequenced on the Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA). The Illumina sequenced data were filtered to remove reads with low-quality bases and adapters using Trimmomatic v0.36 [60] with default settings. Pilon v1.21 [61] was implemented to correct the Nanopore assembly with Illumina reads through three rounds with default settings.

Gene prediction and functional annotation

First, the repetitive elements of the assembled genome were identified and masked. Repetitive elements of the assembled genome were classified into families with five rounds of RepeatModeler v2.0.1 analysis with default settings, followed by genome masking with RepeatMasker v4.1.0 ( with default settings. Second, genes were automatically predicted based on our RNA-Seq data of C. montrouzieri in different life stages and diet treatments (SRA accession: SRR2971112, SRR2971116, SRR6981477, SRR8325176, SRR8325159). RNA-Seq reads were mapped to the assembled genome sequence using HISAT2 v2.1.0 [62]. Augustus 3.0.3 [63] was trained with the mapped data according to the Braker2 pipeline [24], and further used to predict genes in the genome sequences ab initio. BUSCO v4.1.2 [23] with the Insecta set of the OrthoDB v10.1 database was then used to assess the quality of the predicted gene set.

The protein sequences of the automatically predicted genes were subjected to similarity searches against the NCBI NR Hexapoda subset and UniProtKB Swiss-Prot databases using BLASTp with a cutoff E-value of 10− 5. Only the longest protein isoform for each gene was used as a query. Protein domains within genes were searched against the Pfam v32 database using InterProScan v5 [64] with a cutoff E-value of 10− 5. Sequences were also mapped to the GO database [65, 66] and KEGG reference pathways [67] using eggNOG-mapper [68] with a cutoff E-value of 10− 5.

In addition, specific genes with immune functions of beetles were identified based on searches against the databases UniProtKB/Swiss-Prot, Pfam v32 or KEGG pathway. These genes are classified into three major roles in immunity: recognition, signaling cascade and response [12, 42, 69]. Recognition genes included c-type lectin (containing the Pfam protein domain: ID:PF00059), peptidoglycan recognition protein (PGRP, identified by Swiss-Prot annotation) and gram-negative binding protein (GNBP, identified by Swiss-Prot annotation). Signaling cascade genes included those in the Toll and IMD pathway (KEGG: map 04624) and JAK/STAT pathway (KEGG: map 04630). Immune response genes included those encoding antimicrobial peptides, e.g., attacin (containing the Pfam domain: ID: PF03769), defensin (PF01097), coleoptericin (PF06286), thaumatin (PF00314) and apolipophorin (PF07464), and lysozymes, e.g., c-type and i-type lysozyme (identified by Swiss-Prot annotation). In addition, a putative antimicrobial gene, cell wall hydrolase (containing the Pfam domain: ID: PF07486), in ladybirds was also included in the analyses [70].

Orthology search and gene family evolution

OrthoFinder v2 [71] was used to identify orthologous genes by retrieving the protein sequences of the C. montrouzieri and the other published Coleoptera genomes with default settings. A total of nine gene sets of Coleoptera genomes predicted from RNA-Seq data were used (Table 1), and their longest protein isoforms were extracted as input of OrthoFinder. Protein domains within genes were searched against the Pfam v32 database using InterProScan v5 with a cutoff E-value of 10− 5. Information of protein domain was subsequently assigned to the orthogroups using KinFin [72]. Furthermore, these orthogroups were used as input for CAFE v4.1 [30] to assess gene family contraction and expansion dynamics using the birth/death parameter (λ). The species tree used in CAFE was adapted from a recently published Coleoptera phylogeny [6]. In each branch, orthologous groups with p-values < 0.05 were considered significant expansions or contractions.

Feeding experiment: diet-specific life history traits and transcriptome

Diet materials for C. montrouzieri were selected based on common protein sources used for insect diets. The selected protein sources covered whole bodies of seven invertebrate species, the eggs of three invertebrate species, two types of vertebrate products and one type of plant tissue (Table 3). Most solid materials were first dried in an oven (60 ± 1 °C) for 24 h and then ground to a powder using a kitchen blender. For materials that could not be directly fed to ladybirds, media were made with 2.5 g of protein source, 1.5 g of sucrose, 0.23 g of agar, and 19 mL of distilled water; sucrose was included as a feeding stimulant and agar as a thickener. A preliminary blank control treatment in which C. montrouzieri was offered the medium without a protein source from the 2nd instar onwards showed that none of the larvae developed to the 4th instar. We also set up a starvation treatment with no food or water provided by the 2nd instar onwards.

The developmental traits of C. montrouzieri were investigated from the 2nd instar onwards because we observed that some first instars died from non-nutritional factors (e.g., they were stuck in the medium). Before the treatments, all larvae were fed mealybugs in communal cultures. Thereafter, 52 to 137 2nd-instar larvae of C. montrouzieri (< 24 h after molting) derived from the population at Sun Yat-sen University [55] were placed individually in plastic Petri dishes (diameter: 5 cm, height: 2 cm) for the different diet treatments. All diets were offered ad libitum and replenished daily. The survival and development of C. montrouzieri larvae and pupae were monitored daily. Newly emerged adults were weighed. All feeding experiments were performed in a climatic chamber at 27 ± 1 °C with an 80 ± 1% relative humidity (RH) and a 14:10 (L:D) h photoperiod. A Kolmogorov-Smirnov test indicated that survival time of larvae and adult weight were normally distributed and therefore could be analyzed using a one-way analysis of variance (ANOVA). As a Levene test indicated homoscedasticity, the means were separated using Tukey tests. In all tests, p values below 0.05 were considered significant. All data were analyzed using SPSS 17.0 (SPSS Inc.).

Two 4th-instar larvae (< 24 h after molting) of C. montrouzieri from each diet treatment in the above life history experiment were randomly collected for transcriptome analysis. After ~ 12 h of starvation, the total RNA of each individual was extracted using TRIzol reagent (CWBIO, Beijing, China). RNA quality and quantity were determined using a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, US). Only RNA samples with a 260/280 ratio from 1.8 to 2.0, a 260/230 ratio from 2.0 to 2.5 and an RNA integrity number (RIN) greater than 8.0 were used for sequencing. Sequencing was performed on the Illumina HiSeq 2500 platform, generating 2 × 125 bp reads. Adaptors and low-quality sequences were removed using the default settings for Trimmomatic v0.36.

Transcript assembly and abundance estimation were performed using the TopHat2 + Cufflinks method [73]. The coefficient of determination (r2) from Pearson’s correlation analysis was used to analyze the relationship of each sample pair based on fragments per kilobase of transcript per million mapped reads (FPKM) values. The regulation of gene expression was tested using Cuffdiff in Cufflinks, with a false discovery rate (FDR) < 0.05 for defining differentially expressed genes (DEGs). We used natural prey as a control to test for transcriptional regulation when ladybirds were fed unnatural diets. Thus, analyses of DEGs were performed only in the comparisons of the 12 factitious prey or artificial diet treatments against the MEALYBUG treatment. We investigated which GO terms and KEGG pathways the DEGs were involved in and evaluated statistical significance of GO and KEGG enrichment by hypergeometric distribution testing.

Availability of data and materials

The raw and assembled sequenced data of C. montrouzieri were deposited in the NCBI BioProject: PRJNA626074. The RNA-Seq raw data of different diet treatments of C. montrouzieri were deposited in the Sequence Read Archive (SRA) repository of the NCBI under accession Nos. SRR8325159 - SRR8325188 of BioProject PRJNA509782.


  1. Eyun SI, Soh HY, Posavi M, Munro JB, Hughes DST, Murali SC, et al. Evolutionary history of chemosensory-related gene families across the Arthropoda. Mol Biol Evol. 2017;34:1838–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Pearce SL, Clarke DF, East PD, Elfekih S, Gordon KHJ, Jermiin LS, et al. Genomic innovations, transcriptional plasticity and gene loss underlying the evolution and divergence of two highly polyphagous and invasive Helicoverpa pest species. BMC Biol. 2017;15:63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Cheng T, Wu J, Wu Y, Chilukuri RV, Huang L, Yamamoto K, et al. Genomic adaptation to polyphagy and insecticides in a major East Asian noctuid pest. Nat Ecol Evol. 2017;1:1747–56.

    Article  PubMed  Google Scholar 

  4. Schoville SD, Chen YH, Andersson MN, Benoit JB, Bhandari A, Bowsher JH, et al. A model species for agricultural pest genomics: the genome of the Colorado potato beetle, Leptinotarsa decemlineata (Coleoptera: Chrysomelidae). Sci Rep. 2018;8:1931.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  5. McKenna DD, Scully ED, Pauchet Y, Hoover K, Kirsch R, Geib SM, et al. Genome of the Asian longhorned beetle (Anoplophora glabripennis), a globally significant invasive species, reveals key functional and evolutionary innovations at the beetle-plant interface. Genome Biol. 2016;17:227.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. McKenna DD, Shin S, Ahrens D, Balke M, Beza-Beza C, Clarke DJ, et al. The evolution and genomic basis of beetle diversity. P Natl Acad Sci U S A. 2019;116:24729–37.

    Article  CAS  Google Scholar 

  7. Hazzouri KM, Sudalaimuthuasari N, Kundu B, Nelson D, Al-Deeb MA, Le Mansour A, et al. The genome of pest Rhynchophorus ferrugineus reveals gene families important at the plant-beetle interface. Commun Biol. 2020;3:323.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Seppey M, Ioannidis P, Emerson BC, Pitteloud C, Robinson-Rechavi M, Roux J, et al. Genomic signatures accompanying the dietary shift to phytophagy in polyphagan beetles. Genome Biol. 2019;20:98.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Vogel H, Muller A, Heckel DG, Gutzeit H, Vilcinskas A. Nutritional immunology: diversification and diet-dependent expression of antimicrobial peptides in the black soldier fly Hermetia illucens. Dev Comp Immunol. 2018;78:141–8.

    Article  CAS  PubMed  Google Scholar 

  10. Vilcinskas A. Evolutionary plasticity of insect immunity. J Insect Physiol. 2013;59:123–9.

    Article  CAS  PubMed  Google Scholar 

  11. Login FH, Balmand S, Vallier A, Vincent-Monegat C, Vigneron A, Weiss-Gayet M, et al. Antimicrobial peptides keep insect endosymbionts under control. Science. 2011;334:362–5.

    Article  CAS  PubMed  Google Scholar 

  12. Vogel H, Schmidtberg H, Vilcinskas A. Comparative transcriptomics in three ladybird species supports a role for immunity in invasion biology. Dev Comp Immunol. 2017;67:452–6.

    Article  CAS  PubMed  Google Scholar 

  13. Pang H, Tang XF, Booth RG, Vandenberg N, Forrester J, Mchugh J, et al. Revision of the Australian Coccinellidae (Coleoptera). Genus Novius Mulsant of tribe Noviini. Ann Zool. 2020;70:1–24.

    Article  Google Scholar 

  14. Caltagirone LE, Doutt RL. The history of the vedalia beetle importation to California and its impact on the development of biological control. Annu Rev Entomol. 1989;34:1–16.

    Article  Google Scholar 

  15. De Clercq P, Mason PG, Babendreier D. Benefits and risks of exotic biological control agents. BioControl. 2011;56:681–98.

    Article  Google Scholar 

  16. Kairo MKT, Paraiso O, Gautam RD, Peterkin DD. Cryptolaemus montrouzieri (Mulsant) (Coccinellidae: Scymninae): a review of biology, ecology, and use in biological control with particular reference to potential impact on non-target organisms. CAB Reviews. 2013;8:1–20.

    Article  Google Scholar 

  17. Maes S, Antoons T, Gregoire JC, De Clercq P. a semi-artificial rearing system for the specialist predatory ladybird Cryptolaemus montrouzieri. BioControl. 2014;59:557–64.

    Article  Google Scholar 

  18. Xie J, Wu H, Pang H, De Clercq P. An artificial diet containing plant pollen for the mealybug predator Cryptolaemus montrouzieri. Pest Manag Sci 2016; 73:541–545.

  19. Maes S, Grégoire J-C, De Clercq P. Prey range of the predatory ladybird Cryptolaemus montrouzieri. BioControl. 2014;59:729–38.

    Article  Google Scholar 

  20. Surwase S, Shetgar S, Khandare R, Magar S, Nalwandikar P. Biology of Cryptolaemus montrouzieri Mulsant on mealy bugs and aphids. J Entomol Res. 2016;40:95–100.

    Article  Google Scholar 

  21. Li H-S, Pan C, De Clercq P, Ślipiński A, Pang H. Variation in life history traits and transcriptome associated with adaptation to diet shifts in the ladybird Cryptolaemus montrouzieri. BMC Genomics. 2016;17:281.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  22. Zhang L, Li S, Luo J, Du P, Wu L, Li Y, et al. Chromosome-level genome assembly of the predator Propylea japonica to understand its tolerance to insecticides and high temperatures. Mol Ecol Resour. 2019.

  23. Waterhouse RM, Seppey M, Simao FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2018;35:543–8.

    Article  CAS  PubMed  Google Scholar 

  24. Hoff KJ, Lange S, Lomsadze A, Borodovsky M, Stanke M. BRAKER1: unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS. Bioinformatics. 2016;32:767–9.

    Article  CAS  PubMed  Google Scholar 

  25. Thomas GWC, Dohmen E, Hughes DST, Murali SC, Poelchau M, Glastad K, et al. Gene content evolution in the arthropods. Genome Biol. 2020;21:15.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Fallon TR, Lower SE, Chang CH, Bessho-Uehara M, Martin GJ, Bewick AJ, et al. Firefly genomes illuminate parallel origins of bioluminescence in beetles. eLife. 2018;7:e36495.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Cunningham CB, Ji L, Wiberg RA, Shelton J, McKinney EC, Parker DJ, et al. The genome and methylome of a beetle with complex social behavior, Nicrophorus vespilloides (Coleoptera: Silphidae). Genome Biol Evol. 2015;7:3383–96.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Richards S, Gibbs RA, Weinstock GM, Brown SJ, Denell R, Beeman RW, et al. The genome of the model beetle and pest Tribolium castaneum. Nature. 2008;452:949–55.

    Article  CAS  PubMed  Google Scholar 

  29. Keeling CI, Yuen MM, Liao NY, Docking TR, Chan SK, Taylor GA, et al. Draft genome of the mountain pine beetle, Dendroctonus ponderosae Hopkins, a major forest pest. Genome Biol. 2013;14:R27.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Han MV, Thomas GW, Lugo-Martinez J, Hahn MW. Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. Mol Biol Evol. 2013;30:1987–97.

    Article  CAS  PubMed  Google Scholar 

  31. Farrell BD. "Inordinate fondness" explained: why are there so many beetles? Science. 1998;281:555–9.

    Article  CAS  PubMed  Google Scholar 

  32. Zhang SQ, Che LH, Li Y, Dan L, Pang H, Slipinski A, et al. Evolutionary history of Coleoptera revealed by extensive sampling of genes and species. Nat Commun. 2018;9:205.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Giorgi JA, Vandenberg NJ, McHugh JV, Forrester JA, Ślipiński SA, Miller KB, et al. The evolution of food preferences in Coccinellidae. Biol Control. 2009;51:215–31.

    Article  Google Scholar 

  34. Seago AE, Giorgi JA, Li J, Slipinski A. Phylogeny, classification and evolution of ladybird beetles (Coleoptera: Coccinellidae) based on simultaneous analysis of molecular and morphological data. Mol Phylogenet Evol. 2011;60:137–51.

    Article  PubMed  Google Scholar 

  35. Magro A, Lecompte E, Magne F, Hemptinne JL, Crouau-Roy B. Phylogeny of ladybirds (Coleoptera: Coccinellidae): are the subfamilies monophyletic? Mol Phylogenet Evol. 2010;54:833–48.

    Article  CAS  PubMed  Google Scholar 

  36. Grenier S. In vitro rearing of entomophagous insects- past and future trends: a minireview. B Insectol. 2009;62:1–6.

    Google Scholar 

  37. Coudron TA, Yocum GD, Brandt SL. Nutrigenomics: a case study in the measurement of insect response to nutritional quality. Entomol Exp Appl. 2006;121:1–14.

    Article  CAS  Google Scholar 

  38. Yocum GD, Coudron TA, Brandt SL. Differential gene expression in Perillus bioculatus nymphs fed a suboptimal artificial diet. J Insect Physiol. 2006;52:586–92.

    Article  CAS  PubMed  Google Scholar 

  39. Zou DY, Coudron TA, Liu CX, Zhang LS, Wang MQ, Chen HY. Nutrigenomics in Arma chinensis: transcriptome analysis of Arma chinensis fed on artificial diet and Chinese oak silk moth Antheraea pernyi pupae. PLoS One. 2013;8:e60881.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Zou DY, Coudron TA, Zhang LS, Gu XS, Xu WH, Liu XL, et al. Performance of Arma chinensis reared on an artificial diet formulated using transcriptomic methods. Bull Entomol Res. 2018;109:24–33.

    Article  PubMed  CAS  Google Scholar 

  41. Coudron TA, Brandt SL, Hunter WB. Molecular profiling of proteolytic and lectin transcripts in Homalodisca vitripennis (Hemiptera: Auchenorrhyncha: Cicadellidae) feeding on sunflower and cowpea. Arch Insect Biochem. 2007;66:76–88.

    Article  CAS  Google Scholar 

  42. Vilcinskas A, Mukherjee K, Vogel H. Expansion of the antimicrobial peptide repertoire in the invasive ladybird Harmonia axyridis. P Roy Soc B-Biol Sci. 2013;280:20122113.

    Google Scholar 

  43. Sun Y-X, Hao Y-N, Riddick EW, Liu T-X. Factitious prey and artificial diets for predatory lady beetles: current situation, obstacles, and approaches for improvement: a review. Biocontrol Sci Tech. 2017;27:601–19.

    Article  Google Scholar 

  44. Riddick EW. Benefits and limitations of factitious prey and artificial diets on life parameters of predatory beetles, bugs, and lacewings: a mini-review. BioControl. 2008;54:325–39.

    Article  Google Scholar 

  45. Alaux C, Ducloz F, Crauser D, Le Conte Y. Diet effects on honeybee immunocompetence. Biol Lett. 2010;6:562–5.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Cotter SC, Simpson SJ, Raubenheimer D, Wilson K. Macronutrient balance mediates trade-offs between immune function and life history traits. Funct Ecol. 2010;25:186–98.

    Article  Google Scholar 

  47. Srygley RB, Lorch PD. Weakness in the band: nutrient-mediated trade-offs between migration and immunity of Mormon crickets, Anabrus simplex. Anim Behav. 2011;81:395–400.

    Article  Google Scholar 

  48. Jacobs CG, Steiger S, Heckel DG, Wielsch N, Vilcinskas A, Vogel H. Sex, offspring and carcass determine antimicrobial peptide expression in the burying beetle. Sci Rep. 2016;6:25409.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Baumann P. Diversity of prokaryote–insect associations within the Sternorrhyncha (psyllids, whiteflies, aphids, mealybugs). In: Bourtzis K, Miller T, editors. Insect Symbiosis, vol. 2: CRC Press; 2006. p. 23–46.

  50. Brownlie JC, Johnson KN. Symbiont-mediated protection in insect hosts. Trends Microbiol. 2009;17:348–54.

    Article  CAS  PubMed  Google Scholar 

  51. McLean AHC. Cascading effects of defensive endosymbionts. Curr Opin Insect Sci. 2019;32:42–6.

    Article  PubMed  Google Scholar 

  52. Haine ER. Symbiont-mediated protection. Proc Biol Sci. 2008;275:353–61.

    PubMed  Google Scholar 

  53. Piel J, Hui D, Wen G, Butzke D, Platzer M, Fusetani N, et al. Antitumor polyketide biosynthesis by an uncultivated bacterial symbiont of the marine sponge Theonella swinhoei. P Natl Acad Sci U S A. 2004;101:16222–7.

    Article  CAS  Google Scholar 

  54. Schmidtberg H, Shukla SP, Halitschke R, Vogel H, Vilcinskas A. Symbiont-mediated chemical defense in the invasive ladybird Harmonia axyridis. Ecol Evol. 2019;9:1715–29.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Li H-S, Heckel G, Huang Y-H, Fan W-J, Slipinski A, Pang H. Genomic changes in the biological control agent Cryptolaemus montrouzieri associated with introduction. Evol Appl. 2019;12:989–1000.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Milligan BG. Total DNA isolation. In: Hoelzel AR, editor. Molecular genetic analysis of populations. Oxford: Oxford University Press; 1988. p. 29–64.

    Google Scholar 

  57. De Coster W, D'Hert S, Schultz DT, Cruts M, Van Broeckhoven C. NanoPack: visualizing and processing long-read sequencing data. Bioinformatics. 2018;34:2666–9.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  58. Koren S, Walenz BP, Berlin K, Miller JR, Bergman NH, Phillippy AM. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 2017;27:722–36.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Vaser R, Sovic I, Nagarajan N, Sikic M. Fast and accurate de novo genome assembly from long uncorrected reads. Genome Res. 2017;27:737–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Walker BJ, Abeel T, Shea T, Priest M, Abouelliel A, Sakthikumar S, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS One. 2014;9.

  62. Kim D, Landmead B, Salzberg SLHISAT. A fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Stanke M, Waack S. Gene prediction with a hidden Markov model and a new intron submodel. Bioinformatics. 2003;19:Ii215–25.

    Article  PubMed  Google Scholar 

  64. Jones P, Binns D, Chang HY, Fraser M, Li WZ, McAnulla C, et al. InterProScan 5: genome-scale protein function classification. Bioinformatics. 2014;30:1236–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. The Gene Ontology C. The Gene Ontology resource: 20 years and still GOing strong. Nucleic Acids Res. 2019;47:D330–8.

    Article  CAS  Google Scholar 

  66. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Gene Ontol Consortium Nat Genet. 2000;25:25–9.

    CAS  Google Scholar 

  67. Kanehisa M, Goto SKEGG. Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Jaime HC, Kristoffer F, Pedro CL, Damian S, Juhl JL, vM C, et al. Fast genome-wide functional annotation through Orthology assignment by eggNOG-mapper. Mol Biol Evol. 2017;34:2115–22.

    Article  CAS  Google Scholar 

  69. Altincicek B, Knorr E, Vilcinskas A. Beetle immunity: identification of immune-inducible genes from the model insect Tribolium castaneum. Dev Comp Immunol. 2008;32:585–95.

    Article  CAS  PubMed  Google Scholar 

  70. Li H-S, Tang X-F, Huang Y-H, Xu Z-Y, Chen M-L, Du X-Y, et al. Horizontally acquired antibacterial genes associated with adaptive radiation of ladybird beetles. BMC Biol. 2021;19:7.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Emms DM, Kelly S. OrthoFinder: phylogenetic orthology inference for comparative genomics. Genome Biol. 2019;20:238.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Laetsch DR, Blaxter ML. KinFin: software for taxon-aware analysis of clustered protein sequences. G3-Genes Genom Genet. 2017;7:3349–57.

    CAS  Google Scholar 

  73. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


We would like to thank Li-Jun Ma (Sun Yat-sen University) for assistance with the experiments.


This work was supported by the National Natural Science Foundation of China (Grant No. 31702012, 31970439), Science and Technology Planning Project of Guangdong Province, China (Grant No. 2017B020202006), National Key R&D Program of China (Grant No. 2017YFD0201000), Basal Research Fund of Sun Yat-sen University (Grant No. 18lgpy50), Science and Technology Planning Project of Guangzhou, China (Grant No. 201904020041) and Open Project of the State Key Laboratory of Biocontrol (Grant No. 2018–04).

Author information

Authors and Affiliations



LHS and PH designed the study. LHS, HYH, CML, RZ and QBY performed the laboratory work. LHS and HYH analyzed the data. LHS, DCP, HG and PH wrote the manuscript. All authors gave final approval for publication.

Corresponding author

Correspondence to Hong Pang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary Figure S1.

Results of CAFE analysis inferring the size change of gene families in the genome of nine tested Coleoptera species. This summary tree shows the number of expanded (red) and contracted (green) families. The species’ ultrametric tree was adapted from Mckenna et al. [6]. Abbreviations of the tested species are defined in Table 1. Table S3. Comparison of life history traits of Cryptolaemus montrouzieri subjected to different diet treatments and starvation. Means (± SE) within a column followed by the same letter are not significantly different (P > 0.05). Figure S2. Heatmap of adjusted P values (Q) in Gene Ontology (GO) enrichment analysis for the transcriptome comparisons of alternative diets versus the natural prey of C. montrouzieri larvae. C: Cellular component; F: Molecular function; P: Biological process. The number of background genes of each GO term is shown in brackets. Enrichment with Q < 0.05 is marked with an asterisk.

Additional file 2: Table S1.

Functional annotation of the gene set predicted from the genome of Cryptolaemus montrouzieri. Table S2. Information on the orthologous groups identified by OrthoFinder, number of genes in each genome/orthologous group and the p-value in the branch of Cryptolaemus montrouzieri estimated by CAFE. Table S4. Details of transcriptome profiling of Cryptolaemus montrouzieri under different diet treatments. Expression levels of genes of the 12 factitious prey or artificial diet treatments were compared to those of the mealybug treatment. Values of fragments per kilobase of transcript per million mapped reads (FPKM) of each gene/treatment and log2(fold change) and adjusted p of each gene/comparison are shown.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, HS., Huang, YH., Chen, ML. et al. Genomic insight into diet adaptation in the biological control agent Cryptolaemus montrouzieri. BMC Genomics 22, 135 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: