Mitogenomics and phylogenetics of twelve species of African Saturniidae (Lepidoptera)

Specimen collection, morphological identification, DNA extraction

Adult and caterpillar specimens were collected at several locations in South Africa and Namibia, between January and March 2020 (Fig. 1, Table 1) under permits approved by the Ezemvelo KZN Permits Office (OP408/2020), Cape Nature (CN44-59-13497) and National Commission on Research Science and Technology (AN20190911). Specimens were euthanized by freezing within a few hours of collection, whenever field conditions allowed. Species identification of adults and caterpillars was made by R. Oberprieler based on photographs of reference specimens using current literature (Pinhey, 1972; Oberprieler, 1995; Staude et al., 2016; Staude et al., 2020). Legs from each adult and caterpillar were excised and stored in 100% ethanol at −20 °C until DNA extraction. Total DNA was extracted from one leg from each adult and caterpillar using a standard phenol-chloroform protocol (Sambrook & Russell, 2012).

DNA barcoding

The standard COI barcoding region was amplified using the primer pair LepF/LepR in a total reaction volume of 5 µL containing 2.5 µL of Qiagen Multiplex PCR Kit (QIAGEN), 0.5 µL of each primer, 0.5 µL of Milli-Q water and 1 µL of template DNA. PCR amplifications were performed with initial denaturation at 95 °C for 15 min; 35 cycles of 94 °C for 30 s, 50 °C for 90 s; followed by a final extension at 72 °C for 10 min. Sequencing reactions were performed using LepR with the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA, USA), at the Central Analytical Facilities of Stellenbosch University. The sequences were queried against the BOLD Systems (https://www.boldsystems.org/) and GenBank (https://blast.ncbi.nlm.nih.gov/) for confirmation of species identification on February 23, 2022.

Sequencing, assembly, and annotation of mitogenomes

Specimens were individually sequenced using the Ion Torrent™ S5™ platform (ThermoFisher Scientific, Waltham, MA, USA) available at the Central Analytical Facilities of Stellenbosch University, South Africa. Sequence libraries were prepared using the Ion Xpress™ Plus gDNA Fragment Library Kit (ThermoFisher Scientific, Waltham, MA, USA), according to the protocol MAN0009847, REV J.0. Libraries were pooled and sequenced using the Ion 540™ Chef Kit (ThermoFisher Scientific, Waltham, MA, USA). The Ion Torrent reads were quality trimmed with a 30-base sliding-window at an average threshold of Q16. The remaining sequencing reads were then filtered for read length, where any reads shorter than 25 bases were removed from the data. Reference-based assembly of the next-generation sequencing (NGS) reads was performed using mitogenomes of the closest Saturniidae relatives available on GenBank (Table S1). NGS reads were assembled and mapped to the reference sequences on Geneious Prime v2021.1.1 (https://www.geneious.com/) using the “Map to reference” function, under the low to medium sensitivity option and fine-tuning up to five iterations. The open reading frames of the expected 13 mitochondrial protein-coding genes (PCGs) were identified with Geneious Prime using the invertebrate mitochondrial genetic code. Transfer RNA (tRNA) genes were located with ARWEN software (http://130.235.244.92/ARWEN/) (Laslett & Canbäck, 2008), under the default composite metazoan mitochondrial genetic code. The two ribosomal RNA (rRNA) genes and the non-coding AT-rich region were annotated by manual comparison to other Saturniidae mitogenomes available on GenBank. The raw sequence data were deposited on GenBank (PRJNA796275), as well as the 12 assembled and annotated mitogenomes (OL912807 to OL912817/ and OL912953).

Mitogenomic and phylogenetic analyses

Mitogenome nucleotide composition and biases [(AT-skew = (A−T)/(A+T) and GC-skew = (G−C)/(G + C)] were calculated in Geneious Prime. Gene overlapping regions and intergenic spaces were counted manually. Non-synonymous (Ka) and synonymous (Ks) substitution rates in PCGs were calculated in DnaSP6 (Rozas et al., 2017), and relative synonymous codon usage was calculated in MEGA X (Kumar et al., 2018), with the invertebrate mitochondrial genetic code in both cases. The presence of tandem repeat elements in the AT-rich region was searched using the Repeat Finder plug-in available on Geneious Prime.

The phylogenetic positions of the 12 African Saturniidae species were inferred in the context of other sequences available on GenBank for the family as of October 2021 (Table S2). The final dataset (n = 44) included 35 species in 19 genera and five tribes, with Bombycidae (Bombyx mori, Bombyx mandarina and Rondotia menciana), Sphingidae (Manduca sexta and Sphinx morio) and non-Bombycoidea (Biston panterinaria, Phthonandria atrilineata, Protantigius superans and Spindasis takanonis) as outgroups. Phylogenetic analyses were performed using the 13 PCGs (all codon positions). Individual PCG sequences were extracted from each mitogenome and stop codons were removed manually. Translation alignments were performed separately for each PCG using the MAFFT algorithm in Geneious Prime and then concatenated in a single alignment for each specimen. Poorly aligned regions and alignment gaps were removed using GBlocks v0.91b (Castresana, 2000). Two Bayesian inference (BI) and one maximum likelihood (ML) methods were used for phylogenetic reconstruction. The best partitioning scheme and model for MrBayes v3.1.2 (Huelsenbeck & Ronquist, 2001) and ML were determined using the edge-linked greedy strategy (Lanfear et al., 2012) in PartitionFinder v2.1.1 (Lanfear et al., 2017) built-in Anaconda programme (https://repo.anaconda.com/), and GTR+I+G was considered the best model. BI in MrBayes was performed under the GTR+I+G nucleotide substitution model selected with PartitionFinder. Analyses were performed as follows: two independent runs of four heat chains, ten million generations run simultaneously, resampling every 1000 generations. The first 25% of trees were discarded as burn-in, and the decision criterion for the convergence of the two runs was set as an average split frequency of ≤ 0.01. BI analyses were also performed under the site-heterogeneous mixture model CAT+GTR in PhyloBayes MPI in XSEDE v1.8c (Lartillot et al., 2013) to minimize the effect of mitochondrial compositional heterogeneity on phylogenetic reconstructions (Lartillot, Brinkmann & Philippe, 2007; Cai et al., 2020). Constant sites were removed from the alignment, and the minimum number of cycles was set to 30,000 with the burn-in set to 1000. The “maxdiff” was set to 0.3, and the minimum effective size was set to 50. Nodal support in the MrBayes and PhyloBayes trees was estimated as Bayesian posterior probabilities (BPP). ML analyses were performed in IQ-TREE v1.6.12 (Nguyen et al., 2015) under the GTR+I+G substitution model selected by PartitionFinder. Branch supports were determined using 1000 replicates for both the ultrafast bootstrapping (UFBoot) and the SH-aLRT branch test (Guindon et al., 2010; Hoang et al., 2018). MrBayes, PhyloBayes and IQ-TREE analyses were run on the CIPRES Science Gateway Portal (Miller, Pfeiffer & Schwartz, 2010). The final trees were drawn using FigTree v1.4.4 (https://bio.tools/FigTree).

DNA barcoding

The sequence queries against BOLD Systems confirmed the morphological identification for all specimens used for mitogenome sequencing with sequence similarity between 99.44 and 100% (Table S3). The BLASTn queries against GenBank produced no matches with high identity (>98%) except for B. alcinoe and N. cytherea due to the absence of sequences for these species on the database.

Comparative mitogenomics of Saturniidae

The Ion Torrent runs generated between 11,321,643 (V. ducalis) and 21,284,325 (E. bauhiniae) reads, with an average size of 184 bp. Total mapped reads ranged from 16,637 (H. smilax) to 191,295 (H. dyops). Average sequence coverage varied between 197× (H. smilax) and 1810× (H. dyops) (Table 2). The total length of the new mitogenomes ranged from 15,218 bp (V. ducalis) to 15,363 bp (N. cytherea), in line with those of other Saturniidae (Kim et al., 2017a; Kim et al., 2017b; Langley et al., 2020; Liu et al., 2021; Chen et al., 2021). All mitogenomes generated in this study are identical in gene content and organization and include the 13 PCGs, two rRNA genes and 22 tRNA genes typical of Metazoa. Nine PCGs (ATP6, ATP8, CYTB, COI, COII, COIII, ND2, ND3 and ND6) and 14 tRNAs are located on the major (J) strand, and four PCGs (ND1, ND4, ND4L and ND5), eight tRNAs and the two rRNAs are located on the minor (N) strand (Table S4). Gene order is conserved across all species and identical to that of other Saturniidae, with the reciprocal translocation of tRNA^Met and tRNA^Gln relative to tRNA^Ile (M-I-Q) as the only difference to the hypothetical ancestral arrangement for the non-Ditrysia lineage Hepialoidea of Lepidoptera (I-Q-M) (Cao et al., 2012) (Fig. 2).

All tRNAs have the typical cloverleaf structure except for tRNA^Ser1 in B. alcinoe, E. bauhiniae, G. belina, G. tyrrhea, G. maja, H. apollonia, H. dyops, H. smilax, L. delegorguei, V. ducalis and V. grimmia, in that the dihydrouridine (DHU) arm is absent, as occurs in many Metazoa (Cameron, 2014) (Fig. 3). In contrast, the DHU arm of tRNA^Ser1 is present in E. trimenii, N. cytherea and N. wahlbergii. The length of the tRNAs ranges from 57 bp (tRNA^Phe) in N. wahlbergii to 74 bp (tRNA^Trp) in H. apollonia. The location and average size of 16S rRNA (1,357 bp; between tRNA^Leu1 and tRNA^{V al}) and 12S rRNA (760 bp; be-tween tRNA^{V al} and AT-rich region) across the 12 species is in line with the average size and position of the two genes in other Saturniidae. No tandem repeat elements were detected in the AT-rich regions of the 12 species.

The new sequences have the high AT content characteristic of insect mitogenomes, ranging from 78.80% in H. apollonia to 81.00% in E. trimenii (Table S5). The AT content of the total PCGs follows the same trend and varies from 77.40% in H. apollonia to 79.90% in E. trimenii. Average AT content across the 12 species is highest in ATP8 (92.70%) and lowest in COI (69.30%) (Fig. 4), in line with previous reports on Saturniidae (Jiang et al., 2009; Singh et al., 2017; Chen et al., 2021). AT-skew is negative for PCGs on the J-strand across all species and varies from −0.22 in ND3 to −0.02 in ATP8. All four PCGs on the N-strand have a positive AT-skew, ranging from 0.14 in ND5 to 0.24 in ND1 and ND4L. Additionally, most PCGs have a slightly negative GC-skew, ranging from −0.69 (ATP8) to −0.02 (COI), except for COIII = 0.00 in L. delegorguei.

All 13 PCGs initiate with standard ATN start codons except COI and COII, which start with the alternative codons CGA and GTG, respectively. CGA as the start codon for COI has been reported to be highly conserved in Lepidoptera (Margam et al., 2011), and GTG was reported as the start codon for COII in A. assamensis, E. pyretorum, S. boisduvalii and S. ricini (Jiang et al., 2009; Kim et al., 2012; Singh et al., 2017; Kim et al., 2017a; Kim et al., 2017b). PCGs terminate with complete TAA or TAG codons or with incomplete TA- or T-, which are presumed to be completed by posttranscriptional modifications such as polyadenylation (Salinas-Giegé, Giegé & Giegé, 2015) (Fig. 5, Table S6). The most frequent amino-acid codons are Leu, Ile, Phe and Ser across all species, whereas Cys, Asp, Arg, Glu and Gln are rare (Fig. 6A). The most frequently used codons are AT-rich, and this feature is reflected in the total mitogenome nucleotide bias towards A and T. ATP8 has the largest number of different codons compared to all other genes across the 12 species (Fig. 6B). Relative synonymous codon usage (RSCU) is higher than 1.0 for all codons and highest for Leu1 (Fig. 7). Average Ka/Ks was found to be less than 1.0 in all PCGs across all species, indicating purifying or stabilizing selection, and to be highest for ATP8 (Ka/Ks = 0.31) (Fig. 8).

Phylogenetic position of African Saturniidae in the family

The phylogeny of the family Saturniidae has been reconstructed using morphological, nuclear and mitochondrial sequence data. Regier et al. (2008) analysed four protein-coding nuclear genes using parsimony and ML methods and included five tribes present in Africa (Attacini, Bunaeini, Micragonini, Saturniini and Urotini). The recovered tribal structure strongly supported Attacini as sister group of Saturniini, which together formed the sister group of Urotini, Bunaeini and Micragonini, although with Urotini paraphyletic in regard to Bunaeini and Micragonini. More recently, the phylogeny of Saturniidae was reconstructed using a phylomorphology approach based on anchored-hybrid-enriched (AHE) loci followed by geometric morphometrics of hindwings (Rubin et al., 2018). The tree also supported a sister relationship between Attacini and Saturniini, the African tribes Bunaeini, Urotini and Micragonini together forming a monophylum and the sister group of Attacini+Saturniini, with Urotini again paraphyletic in regard to Bunaeini. In addition to phylogeny estimates based on nuclear genes and phylomorphology, the relationships among Saturniidae have been reconstructed using mitogenomic data (Sima et al., 2013; Kim et al., 2018; He, Wang & Chen, 2017; Zhong et al., 2017; Langley et al., 2020; Chen et al., 2021; Liu et al., 2020; Liu et al., 2021; Zhao et al., 2021). Although these studies focused on Asian species in the tribes Attacini and Saturniini, reciprocal monophyly and sister relationships were largely consistent, except for the positions of Cricula trifenestrata and Neoris haraldi, which were not consistently recovered, likely due to low sequence quality. The addition of the African tribe Bunaeini to the mitochondrial phylogeny did not challenge the sister relationship between Attacini and Saturniini, which together formed the sister group of Bunaeini (Langley et al., 2020; Chen et al., 2021), in agreement with the non-mitochondrial phylogeny estimates (Regier et al., 2008; Rubin et al., 2018). Thus, mitochondrial phylogeny estimates of Saturniidae have included only the tribes Attacini, Bunaeini and Saturniini, and African species were severely under-represented compared to Asian species (Langley et al., 2020; Chen et al., 2021). To contribute to the filling of this gap, we assessed the phylogenetic position of 12 African species of Saturniidae including the first representatives of the tribes Eochroini and Micragonini. Due to the suspected low quality of the mitogenomes of C. trifenestrata and N. haraldi, we performed phylogeny estimates with and without these sequences.

Phylogenetic trees excluding C. trifenestrata and N. haraldi

The ML and MrBayes trees excluding C. trifenestrata and N. haraldi (Fig. 9) displayed the same phylogenetic structure, recovering Bunaeini as a basal lineage, a sister relationship between Micragonini and Eochroini and the sister relationship between Saturniini and Attacini as found in previous analyses based on nuclear genes (Regier et al., 2008), phylomorphology (Rubin et al., 2018) and mitogenome data (Langley et al., 2020; Chen et al., 2021) (Fig. 9). In contrast, the PhyloBayes tree (Fig. 10) recovered a topology different from those of ML and MrBayes and of previous phylogenetic hypotheses in that Saturniini appear non-monophyletic, and several nodes have lower statistical support. Therefore, our discussion will focus on the topology recovered by the ML and MrBayes methods.

Eochroa trimenii, traditionally classified in Bunaeini (e.g., Bouvier, 1936; Pinhey, 1972; Staude et al., 2016; Staude et al., 2020), was not recovered in this tribe but as sister group of Micragonini, albeit with full nodal support only in the MrBayes tree. Its exclusion from Bunaeini had earlier been proposed by Oberprieler (1997) and was formalized by Cooper (2002), who placed it in a separate tribe, Eochroini. This result cannot be compared with previous mitochondrial and non-mitochondrial phylogeny estimates, because Eochroa is here included in a phylogenetic analysis for the first time. Although our result does not discount the possible inclusion of Eochroa in Micragonini, such a placement is repudiated by the fact that Eochroa does not possess the synapomorphic characters of Micragonini, as outlined by Oberprieler & Nässig (1994) and Oberprieler (1997). Furthermore, its relationship to the genus Usta, which appears in a similar position as sister group of Micragonini in the analyses of Regier et al. (2008) and Rubin et al. (2018), needs to be investigated. The ML and MrBayes trees recovered the monophyly of the tribes Bunaeini and Micragonini with full nodal support, and the results were broadly in agreement with previous phylogeny estimates (Rubin et al., 2018; Chen et al., 2021).

ML and MrBayes outperformed PhyloBayes in that the broad tribal structure of the first two methods is congruent with those of non-mitochondrial phylogeny estimates (Regier et al., 2008; Barber et al., 2015; Rubin et al., 2018), except that the tribes Micragonini and Eochroini are not recovered as sister group of Bunaeini but of Attacini+Saturniini instead. Different phylogenetic reconstruction methods (i.e., Bayesian inference vs ML) may impact tree topology, especially with regards to the order of deeper nodes. For example, a study of the high-level phylogeny of the Coleoptera inferred with mitogenome sequences showed that different inference models (ML and Bayesian inference) can yield inconsistent topologies for the same data (Yuan et al., 2016). Additionally, our analyses did not include all African Saturniidae lineages (missing in particular the various genera currently classified in Urotini), and this incomplete sampling may also not be sufficient for a well-supported recovery of the tribal structure of the family. For example, mitogenomic phylogeny estimates of Lepidoptera (Timmermans, Lees & Simonsen, 2014) and Hymenoptera (Mao, Gibson & Dowton, 2015) have failed to resolve some of the relationships with confidence, likely due to limited taxon sampling. The inclusion of representatives of Urotini will allow further insights into the mitochondrial phylogenetic reconstruction of Saturniidae, as this tribe was non-monophyletic in the phylogeny estimates based on nuclear genes (Regier et al., 2008) and phylomorphology (Rubin et al., 2018).

Phylogenetic trees including C. trifenestrata and N. haraldi

The phylogenetic analyses including C. trifenestrata and N. haraldi (Figs. S1–S2) resulted in lower support for the monophyly of Saturniini, and N. haraldi was recovered in a basal position relative to all other Saturniidae. The genus Neoris was not represented in the phylogeny estimates derived from nuclear genes (Regier et al., 2008) or phylomorphology (Rubin et al., 2018); therefore, it is not possible to establish comparisons. The monophyly of Attacini had low support in the ML tree with Rhodinia fugax, traditionally classified in Saturniini but morphologically similar to Attacini, as the basal-most member of Attacini, but the Bayesian trees more strongly and consistently supported the monophyly of this tribe including R. fugax. Cricula was recovered as sister group of Antheraea in non-mitochondrial phylogeny estimates (Regier et al., 2008; Rubin et al., 2018), but the position of C. trifenestrata was recovered inconsistently in our trees as well as in a previous work, where it formed a group diverged from the other Saturniini along with N. haraldi (Langley et al., 2020). In the present study, we detected several issues in the mitogenomes of C. trifenestrata and N. haraldi, including cases of shifts in the reading frame caused by artefactual single nucleotide indels that we corrected, but it is possible that other sequencing errors evaded our curation. Therefore, the positions of N. haraldi and C. trifenestrata seem be the result of substandard sequencing quality of these publicly available mitogenomes.