The squeaking silkmoths Rhodinia is a genus of wild silkmoths in the family Saturniidae of Lepidoptera that was first described by Otto Staudinger in 1892 (Kitching et al., Reference Kitching, Rougerie, Zwick, Hamilton, St Laurent, Naumann, Ballesteros and Kawahara2018). Up to date, at least 13 Rhodinia species have been recorded in both the subtropical and temperate Palaearctic regions of Asia (Ratnasingham and Hebert, Reference Ratnasingham and Hebert2020). The name of squeaking silkmoths came from their squeaking larvae; the caterpillars can squeak loudly to discourage potential predators from attacking them when touched (Coppens, Reference Coppens2020). Rhodinia fugax Butler 1877 (fig. 1) is one such species native to Northern China, Japan, Korea, and Russia. Rhodinia fugax is a univoltine insect that enters an egg diapause in the winter, and its larvae can feed on a wide variety of plants including Quercus and Castanea (Fagaceae), Pterocarya (Juglandaceae), Ulmus (Ulmaceae), and Salix viminalis (Salicaceae). This species spins a unique green cocoon that is shaped like a purse and has a small slit on the top that can split open which functions as an escape hole for the emerging moth (Coppens, Reference Coppens2020). In sericulture, Rhodinia species possess the potential significance (Mitamura, Reference Mitamura1998).
The life cycle of R. fugax. Photographed by Zheng-Yang Liu (Saturniidae@qq.com). (1–4) The first to fourth instar larvae. (5–7) The frontal, lateral, and ventral views of fifth instar larvae. (8–9) A cocooning larva and some cocoons. (10–11) The female moth and its antennae. (12) The male moth. (13–14) Mating moths and a moth laying eggs.
The mitochondrion is a fundamental eukaryotic organelle that functions in energy production. The mitochondrion has its genome that governs maternal inheritance and has several unique features such as faster evolution rate, low or absence of homologous recombination, evolutionary conserved gene products, and richness in genetic polymorphism (Cameron, Reference Cameron2014). Over the past few decades, insect mitochondrial genomes (mitogenomes) have been widely used as an informative marker for diverse evolutionary studies including molecular evolution, phylogenetics, and population genetics (Cameron, Reference Cameron2014). To date, nearly 600 Lepidopteran insect mitogenomes have been sequenced (https://www.ncbi.nlm.nih.gov/). In the majority of Lepidopteran insects, the mitogenome is a double-stranded, circular molecule, varying in size between 15 and 18 kb, and contains 37 genes including 13 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and an A + T-rich region (Zardoya, Reference Zardoya2020). As for Saturniidae species, 19 mitogenomes have been available in GenBank (table 1); however, no complete mitogenome is available for Rhodinia genus.
List of species considered for comparative mitogenomics study
In the present study, we report the complete mitogenome of R. fugax and compare the resulting mitogenome with other Saturniidae species. We also perform a phylogenetic analysis based on the available complete mitogenome sequences of Saturniidae to provide insight into the phylogenetic position of this wild silkmoth. The mitogenome of R. fugax presented in this study represents the first mitogenome from Rhodinia genus that is helpful for insight into the phylogenetic relationship of Saturniidae species.
Materials and methods
Insect, DNA extraction, and species identification
The wild larvae of R. fugax were collected at the Experimental Field of Sericultural Institute of Jilin Province (N43°47′31.45″; E126°32′28.74″), Jilin, China. They were then reared on oak trees (Quercus liaotungensis) in the field until cocooning. The partial body of one male pupa was used to extract total DNA. For species identification, the partial COI fragment corresponding to the DNA barcode was PCR amplified with insect universal primer COI-J1483 and COI-N2194 (Simon et al., Reference Simon, Buckley, Frati, Stewart and Beckenbach2006), and the resulting sequence was submitted to NCBI for blast search. The partial COI fragment of 658 bp corresponding to DNA barcode showed a 99% sequence identity with R. fugax shaanxiana voucher barcode SNB 1333 (GenBank accession no. GU702872) released by International Barcode of Life, indicating that the target insect belongs to R. fugax species.
Mitogenome amplification and DNA sequencing
The whole R. fugax mitogenome was amplified in four overlapping fragments with two universal primer pairs (Simon et al., Reference Simon, Buckley, Frati, Stewart and Beckenbach2006) and two specific primer pairs (table 2). In the first step, two fragments of less than 1300 bp were successfully amplified with insect universal primer pairs as described (Simon et al., Reference Simon, Buckley, Frati, Stewart and Beckenbach2006). The resulting PCR products were sequenced by the Sanger method with the PCR primers. Then, two long overlapping fragments were amplified with two specific primer pairs based on the obtained sequences of the two fragments with the universal primer pairs. The two long fragments are overlapped in COI and ND5-ND4 by 318 and 749 bp, respectively (table 2). PCR reaction for the two long fragments was performed in a 25 μl volume with 0.25 U of LA Taq (TaKaRa Co., Dalian, China), 1 μl (about 20 ng) of genomic DNA, 2.5 μl 10 × LA Taq buffer (Mg2+ plus), 100 μM dNTPs, and 10 pmol each primer. PCR amplification was performed under the following procedure: 96°C for 2 min, followed by 35 cycles of 10 s at 97°C and 9 min at 55°C, with a subsequent 10 min final extension at 72°C. The PCR products were resolved by electrophoresis in 1.0% agarose gel. After purification with TIANgel Midi Purification Kit (TIANGEN, Beijing, China), the two long PCR fragments were mixed equally and sequenced on Illumina Hiseq platform (PE150) by Frasergen Bioinformatics Co. (Wuhan, China).
Primers used for PCR amplification
F, forward; R, reverse.
Sequence assembly and genome annotation
A total of 2 G raw paired reads obtained was uploaded to Galaxy Europe (https://usegalaxy.eu/) (Jalili et al., Reference Jalili, Afgan, Gu, Clements, Blankenberg, Goecks, Taylor and Nekrutenko2020). Sequencing data were assembled using online assembly tools (Shovill 1.0.4, ABySS 2.2.3 and VelvetOptimiser 2.2.6) in Galaxy Europe (Jalili et al., Reference Jalili, Afgan, Gu, Clements, Blankenberg, Goecks, Taylor and Nekrutenko2020). The R. fugax mitogenome was completed through the assembly of two overlapping sequences via the alignment of neighboring fragments using Clustal X (Thompson et al., Reference Thompson, Gibson, Plewniak, Jeanmougin and Higgins1997). The minimum depth-of-coverage was ×100. The PCGs were determined by sequence comparison with the known Saturniidae mitogenomes available in GenBank. The 5′ ends of PCGs were assumed to be at the first legitimate in-frame start codon (ATN, GTG, TTG, GTT) in an open reading frame that was not located within an upstream gene encoded on the same strand, and the 3′ ends were inferred to be at the first in-frame stop codon encountered downstream of the start codon. Both tRNAscan-SE search server and MITOS web server were used to identify the tRNA genes and their secondary structures with default parameters (Lowe and Eddy, Reference Lowe and Eddy1997; Bernt et al., Reference Bernt, Donath, Jühling, Externbrink, Florentz, Fritzsch, Pütz, Middendorf and Stadler2013). The combination of these two methods could identify all of the tRNAs. The lrRNA (large ribosomal RNA) gene was annotated to extend to the boundaries of the flanking tRNALeu(CUN) and tRNAVal. The 3′ end of the srRNA (small ribosomal RNA) gene was annotated to be adjacent to the start of tRNAVal, while the 5′end was determined via comparing with orthologous sequences of other known Saturniidae mitogenomes. The entire A + T-rich region was subjected to a search for the tandem repeats using the Tandem Repeats Finder program (Benson, Reference Benson1999). The overlapping regions and intergenic spacers between genes were counted manually. Finally, the complete annotated file of the mitogenome was prepared using the NCBI Sequin tool (http://www.ncbi.nlm.nih.gov/Sequin/) and the sequin file was submitted to GenBank. The linear genome map was drawn by the program OGDRAW (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html) (Greiner et al., Reference Greiner, Lehwark and Bock2019).
Comparative genomics analysis
A total of seven Saturniidae mitogenomes were used to perform a comparative analysis with R. fugax (table 1). Comparative analysis was carried out to find out similarities and differences of R. fugax with other Saturniidae insects in terms of characteristics such as genome size, organization, structure, gene content, and arrangement. The nucleotide composition of sequences of the whole mitogenome, concatenated and individual PCGs, tRNAs, rRNAs, and control region was calculated using MEGA X (Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018). The base composition skewness was calculated for all the regions of mitogenome using formulas AT skewness = [A − T]/[A + T] and GC skewness = [G − C]/[G + C] (Perna and Kocher, Reference Perna and Kocher1995). Multiple sequence alignment was carried out using Clustal X (Thompson et al., Reference Thompson, Gibson, Plewniak, Jeanmougin and Higgins1997).
Phylogenetic analysis
A total of 25 mitogenomes were used to reconstruct the phylogenetic relationships within Bombycoidea (table 1). Three Bombycidae and two Sphingidae species were included (Yukuhiro et al., Reference Yukuhiro, Sezutsu, Itoh, Shimizu and Banno2002; Cameron and Whiting, Reference Cameron and Whiting2008; Pan et al., Reference Pan, Yu, Xia, Dai, Liu, Lu, Zhang and Xiang2008; Kim et al., Reference Kim, Choi and Kim2013; Kong and Yang, Reference Kong and Yang2015), and five non-Bombycoidea species Protantigius superans, Biston panterinaria, Phthonandria atrilineata, Spindasis takanonis, and Drosophila melanogaster served as the outgroups (de Bruijn, Reference de Bruijn1983; Yang et al., Reference Yang, Wei, Hong, Jiang and Wen2009, Reference Yang, Xue and Han2013; Kim et al., Reference Kim, Kang, Jeong, Kim and Kim2011). Each PCG was aligned on the RevTrans 2.0b Server (http://www.cbs.dtu.dk/services/RevTrans-2.0/web/) using the MAFFT algorithm and the invertebrate mitochondrial genetic code (Wernersson and Pedersen, Reference Wernersson and Pedersen2003). The well-aligned blocks were selected using GBlocks 0.91b software and concatenated to generate PCG alignment (Castresana, Reference Castresana2000). Substitution model selection was also conducted based on the lowest BIC scores (Bayesian Information Criterion) using MEGA X (Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018). Bayesian inference (BI) and maximum-likelihood (ML) methods were conducted using MrBayes ver. 3.1.2 and IQ-TREE 1.6.12, respectively, under the GTR + G + I nucleotide substitution model (Huelsenbeck and Ronquist, Reference Huelsenbeck and Ronquist2001; Nguyen et al., Reference Nguyen, Schmidt, von Haeseler and Minh2015). For the BI analysis using Markov Chains and Monte Carlo methods, two independent runs of four incrementally heated chains (one cold chain and three hot chains) were simultaneously run for one million generations, with tree sampling conducted at every 100 generations; the first 25% of the sampled trees were discarded as burn-in. An average split frequency of <0.01 was the decision criterion for convergence of the two simultaneous runs. The confidence values used for the BI tree were the Bayesian posterior probabilities. Support for the inferred ML tree was inferred by bootstrapping with 10,000 replicates.
Results and discussion
General characterization of R. fugax mitogenome
The entire R. fugax mitogenome is 15,334 bp long, well within the range detected in Saturniidae species, ranging from 15,236 bp in Actias selene (Liu et al., Reference Liu, Zhu, Dai, Wei and Liu2012) to 15,575 bp in Antheraea proylei (Yang et al., Reference Yang, Zhang, Chen, Chen, Li and Liu2019). The mitogenome presents a typical gene content observed in metazoan mitogenomes: 13 PCGs, 22 tRNA genes, two rRNAs, and a non-coding A + T-rich region (fig. 2 and table 3). In the major strand, nine PCGs and 14 tRNAs are encoded, whereas in the minor strand, four PCGs, eight tRNAs, and two rRNA genes are encoded. The gene order and orientation of this mitogenome are identical to the completely sequenced Saturniidae species.
Mitogenome arrangement of R. fugax with respect to other selected species. tRNAs are denoted as one-letter symbols according to the IUPAC-IUB single-letter amino acid codes. The one-letter symbols L, L*, S, and S* denote tRNALeu(CUN), tRNALeu(UUR), tRNASer(AGN), and tRNASer(UCN), respectively.
Gene organization of the R. fugax mitogenome
F, forward; R, reverse. IGS, intergenic spacer region; OL, overlapping region; T-tRNA, incomplete termination codon T adjacent to tRNAs.
The nucleotide composition of the whole mitogenome exhibits an 80.57% AT content and 19.43% GC content. The enhanced AT content has been observed in many lepidopteran insects (Kim et al., Reference Kim, Kang, Jeong, Kim and Kim2011; Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017). The A + T-rich region exhibits the maximal AT content (91.62%), followed by rRNAs, tRNAs, and PCGs (table 4). The individual PCGs also exhibit variation in the corresponding AT content. For instance, ATP8, ND6, and ND4L had higher AT content while COI and COIII had lower values. This variation pattern has been detected in the other Saturniidae insects, such as Eriogyna pyretorum and Antheraea assamensis (Jiang et al., Reference Jiang, Hong, Yu, Li, Yang, Liu and Wei2009; Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017). The AT skewness of the major strand of this mitogenome is slightly negative (−0.007), as found in Samia cynthia ricini (−0.006), E. pyretorum (−0.031), A. assamensis (−0.018), Antheraea pernyi (−0.021), A. selene (−0.023), Gonimbrasia belina (−0.008). However, the AT skewness is slightly positive in Cricula trifenestrata (0.011).
Nucleotide composition of sequences of whole mitochondrial genome, PCGs, tRNAs, rRNAs, and the A + T-rich region
The mitogenome of R. fugax contains a total of 138 bp intergenic spacer sequences, which are spread over 16 regions ranging in size from 1 to 54 bp. The largest intergenic spacer sequence is present between tRNAGln and ND2, with an A + T content of 94.44%, as found in many other Saturniidae insects such as S. cynthia ricini, E. pyretorum, and A. selene (Jiang et al., Reference Jiang, Hong, Yu, Li, Yang, Liu and Wei2009; Kim et al., Reference Kim, Park, Kim, Kang, Kim, Jin, Han and Kim2012; Liu et al., Reference Liu, Zhu, Dai, Wei and Liu2012). This region displays limited sequence conservation between even closely related lepidopteran groups such as within Bombycoidea or between bombycoids (Cameron and Whiting, Reference Cameron and Whiting2008), indicating that it would imply no functional significance or might not serve as another origin of replication (Hong et al., Reference Hong, Lee, Jo, Park, Kim, Hwang, Jin, Kang, Kim, Han and Kim2008). In the spacer tRNASer(AGN)-ND1, a conserved motif ‘ATACTAA’ is also found in the R. fugax mitogenome, which has been found to be present in most insect mitogenomes even if the nucleotide sequence can be quite divergent (Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017). This 7 bp motif is implicated in mitochondrial transcription termination where the ‘ATACTAA’ motif is essential as a recognition site for mtTERM (Cameron and Whiting, Reference Cameron and Whiting2008).
This genome has a total of 23 bp overlapping sequences varying from 1 to 8 bp in size in seven regions. The longest overlap is 8 bp long, located between tRNATrp and tRNACys, similarly sized overlaps are also observed in sequenced Saturniidae species S. cynthia ricini, E. pyretorum, and A. assamensis (Jiang et al., Reference Jiang, Hong, Yu, Li, Yang, Liu and Wei2009; Kim et al., Reference Kim, Park, Kim, Kang, Kim, Jin, Han and Kim2012; Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017). The 7 bp overlap with the reading frame involving the ATP8/ATP6 genes is also found, which is common to other sequenced Saturniidae mitogenomes (Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017).
The 13 PCGs in the mitogenome of R. fugax span 11,226 bp constituting around 73.21% of the total mitogenome. The putative start codons of PCGs are those previously known for metazoan mitogenomes with the typical ATN codon (six with ATG, three with ATT, two with ATA, and one with ATC), with the only exception represented by the CGA start codon observed in COI. Further, three genes COI, COII, and ND5 show incomplete termination codon (only T) instead of typical termination codon pattern TAA or TAG (table 3). Such an incomplete termination codon can become complete (TAA) by posttranslational modifications such as polyadenylation that occur during the mRNA maturation process (Taanman, Reference Taanman1999).
The typical set of 22 tRNA genes are interspersed in the mitogenome of R. fugax, ranging in size from 64 to 71 nucleotides. The anticodons of 22 tRNA genes are all identical to their counterpart tRNA genes of the sequenced lepidopterans (Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017). All tRNAs present the typical clover-leaf secondary structure except for tRNASer(AGN) and tRNAAsn. The R. fugax tRNASer(AGN) lacks a stable DHU arm (Supplementary fig. 1), as observed in other lepidopteran species (Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017). The tRNAAsn lacks a stable TΨC arm, exhibiting a unique feature among lepidopteran mtDNAs. A total of 27 mismatched base pairs occur in the R. fugax mitochondrial tRNA genes, 20 of them are G–U pairs, which form a weak bond. The mismatches are scattered among 18 of the 22 tRNA genes (table 5), leading to changes in tRNA structures. Mismatches observed in tRNAs can be corrected through RNA-editing mechanisms that are well known for arthropod mtDNA (Lavrov et al., Reference Lavrov, Brown and Boore2000).
The mismatches in the tRNAs of R. fugax mitogenome
Characterization of R. fugax A + T-rich region
The A + T-rich region of R. fugax spans 370 bp and have several features common to the Bombycidea insects (fig. 3), including the motif ATAGA followed by an 18 bp poly-T stretch, microsatellite AT repeat sequence preceded by the ATTTA motif, and a poly-A element immediately upstream tRNAMet (Sima et al., Reference Sima, Chen, Yao, Li, Liu, Jin, Wang, Su, Li and Liu2013). The ATAGA motif has been suggested to be the precise position of the origin of replication, and the poly-T sequence has been suggested to function as a structural signal for the recognition of proteins in the initiation of replication for minor-strand mtDNA in Bombyx mori (Saito et al., Reference Saito, Tamura and Aotsuka2005). An 18 bp poly-A element immediately upstream tRNAMet and intermittently interrupted by T was also found (fig. 3). The microsatellite AT repeat sequence preceded by the ATTTA motif often seen in other sequenced Bombycidea insects is not found in the R. fugax A + T-rich region. The presence of multiple tandem repeat elements is a characteristic of the insect A + T-rich region. For example, a repeat element of 38 bp has been identified in the A + T-rich region of A. pernyi and Antheraea roylei (Arunkumar et al., Reference Arunkumar, Metta and Nagaraju2006; Liu et al., Reference Liu, Li, Pan, Dai, Zhu, Lu and Xiang2008), and another repeated element of 126 bp has been reported in the A + T-rich region of B. mori and Bombyx mandarina (Yukuhiro et al., Reference Yukuhiro, Sezutsu, Itoh, Shimizu and Banno2002; Pan et al., Reference Pan, Yu, Xia, Dai, Liu, Lu, Zhang and Xiang2008). However, such tandem repeats were not detected in the A + T-rich region of R. fugax mitogenome.
The features present in the A + T-rich region of R. fugax. The bolded and colored nucleotides sequentially indicate the ATATG motif (red), poly-T stretch (blue), ATTTA sequence (green), and poly-A stretch (purple).
Phylogenetic relationships
To date, the phylogenetic relationship of Rhodinia species remains uncertain. Bouvier (Reference Bouvier1936) divided Saturniinae into section Saturnitae and section Attacitae that consisted of tribe Attaciicae and tribe Rhodiicae, and Rhodinia and Pararhodia were treated as the genera of the tribe Rhodiicae (Bouvier, Reference Bouvier1936). In Michener's (Reference Michener1952) influential revision, the tribe Rhodiicae including Rhodinia and Pararhodia was reclassified into Saturniini (Michener, Reference Michener1952). Some researchers thought that Rhodinia had a closer relationship with Attacini than Saturniini and even belonged to Attacini (Regier et al., Reference Regier, Grant, Mitter, Cook, Peigler and Rougerie2008). However, the NCBI taxonomy browser (https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/), Encyclopedia of Life (https://eol.org/), Wilispecies (https://species.m.wikimedia.org/wiki/Rhodinia/), and iNaturalist (https://www.inaturalist.org/) still place R. fugax into Saturniini tribe. In the present study, BI tree and ML tree were reconstructed using the nucleotide sequences of 13 PCGs from 25 mitogenomes (fig. 4). The resulting tree topologies were congruent and nodal support values varied little for the two different analyses. In two phylogenetic trees, the used Saturniidae species were split into Saturniini, Attacini, and Bunaeini, as found based on the traditional morphological classification and previous molecular works (Bouvier, Reference Bouvier1936; Regier et al., Reference Regier, Grant, Mitter, Cook, Peigler and Rougerie2008; Jiang et al., Reference Jiang, Hong, Yu, Li, Yang, Liu and Wei2009; Chen et al., Reference Chen, Li, Chen, Wang, Li, Xia, Zeng, Li, Liu and Qin2014; Singh et al., Reference Singh, Kabiraj, Sharma, Chetia, Mosahari, Neog and Bora2017; Kim et al., Reference Kim, Kim, Jeong and Kim2018; Rubin et al., Reference Rubin, Hamilton, McClure, Chadwell, Kawahara and Barber2018; Wang et al., Reference Wang, Chen, Gu, Wang and Zwick2019; Langley et al., Reference Langley, van der Westhuizen, Morland and van Asch2020). Our phylogenetic analyses indicated R. fugax belongs to Attacini, instead of Saturniini, consistent with the morphological data and some previous studies (Bouvier, Reference Bouvier1936; Rubin et al., Reference Rubin, Hamilton, McClure, Chadwell, Kawahara and Barber2018). Although Neoris haraldi is grouped into the tribe Saturniini (http://www.sp2000.org), it has been recovered as basal to Attacini and Saturniini in a previous study (He et al., Reference He, Wang and Chen2017); our phylogenetic reconstructions within the family Saturniidae further support it as basal to Attacini, Saturniini, and Bunaeini.
Phylogeny of the superfamily Bombycoidea based on 13 PCGs. The tree with maximum likelihood (ML) method was shown. The values near each node specify bootstrap percentages of 10,000 replicates by ML analysis and Bayesian posterior probabilities in percent by Bayesian inference (BI) analysis, respectively.
Conclusion
This study for the first time reports the mitogenome of R. fugax and infers the phylogenetic relationships within Saturniidae based on available mitogenomes. Our results showed that R. fugax harbors an identical genomic organization with known Saturniidae mitogenomes. Our findings indicated that the A + T-rich region of R. fugax contains no tandem repeat elements and the microsatellite AT repeat sequence preceded by the ATTTA motif. The phylogenetic analyses based mitogenome showed that R. fugax belongs to Attacini instead of Saturniini. To our knowledge, this mitogenome presented in this study is the first mitogenome from Rhodinia genus that will provide insights into the phylogenetic relationship within Saturniidae.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485321000808
Acknowledgements
We like to thank Zheng-Yang Liu (Saturiidae@qq.com) for taking photographs of Rhodinia fugax. This work was supported by the National Natural Science Foundation of China (Grant no: 31972617) and the China Agriculture Research System (CARS-18).
Conflict of interest
None.
