Evolutionary genomics of camouflage innovation in the orchid mantis

Sample collection

Captive breeding individuals of H. coronatus (Mantodea, Hymenopodidae) hatched from the same ootheca that was collected from the Xishuangbanna rainforest, Yunnan Province, China in 2018. Individuals of D. lobata (Mantodea, Deroplatyidae) were collected from a captive breeding center in Beijing, China in 2018. All individuals were housed in semitransparent cages (7 cm × 7 cm × 9 cm) in an incubator with a temperature range of 26–27 °C, 80% relative humidity, and a 16 h/8 h light/dark photoperiod. Larvae from the first to third instars were fed fruit flies and switched to crickets from the fourth instar onward. To ensure sufficient output of nucleic acids for Nanopore sequencing and Illumina paired-end sequencing, adult female H. coronatus and D. lobata were used for nucleic acid preps. DNA was extracted using the DNeasy Blood and Tissue Kit following the manufacturer’s protocol (Qiagen, Valencia, CA, USA). The quantity and quality of DNA were checked by a Qubit 2.0 (Thermo Fisher Scientific) and agarose gel electrophoresis, respectively. RNA was extracted using TRIzol reagent (Thermo Fisher Scientific).

All sample collections and animal experimental procedures were approved by the Institutional Animal Care and Use Committee of the Institute of Zoology, Chinese Academy of Sciences. All efforts were made to minimize the suffering of the animals.

Genome sequencing, Hi-C sequencing, and genome assembly

A combination of Nanopore sequencing, Illumina sequencing and Hi-C sequencing was used to generate mantis genome assemblies. Paired-end libraries with insert sizes of 400 bp were constructed and sequenced using the Illumina HiSeq X-Ten platform following the standard manufacturer’s protocol (San Diego, USA), yielding 155.6 Gb of clean reads for H. coronatus and 286.3 Gb for D. lobata, followed by quality filtering (Supplementary Table 9). The high-quality reads were used for genome size estimation by the K-mer method. 1D libraries were constructed according to the manufacturer’s instructions (SQK-LSK109, Oxford Nanopore Technologies, UK) and sequenced on the Nanopore PromethION platform (Oxford Nanopore Technologies, UK). In total, 297.60 Gb of subreads for H. coronatus and 258.05 Gb for D. lobata were generated with a high sequencing depth (approximately, 66–103×) (Supplementary Table 9). To further improve the continuity of the assembled genomes and anchor the assemblies into chromosomes, Hi-C sequencing was performed to order and orient the contigs, as well as to correct misjoined sections and merge overlaps. To construct Hi-C libraries, tissue samples were fixed immediately in formaldehyde solution once isolated. DNA was digested with the restriction enzyme DpnII overnight. Hi-C libraries with a mean size of 350 bp were constructed using NEBNext Ultra enzymes and Illumina-compatible adaptors and sequenced using the Illumina HiSeq 2500 platform. All libraries were quantified by a Qubit 2.0 (Thermo Fisher Scientific), and insert size was checked using an Agilent 2100.

Contigs were assembled by NextDenovo software (v2.0-beta.1, github.com/Nextomics/NextDenovo), and the assembled contig-level genomes were polished by NextPolish (v1.0.5)⁴⁴ and Pilon (v1.22)⁴⁵. Then, the contigs were anchored into chromosomes by Hi-C sequencing reads through the Juicer (v1.5)⁴⁶ and 3D-DNA (v180922)⁴⁷ software workflows. To further improve the chromosome-scale assembly, it was subjected to manual review and refinement using Juicebox Assembly Tools (github.com/theaidenlab/juicebox). Finally, genome quality was estimated with BUSCO (insecta_odb9, v3.0.2)⁴⁸ and k-mer analysis and by mapping the initial reads back to the assembly.

Transcriptome sequencing and analysis

RNA from whole individuals of H. coronatus (first-, second-, and fifth-instar larvae and adults) was used for Nanopore full-length cDNA sequencing. Meanwhile, samples of the brain, thorax, abdomen, and mid and hind legs of each individual were sequenced using Illumina HiSeq 2500 sequencing. Fifth-instar larvae of D. lobata were used for Nanopore full-length cDNA sequencing and Illumina HiSeq 2500 sequencing. RNA from female and male subadults was extracted and subjected to Illumina HiSeq 2500 sequencing to investigate sex-related DEGs. At least three replicates were used for each group.

