Mitogenomes of Three Satyrid Butterﬂy Species (Nymphalidae: Lepidoptera) and Reconstructed Phylogeny of Satyrinae

: Satyrinae is a 3000-species butterﬂy subfamily of Nymphalidae. The higher-level classiﬁcation of this family is still controversial. In this research, we sequenced the complete mitogenomes of three satyrid butterﬂy species, Hipparchia autonoe , Paroeneis palaearctica , and Oeneis buddha , and studied the phylogeny of Satyrinae with all known complete mitogenomes. The results showed that the lengths of the three satyrid butterﬂy mitogenomes are 15,435 bp ( H. Autonoe ), 15,942 bp ( P. palaearctica ), and 15,259 bp ( O. buddha ). Gene content and arrangement of newly sequenced mitogenomes are highly conserved and are typical of Lepidoptera. These three mitogenomes were found to have a typical set of 37 genes and an A + T-rich region. The tRNA genes in these three mitogenomes showed a typical clover leaf structure, but the stem of tRNASer (AGN) was lacking dihydroacridine. In these three species, the lengths of the A + T-rich regions were different, which led to differences in mitochondrial genome sizes. The characterizations of the three mitogenomes enrich our knowledge on the Lepidopteran mitogenome and provide us genetic information to reconstruct the phylogenetic tree. Finally, the phylogenetic results conﬁrmed the position of the genus Davidina in the subfamily Satyrini, had a closer phylogenetic relationship with Oeneis , and the phylogenetic analysis supported the formation of Oeneis buddha as an independent taxon in Oeneis .


Introduction
Satyrinae is the most diversified butterfly group in Nymphalidae (Lepidoptera). It includes approximately 400 genera and 3000 species. They are distributed in all continents except Antarctica [1,2]. Satyrinae is traditionally divided into nine tribes [3,4]. However, the higher-level classification of this family is still controversial [2,5,6]. Oeneis is a type of butterfly that adapts to cold climates and is distributed at high altitudes. Due to their similarity in morphology, the boundaries at the species level have been very blurred. The early classification of this type of butterfly was mainly based on appearance and some gene fragments [7]. Davidina is a mysterious genus. The genus was mistaken as a Pieridae member early on because of the black veins and the absence of any eye spots on the wings [8]. Kusnezov [9] moved the genus into Satyrinae and later placed it in Satyrini (equivalent to the current Satyrinae) [10]. Similarly, Paroeneis had not yet been sampled in previous molecular phylogeny studies. This genus also lacks eye spots and has limited distribution. Due to the intricacy, it is essential to use mtDNA genome sequences of these species for classification and phylogenetic analysis.
The mtDNA of insects is approximately 16 kb in length. This circular DNA encodes 13 mitochondrial proteins, 22 mitochondrial tRNAs, and two mitochondrial-specific ribosomal RNAs, specifically 12S rRNA and 16S rRNA [11]. Additionally, it contains one non-coding DNA region, the A + T-rich region (or control region), which controls replication and transcription [12]. Due to rapid evolution, cellular abundance, and the absence of introns, mitochondrial sequences can be easily amplified. In addition, they have a Diversity 2021, 13, 468 2 of 18 compact size, maternal inheritance, conservative features in the genetic organization, lack of extensive recombination, and a higher mutation rate than nuclear sequences. It is currently widely used in molecular evolution research, population genetic comparison, species identification, phylogenetic analysis, and evolutionary genomics [13][14][15][16][17][18]. In particular, phylogenetic analysis based on the mitochondrial whole genome demonstrated improved resolution of the inferred phylogenetic tree compared to phylogenetic trees based on partial gene fragments. In previous studies of mitochondrial genomes of Satyrinae, only partial genome sequences of limited taxon sampling have been performed. Therefore, it is necessary to study the phylogeny of Satyrinae according to denser taxon and mitochondrial sequence. To date, complete or near-complete mitogenomes have only been sequenced for 41 Satyrinae species, belonging to the subfamilies Satyrini, Melanitini, and Elymniini [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36].
In the present study, three new sequences of Satyrinae were obtained, seven sequences of Satyrinae were extracted from the published butterfly genome data set, and we reannotated the published mitochondrial genomes of eight satyrid. To better understand the functions of related genes, we analyzed the relative synonymous codon usage (RSCU) and AT skew values of protein coding genes (PCGs) and compared them with those of other Lepidopteran sequences. Furthermore, the phylogeny of Satyrinae and related species was constructed, and the relationships between these taxa were discussed. The divergence time of three species in Satyrinae was evaluated.