For H. coronatus, 222.31 Gb of Nanopore full-length cDNA sequencing data and 573.28 Gb of Illumina read cDNA sequencing data were obtained across different developmental stages and assembled into transcript datasets, yielding 13,817 transcripts (Supplementary Table 9). For D. lobata, 123.60 Gb of Nanopore full-length cDNA sequencing data and 68.21 Gb of Illumina read cDNA sequencing data were obtained from adult individuals and yielded 15,494 transcripts (Supplementary Table 9).

Relative gene expression was measured as transcripts per million reads, and DEGs were identified using DESeq2 (v1.0)⁴⁹ with false discovery rate (FDR) corrected p value (q-value) set to less than 0.05 and |log2-fold-change| set to greater than 1. Gene ontology functional enrichment and KEGG pathway analysis of the DEGs were performed with Goatools and KOBAS⁵⁰, respectively, with a Benjamini‒Hochberg-corrected p value less than 0.05 indicating statistical significance.

Genome annotation

Transposable elements were identified using RepeatMasker (open-4.0.7), RepeatModeler (v1.0.8), and MITE-Hunter (v1.0.8)⁵¹. Gene structures were determined by combining ab initio and homology methods. For ab initio annotation, we used Augustus (v3.2.1)⁵² and GENSCAN (v1.0)⁵³ to analyze the repeat-masked genome. For homolog-based annotation, protein sequences of the fruit fly (Drosophila melanogaster), cockroach (Blattella germanica), honeybee (Apis mellifera), mosquito (Aedes aegypti), and small brown planthopper (Laodelphax striatella) genomes were aligned to mantis genome sequences using BLAST software (v2.3.0)⁵⁴. Together with transcriptomic data, gene sets from these three methods were then integrated by EVidenceModeler software (v1.1.1)⁵⁵. For gene functional annotation, the integrated gene set was aligned against public databases, including KEGG, Swiss-Prot, TrEMBL, COG, and NR, with BLAST (v2.3.0)⁵⁴ and merged with annotations by InterProScan (v4.8) software⁵⁶. Annotation integrity was estimated by comparison with reference genome annotations and BUSCO (v3.0.2)⁴⁸, resulting in 95.5–98.2% completeness according to BUSCO analysis, suggesting the high quality of the annotation.

Identification of orthologous groups and phylogenomic analysis

To cluster families of protein-coding genes, we extracted protein sequences from the genomes of H. coronatus, D. lobata, and 16 other species of Insecta: Zootermopsis nevadensis (Blattodea: Termopsidae), B. germanica (Blattodea: Ectobiidae), D. melanogaster (Diptera: Drosophilidae), Locusta migratoria (Orthoptera: Acrididae), Clitarchus hookeri (Phasmatodea: Phasmatidae), Bombyx mori (Lepidoptera: Bombycidae), A. mellifera (Hymenoptera: Apidae), Ctenocephalides felis (Siphonaptera: Pulicidae), Stenopsyche tienmushanensis (Trichoptera: Stenopsychidae), Folsomia candida (Entomobryomorpha: Isotomidae), Ladona fulva (Odonata: Libellulidae), Pediculus humanus (Phthiraptera: Pediculidae), Photinus pyralis (Coleoptera: Lampyridae), Campodea augens (Diplura: Campodeidae), Rhopalosiphum maidis (Hemiptera: Aphididae), and Ephemera danica (Ephemeroptera: Ephemeridae). Protein sequences showing redundancy caused by alternative splicing variations or premature codons were removed from the protein-coding genes. The protein sequences were aligned reciprocally (i.e., all-vs.-all) using BLASTP programs with an E-value ≤ 1e − 5 and then clustered using OrthoMCL (v2.0.9)⁵⁷. Finally, 22,455 orthogroups were identified among the 18 species.

To reveal the phylogenetic relationships among Mantodea and other insects, the protein sequences of 324 1:1 orthologs from all 18 species were aligned with MUSCLE (v3.8.31)⁵⁸ and then concatenated using in-house Perl scripts. RAxML (v8.2.10)⁵⁹ was used to construct a phylogenetic tree for the superalignment using the GTRGAMMA model with Campodea augens and Folsomia candida as outgroups. We used the best-fitting substitution model as deduced by ProteinModelSelection.pl and 1000 replicates for bootstrap support. The MCMCTree program of the PAML (v4.8)⁶⁰ package was used to determine divergence times with the approximate likelihood calculation method and three dated fossil records (Baissatermes lapideus (145–99.6 Ma), Raphogla rubra (260–251 Ma), and Rhyniella praecursor (412.3–391.9 Ma)) from TimeTree (www.timetree).

Expansion and contraction of gene families

Family expansion/contraction analysis was performed by CAFÉ (v3.1)⁶¹ calculations with the parameters lambda -s and p < 0.01 based on the phylogenetic tree constructed above. Gene expansion and contraction results for each branch of the phylogenetic tree were obtained.