Specimens and DNA Extraction
Specimens of Hipparchia autonoe, Paroeneis palaearctica, and Oeneis buddha were collected from the Qinghai Province of China in 2017. Genomic DNA was extracted from the 95-100% ethanol-preserved muscle tissue of three satyrid species using an EasyPure Genomic DNA Kit according to the manufacturer's instructions ((TransGen Biotech, Beijing, China).

Mitogenome Sequencing, Assembly, and Annotation
The entire mitochondrial genome was sequenced using an Illumina HiSeq2500platform (Biomarker Technologies, Beijing, China). We used Geneious 11.0.2 software to retrieve other Satyrinae species as reference sequences and used the default parameters to qualitatively prune the original sequence [37]. The 13 PCGs were predicted by comparison with the homologous sequence of reference mitogenomes and finding the open reading frames (ORFs) based on the invertebrate mitochondrial genetic code. The locations of 22 tRNAs were identified by using the MITOS WebServer (http://mitos.bioinf.uni-leipzig.de/index. py, accessed on 1 July 2021) [38]. The two ribosomal RNA genes (rrnS and rrnL) and the A + T-rich region were determined by the locations of adjacent genes (trnL1 and trnV) and alignment with the homologous sequences of reference mitogenomes. A Mitochondrial Genome ring map was constructed online utilizing Organellar genome-draw [39].
In this study, we found the genome data of seven Satyrinae species (Coenonympha arcania, Lasiommata maera, Lasiommata megera, Maniola bathseba, Maniola Cecilia, Melanargia galathea, Melanargia ines) from the published butterfly genome data set and selected the mitochondrial genome data we needed from them [40]. Mitochondrial genomes were spliced and annotated using Geneious11.0.2 and MITOS Web server (http://mitos.bioinf. uni-leipzig.de/index.py, accessed on 1 July 2021) [37]. The mitochondrial genome of Lasiommata deidamia was selected as the reference sequence, and mitochondrial reads were captured from the genome-wide data of seven satyrid. Mapping the whole genome data of seven satyrid to the mitochondrial genome of L. deidamia, high-coverage, and continuous mitochondrial reads formed sequence blocks (bins), the individual bins, or connected them to Contigs according to the overlap of bins to replace the original reference. The sequence was used as the target sequence for the next mapping (baiting sequencing), the whole genome data were sequentially and repeatedly mapped to the newly generated Diversity 2021, 13, 468 3 of 18 target sequence to extend the sequence, and, finally, extended to the length of the complete mitochondrial genome [41].

Sequence Analyses
The nucleotide composition and skew, codon usage of PCGs, and relative synonymous codon usage (RSCU) values of each PCG were calculated using PhyloSuite v1.2.1 [42], and tandem repeat units of the A + T-control region were analyzed with Tandem Repeats Finder online server (http://tandem.bu.edu/trf/trf.html, access on 1 July 2021) [43].

Phylogenetic Analysis
A total of 54 mitogenomes of Nymphalidae insects were collected to analyze the phylogenetic relationships. Three newly sequenced specimens and 48 available mitogenomes were selected as ingroups (Tables S2 and S3). Three species, Polyura arja, Polyura Schreiber, and Polyura nepenthe (Nymphalidae) were employed as outgroup taxa. The nucleotide sequences of all 13 PCGs and two rRNA genes were used to elucidate the phylogenetic relationships of this tribe. All the available mitochondrial genomes were downloaded from GenBank for phylogenetic analyses.
The optimal partitioning scheme and nucleotide substitution model for Bayesian inference (BI) and maximum likelihood (ML) phylogenetic analyses with PartitionFinder2.1.1 incorporated into PhyloSuitev1.2.1 [48]. BI analyses was performed using MrBayesv3.2.2 [49]. The ML phylogenetic analysis was conducted by IQ-TREE v1.6.8 [50] using the ultrafast bootstrap (UFB) algorithm with 1000 replicates. Bootstrap support (BS) values were evaluated with 1000 replicates. BI analysis was conducted using four independent Markov chains set to run for 10 million generations with sampling every 1000 generations. The initial 25% of samples were discarded as burn-in and the remaining samples were used to generate a consensus tree and estimate the posterior probabilities (PP).

General Mitogenomic Features
H. autonoe, P. palaearctica, and O. buddha complete mitogenomes were found to be 15,435, 15,942, and 15,259 bp, respectively (Table S1). P. palaearctica Satyrinae mitogenomes were the largest of the 51 currently published Satyrinae genomes, and the other two mitogenome sizes were in the range of the 51 Satyrinae mitogenomes (from 14,675 bp of Maniola bathseba to 15,942 bp of P. palaearctica, Table S2). The genetic composition of the three species was found to be similar to that of other insect mitochondrial DNA, including 13 PCGs, two rRNA genes (16S rrnL and 12S rrnS), 22 tRNA genes, and an AT-rich region (A + T-rich region) ( Figure 1). Among the 37 genes encoded by the mitochondrial genome, 14 were determined to be encoded by the N-strand (minority-strand), including four PCGs, eight tRNA genes, and two rRNA genes; the other 23 genes are encoded by the J-strand. The A + T-rich region was found between the rrnS gene and tRNA Met gene. The arrangement of trn M-trn I-trn Q was different from that of trn I-trn Q-rn M [51], but this arrangement is typical in Lepidopterans.  The A + T contents in the three satyrid mitogenomes were 79.1%, 78.2%, and 79.6%, which were significantly higher than those of G+C and had obvious AT-bias (Table S2). On the whole, the A + T content in each gene component showed a decreasing trend as follows: A + T-rich region > rRNA > tRNA > PCGs. The AT skews of the mitogenomes were −0.017, −0.022, and −0.029, respectively, whereas the AT skew of the PCGs and the first, second, and third codons were all negative, indicating that the T content is higher than that of A ( Table 1). The AT skew values of all three mitogenomes were both within the corresponding values of other Satyrinae species, ranging from −0.055 in Neope muirheadii to −0.017 in H. autonoe (Table S2).

Protein Coding Genes and Codon Usage
The lengths of the 13 PCGs of H. autonoe, P. palaearctica, and O. buddha were 11,205, 11,208, and 11,202 bp, respectively, accounting for 70.3%, 73.5%, and 73.4% of the total sequence. Except for COX1, all PCGs in the H. autonoe mitogenome were found to be initiated by typical ATN codons (COX2, ATP6, COX3, ND4, ND4L, and Cyt b with ATG; ND2, ATP8, and ND6 with ATT; ND3 and ND5 with ATC; ND1 with ATA). The start codons of 10 PCGs (ND1, ND2, ND4, ND4L, COX1, COX2, COX3, ATP8, ATP6, and Cyt b) in the three species were consistent with each other. Unlike that in H. autonoe, in P. palaearctica and O. buddha, ND3 and ND5 were determined to start with an ATT codon. Results showed that COX1 in the three Satyrinae species starts with the alternative starting codon CGA, which is similar to that in all other known Lepidopteran mitogenomes. The termination codons of the three satyrids in this paper also had some similarities; except for the COX2, ND4, and ND5 genes that terminate with a single T, the other PCGs terminate with a complete TAA codon.

The rRNAs and tRNAs
Like other Lepidopteran mitogenomes, the mitogenomes of these two species were found to have two ribosomal RNA genes (rrnL and rrnS) between tRNA Leu (CUN) and tRNA Val and between tRNA Val and the A + T rich region, respectively. The sizes of rrnL gene of H. autonoe, P. palaearctica, and O. buddha mitogenomes were 1358, 1362, and 1362 bp, respectively. The A + T contents of the rrnL genes in H. autonoe, P. palaearctica, and O. buddha were 83.4%, 83%, and 83.6%. The rrnS gene was found to be 498 bp long in H. autonoe, 852 bp in P. palaearctica, and 768 bp in O. buddha, and their A + T contents were 83.7%, 78%, and 86.2%, respectively (Table 1). These results were similar to those of other Satyrinae mitogenomes (Table S2).
H. autonoe, p. palaearctica, and O. buddha mitogenomes were all found to contain 22 tRNA genes. The length of the 22 tRNA genes ranged from 60 bp to 71 bp. All tRNA genes except tRNA Ser (AGN) were determined to have a typical clover structure, whereas tRNA Ser (AGN) lacked the dihydrouridine (DHU) stem arm and formed a simple ring instead.

The rRNAs and tRNAs
Like other Lepidopteran mitogenomes, the mitogenomes of these two species were found to have two ribosomal RNA genes (rrnL and rrnS) between tRNA Leu (CUN) and tRNA Val and between tRNA Val and the A + T rich region, respectively. The sizes of rrnL gene of H. autonoe, P. palaearctica, and O. buddha mitogenomes were 1358, 1362, and 1362 bp, respectively. The A + T contents of the rrnL genes in H. autonoe, P. palaearctica, and O. buddha were 83.4%, 83%, and 83.6%. The rrnS gene was found to be 498 bp long in H. autonoe, 852 bp in P. palaearctica, and 768 bp in O. buddha, and their A + T contents were 83.7%, 78%, and 86.2%, respectively (Table 1). These results were similar to those of other Satyrinae mitogenomes (Table S2).
H. autonoe, p. palaearctica, and O. buddha mitogenomes were all found to contain 22 tRNA genes. The length of the 22 tRNA genes ranged from 60 bp to 71 bp. All tRNA genes except tRNA Ser (AGN) were determined to have a typical clover structure, whereas tRNA Ser (AGN) lacked the dihydrouridine (DHU) stem arm and formed a simple ring instead.

Intergenic and Overlapping Spacer Regions
Although the mitogenomes of animals comprise a tightly arranged structure without introns, there are still multiple intergenic spacer sequences except for in the rich region.

Intergenic and Overlapping Spacer Regions
Although the mitogenomes of animals comprise a tightly arranged structure without introns, there are still multiple intergenic spacer sequences except for in the rich region. The gene overlapping sequences of the three mitogenomes had a total length of 33 bp at 10 locations in H. autonoe, 36 bp at 11 locations in P. palaearctica, and 36 bp at 11 locations in O. buddha. One was found to be composed of "AAGCCTTA" at the tRNA Trp and tRNA Cys junction ( Figure 4A); the second was a shorter sequence of "TCTAA" located at the COX1 and tRNA Leu (CUN) junction ( Figure 4C). Like previously described Lepidopteran mitogenomes, the overlapping sequences (a 7-bp overlap) in the three species were found to be present between ATP8 and ATP6 genes ( Figure 4B).

A + t-Rich Region
The lengths of the A + t-rich region of the three Satyrinae mitochondrial genomes were as follows: H. autonoe = 896 bp, P. palaearctica = 1039 bp, and O. buddha = 452 bp, located between srRNA and tRNAIle ( Figure 1 and Table 1). The A + t content was 94.6% in H. autonoe, 61.8% in P. palaearctica, and 91.4% in O. buddha. The A + t content of P. palaearctica was the smallest of the 49 Satyrinae A + T contents published to date. locations in O. buddha. One was found to be composed of "AAGCCTTA" at the tRNA Trp and tRNA Cys junction ( Figure 4A); the second was a shorter sequence of "TCTAA" located at the COX1 and tRNA Leu (CUN) junction ( Figure 4C). Like previously described Lepidopteran mitogenomes, the overlapping sequences (a 7-bp overlap) in the three species were found to be present between ATP8 and ATP6 genes ( Figure 4B).

Phylogenetic Analyses
This study was based on the aforementioned three newly acquired species of Satyrinae, combined with the known complete mitochondrial genome sequences of 48 other Satyrinae species, whereas the mitogenomes of Nymphalidae (P. arja, P. schreiber and P. nepenthe) were used as outgroups (Figures 5 and 6). The connected data set was analyzed by Bayesian inference (BI) and maximum likelihood (ML). Satyrinae is one of the most abundant taxa of butterfly insects and is an important model material for many research fields such as ecology, functional genomics, and bionics. However, there are many problems that need to be resolved in research on the system classification of Satyrinae low-level and high-level elements [2][3][4]52]. To further understand their evolutionary relationship, we conducted a preliminary investigation based on mitochondrial data. verified and were supported in our analyses [2,4]. However, it is important to note that only eight of the 13 subtribes have been analyzed. To better deal with the phylogenetic relationships within Satyrinae, more mitochondrial genome data need to be added for analysis.

Divergence Time Estimation of Satyrinae Species
Times of divergence within Satyrinae were estimated using a Bayesian approach in r8s (Figure 7). The divergence time of Satyrinae was approximately 47 MYA. Among them, Melanitini divergence time was approximately 0.5057MYA, whereas Satyrini divergence time was approximately 43.2742 MYA.
The divergence time of Satyrinae was in the Paleogene Period. This is related to the rapid evolution of angiosperms in terrestrial ecosystems during this period. According to the fossil record, the tertiary period was the prime of species evolution and angiosperms began to dominate [56][57][58][59]. At the same time, one of the fastest and most intense global warming events occurred during this period. During this time, high temperatures and a warm ocean created a moist and mild earth environment; except for the deserts, the In our phylogenetic relationship, Lethina is closer to Parargina and Mycalesina, and Satyrina is closer to Melanargia and Maniolina. This kind of relationship is similar to Marín's research results [6]. In this study, we selected 18 species of Lethina. In the obtained phylogenetic relationship, Neope and Ninguta are sister groups of Lethe, and similar relationships have been established in previous studies [3,28].
Davidina is a mysterious genus. The genus was mistaken as a Pieridae member early on because of the black veins and the absence of any eye spots on the wings [8]. Kusnezov [9] moved the genus into Satyrinae and later placed it in Satyrini (equivalent to the current Satyrinae) [10]. Additionally, the latest research proved Davidina is not a local monotypic Chinese endemic genus, as has been previously supposed, but is composed of nine species and has a broad distribution area in the Holarctic region including Europe and America [53]. Similarly, Paroeneis had not yet been sampled in previous molecular phylogeny studies. This genus also lacks eye spots and has limited distribution. In this study, we found that it is very close to Davidina.
Due to the wide-ranging variations, the classification of the genus Oeneis is difficult. There are still many unresolved problems: (1) Wing patterns are quite similar among different species, whereas (2) great variations in wing patterns are noted even within the same species, for example, in Oeneis urda. In addition, in the morphological features, Davidina and Oeneis share common morphological characteristics. Their first thoracic limbs are significantly degenerated and traced, male and female genitalia are similar to those of the genus Oeneis, and the wing venation is also similar to that of the genus Oeneis [9]. In our phylogenetic relationships, O. urda and Davidina have a closer phylogenetic relationship. A similar relationship was established in previous research [53,54]. In the genus Oeneis, the morphological characteristics of male genitalia in O. buddha were classified as an independent group due to the presence of features not found in other Oeneis groups (valve without serrations) [55]. The present phylogeny supports the morphological classification of O. buddha forming an independent group.
Satyrini represents the most diverse tribe in Satyrinae (approximately 2200 species) and is divided into 13 subtribes [6]. In this study, three new mitochondrial genomes of Satyrini were determined, increasing the number of published Satyrini mitogenomes to 51. In previous studies, the close relationships among Parargina, Lethina, Mycalesina, Melanargiina, Satyrina, Coenymphina, Maniolina, and Ypthimina groups have been verified and were supported in our analyses [2,4]. However, it is important to note that only eight of the 13 subtribes have been analyzed. To better deal with the phylogenetic relationships within Satyrinae, more mitochondrial genome data need to be added for analysis.

Divergence Time Estimation of Satyrinae Species
Times of divergence within Satyrinae were estimated using a Bayesian approach in r8s (Figure 7). The divergence time of Satyrinae was approximately 47 MYA. Among them, Melanitini divergence time was approximately 0.5057MYA, whereas Satyrini divergence time was approximately 43.2742 MYA.
The divergence time of Satyrinae was in the Paleogene Period. This is related to the rapid evolution of angiosperms in terrestrial ecosystems during this period. According to the fossil record, the tertiary period was the prime of species evolution and angiosperms began to dominate [56][57][58][59]. At the same time, one of the fastest and most intense global warming events occurred during this period. During this time, high temperatures and a warm ocean created a moist and mild earth environment; except for the deserts, the surface was completely covered by forests and there was a great abundance of plants to provide food for a wide range of plant-feeding insects.
In this study, the divergence time of H. autonoe, P. palaearctica, and O. buddha was approximately 27.8749 MYA (in the Paleogene Oligocene). The three Satyrinae samples were from Qinghai Province in the northeast of the Qinghai-Tibet Plateau, and P. palaearctica and O. buddha are alpine species with special living environments. According to the history of the evolution of the global environment, during the Oligocene (33.7 MYA-23.5 MYA), the temperature was generally low and the Qinghai-Tibet Plateau began to lift up. The Paleozoic era had prosperous angiosperms, providing prerequisites for the evolution of alpine species of butterflies towards high-altitude environments adapted to low temperature, hypoxia, and climate change, as well as rapid adaptation to radiation events that occurred in the Quaternary period. In this study, the divergence time of three satyrids was predicted to coincide with the geological and environmental events mentioned previously herein. Diversity 2021, 13, x FOR PEER REVIEW 16 of 19 Figure 7. Divergence times of Satyrinae: maximum clade credibility tree with median age and 95% confidence interval estimated using a Bayesian uncorrelated relaxed clock implemented in r8s.

Conclusions
In this study, we identified three new mitochondrial whole genomes, including H. autonoe, P. palaearctica, and O. buddha. At the same time, we also re-spliced and annotated However, combined with previous research results, it is shown that the divergence time of Satyrinae predicted in this paper is similar to Espeland et al. [60] and Pena et al. [5] and slightly later than other research results. For example, Wahlberg et [4].
The reasons for the differences between the research results of different scholars and this study may be as follows.
(1) When the phylogenetic tree is reconstructed, the choice of sequence data will inevitably have a certain impact on the research results, because of the difference in the evolution rate of different sequences. This study used mitochondrial genome protein-coding genes for phylogenetic tree reconstruction and molecular clock analysis, while previous studies all used partial mitochondrial genome sequence fragments or mitochondrial partial genes combined with one or two nuclear genes to reconstruct the phylogenetic tree. The molecular clock analysis results are slightly different. (2) The results of the study are affected by the different embedding locations of the fossil correction points and the different setting of the time during molecular clock calculation. (3) The value given by the molecular clock calculation is the average value of a confidence interval, that is, the divergence time of each clade in the research results is not a given value but a time range, and the data results given only provide convenience for intuitive comparison.

Conclusions
In this study, we identified three new mitochondrial whole genomes, including H. autonoe, P. palaearctica, and O. buddha. At the same time, we also re-spliced and annotated the mitochondrial genomes of seven Satyrinae species, which are very similar to the sequence structures of other Satyrinae species. Their mitochondrial genomes are highly conserved in base content and composition, genome size and sequence, protein coding genes and codon usage, and tRNA secondary structure. Based on 13 PCGs and two rRNA phylogenetic analyses of BI and ML, a well-resolved topological structure was obtained, with high support for each branch. The divergence time of Satyrinae is the Paleogene Period. The phylogenetic results confirmed the position of the genus Davidina in the subfamily Satyrini and had a closer phylogenetic relationship with Oeneis. The phylogenetic analysis supported the formation of Oeneis buddha as an independent taxon in Oeneis.
In this study, the number of mitochondrial genomes in Satyrini was increased to 51. Satyrinae is a large subfamily of Nymphalidae with abundant species' resources, but there are few species for which completed mitochondrial genome sequencing has been performed. Therefore, it is urgent to obtain more mitochondrial genome sequences of Satyrinae insects to better solve the phylogenetic relationship of Satyrinae.