Demographic history reconstruction

To trace the demographic history of H. coronatus and D. lobata, we employed PSMC to estimate changes in effective population size using heterozygous sites²⁵, with the following set of parameters: −N 30 –t 15 –r 5 −p 4 + 25 × 2 + 4 + 6. H. coronatus and D. lobata take approximately 8 months to reach adulthood and then mature sexually after 25–35 days, so the generation time (g) was set to 0.75 years. The nucleotide mutation rate (μ) of H. coronatus was estimated to be 1.85 × 10^–9 mutations per site per generation with D. lobata as the reference species for comparison using the following formula: μ = D × g/2T⁶², where D is the observed frequency of pairwise differences between the two species, T is the estimated divergence time, and g is the estimated generation time for H. coronatus and D. lobata. Global mean surface temperature was estimated from benthic d¹⁸O⁶³.

Genome duplication

To investigate the genome evolution of H. coronatus and D. lobata, we searched for genome-wide duplications with B. germanica as an outgroup. We identified different modes of gene duplication as whole-genome duplications (WGDs), tandem duplicates (TDs), proximal duplicates (PDs, gene distance on the same chromosome less than 10), transposed duplicates (TRDs, transposed gene duplications), or dispersed duplicates (DSDs) using DupGen_finder⁶⁴ with default parameters. Five duplication categories were identified in the H. coronatus and D. lobata genomes, including DSD (53.59%), TD (23.12%), TRD (14.12%), PD (8.42%), and WGD (0.76%). Then, the overlap between the various types of duplicated genes and expanded genes was identified, yielding 796 tandem replication genes belonging to expanded gene families, including the Cuticle gene family.

To verify whether the identified tandem gene clusters existed, we mapped the Illumina and Nanopore reads back to these two genome assemblies and examined whether the sequencing depth of each tandem replicated gene was consistent with the sequencing depth of the whole genome. Finally, a uniform sequencing depth of all the duplicated genes consistent with the sequence depth of the whole genome was observed, confirming that these tandem gene clusters were real.

Expansion of the ABC, Cuticle, Trypsin, CYP450, and UGT gene families

Genome-wide protein sequences of H. coronatus, D. lobata, Z. nevadensis, B. germanica, and D. melanogaster were extracted, and the conserved nucleotide binding domain (PF00005.24) and transmembrane domain (PF00664.20) were scanned genome-wide for candidate ABC transporter genes using the Hidden Markov Model (HMM) in R (v3.2.1)⁵². To assign the candidate ABC genes to different subfamilies, multiple alignments of the ABC transporter protein sequences were performed using MUSCLE (v3.8.31)⁵⁸, and the poorly aligned regions and partial gaps were removed with trimAI (gt = 0.5). Then, the alignments were subjected to phylogenetic analysis using RAxML based on the maximum likelihood method with the following parameters: -f a -x 12345 -N 1000 -p 12345 -m PROTGAMMAJTTX. The resulting trees were displayed and edited using FigTree (v1.4.4, github.com/rambaut/figtree/releases). The subfamily assignment of ABC proteins in each species was further confirmed using BLASTP analyses on the NCBI webserver (www.ncbi.nlm.nih.gov/blast). In addition, the same analyses were performed to identify Cuticle (pfam:PF00379), Trypsin (pfam:PF00089), UGT (pfam:PF00201), and CYP450 (pfam:PF00067). We also performed manual annotation of these genes to avoid omissions in their general feature format (GFF) files.

Extraction of insect pigments and color reaction in vitro

To extract the pigment of the orchid mantises, second-instar individuals were crushed in liquid nitrogen and dissolved in acidified methanol (hydrochloric acid at 0.5%)⁶⁵, followed by incubation in a thermostatic shaker for 48 h at 220 rpm and 25 °C. After centrifugation, the supernatant was collected and passed through a 0.22-μm filter. The redox-dependent color changes of the pigments were observed by adding 10 μL of oxidant (NaNO₂ with a concentration of 1%) and reductant (ascorbic acid with a concentration of 1%). In acidified methanol, the redox states of the ommochrome pigments were stable, and they were stored at −20 °C. Five individuals were assigned to each group.

Transmission electron microscopy

The ultrastructure of the mid and hind legs was investigated by transmission electron microscopy. The legs were dissected from fifth- instar H. coronatus and fixed with 3% glutaraldehyde in 0.2 mol/L phosphate buffer (pH 7.2) for 48 h at 4 °C. The legs were then rinsed 3 times with phosphate buffer followed by postfixation in 1% osmium tetroxide for 3 h at 4 °C. The legs were washed several times, for 10 min each time, and placed into a series of ascending concentrations of acetone (50, 70, 80, 90, and 100%) for dehydration. The samples were embedded in Epon 812 for 2 h at room temperature. Then, the samples were trimmed to prepare ultrathin sections. The ultrastructure was captured by a JEM-1200EX transmission electron microscope (TEM, JEOL, Japan).

High-performance liquid chromatography analysis

Individuals of the first, second, and fifth instars were used to determine the concentrations of kynurenine and 3-hydroxykynurenine. The same amount of leg tissue was homogenized directly, added to a methanol solution containing the internal standard, homogenized again, and then centrifuged at 4 °C and 18,000 × g for 20 mins (Microfuge 20 R, Beckman Coulter, Inc., Indianapolis, IN, USA). The supernatant was used to detect the concentrations of kynurenine and 3-hydroxykynurenine according to the protocol of a previous study⁶⁶. The standards of kynurenine (CAS No. 343-65-7) and 3-hydroxykynurenine (CAS No. 484-78-6) were purchased from Sigma‒Aldrich. Five individuals were assigned to each group.

Histological analysis of the petal-like femoral lobes and size measurement of leg tissue

Hematoxylin-eosin staining was performed for structural observation of the leg tissue of H. coronatus. Briefly, the mid and hind legs of H. coronatus were fixed in 4% paraformaldehyde overnight at 4 °C. Then, paraffin-embedded, 4 μm tissue slices were cut and placed on slides for hematoxylin-eosin staining, with hematoxylin staining the nucleus and eosin staining the cytoplasm.

ImageJ2 (imagej.net/) was used to measure the length and area of the femur, tibia and tarsus of T2 and T3 legs of H. coronatus. Briefly, each segment of leg was dissected and taken photos. The scale bar of image was set followed by marking the target region by color threshold to automatically measure the area. The length of the target was measured by straight line.

RNA interference

To validate the function of Arm in the morphological development of legs, double-stranded RNAs (dsRNAs) were synthesized for each gene. The dsRNA sequences are shown in Supplementary Table 10. dsRNA targeting each gene was injected into the right midleg of first-instar individuals at 0.1 μL (660 ng/μL) using a microinjection system, and the second injection was conducted as described above at day 5 after the first injection. Control group individuals were injected with the same amount of PBS. Four individuals were assigned to each group. The mRNA expression levels of these genes were detected at day 3 after the second injection. One individual per group was cultured to molting for observation of the phenotypic outcome.

Real-time quantitative polymerase chain reaction (PCR)

Total RNA was extracted from the femur, tibia and tarsus of T2 and T3 legs using RNAprep Pure Micro Kit (DP420, TIANGEN Biotech (Beijing) Co.,Ltd). RNA quality was determined using a NanoDrop 2000 (Thermo Fisher Scientific Inc.), and RNA of suitable quantity was reverse transcribed using a PrimeScript RT Reagent Kit (RR037A, TaKaRa Biotechnology Co., Ltd.). Quantitative real-time PCR was performed using TB Green Premix Ex Taq II (Tli RNaseH Plus) (RR820A, TaKaRa Biotechnology Co., Ltd.) on an Mx3000P Real-Time PCR System (Agilent Technologies, Inc.). The relative expression level of the Ubx and Arm genes was quantified using GAPDH as the internal reference. Relative fold-change in gene expression was calculated using the delta-delta cycle threshold method⁶⁷. Three biological repeats were assigned to each group. Genes with a |fold-change| of greater than 2 and an adjusted p value less than 0.05 were identified as DEGs.

Enzyme activity of trypsin

Insect guts were isolated for Trypsin detection. Samples were immediately frozen in liquid nitrogen and stored at −80 °C before use. The enzyme activity of trypsin was detected using a commercial kit from Solarbio Technology Company (BC2315; Beijing, China). In this experiment, N-benzoyl-L-arginine-ethylester (BAEE) was used as the substrate. Under the catalysis of Trypsin, the ester bond of BAEE was hydrolyzed to produce a molecule of N-benzoyl-L-arginine (BA) and ethanol. The UV absorption of BA at 253 nm was much higher than that of BAEE. The amount of enzyme that increased by 0.001 per minute was calculated as one unit of activity. Four individuals were assigned to each group.

Statistical analysis

Measurements of continuous variables are expressed as the means ± SDs. Unless otherwise stated, differences between two groups were assessed using Student’s t tests, and differences among three or more groups were assessed using one-way analyses of variance in SPSS v18 (SPSS). Dunnett’s post hoc tests were used to compare treatment groups to controls. Adjusted p < 0.05 was considered statistically significant.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.