Three Complete Mitochondrial Genomes of Orestes guangxiensis, Peruphasma schultei, and Phryganistria guangxiensis (Insecta: Phasmatodea) and Their Phylogeny

Simple Summary Twenty-seven complete mitochondrial genomes of Phasmatodea have been published in the NCBI. To shed light on the intra-ordinal and inter-ordinal relationships among Phasmatodea, more mitochondrial genomes of stick insects are used to explore mitogenome structures and clarify the disputes regarding the phylogenetic relationships among Phasmatodea. We sequence and annotate the first acquired complete mitochondrial genome from the family Pseudophasmatidae (Peruphasma schultei), the first reported mitochondrial genome from the genus Phryganistria of Phasmatidae (P. guangxiensis), and the complete mitochondrial genome of Orestes guangxiensis belonging to the family Heteropterygidae. We analyze the gene composition and the structure of the three mitochondrial genomes. We recover the monophyly of Phasmatodea and show the sister-group relationship between Phasmatodea and Mantophasmatodea after removal of the Embioptera and Zoraptera species. We recover the monophyly of Heteropterygidae and the paraphyly of Diapheromeridae, Phasmatidae, Lonchodidae, Lonchodinae, and Clitumninae. Abstract Insects of the order Phasmatodea are mainly distributed in the tropics and subtropics and are best known for their remarkable camouflage as plants. In this study, we sequenced three complete mitochondrial genomes from three different families: Orestes guangxiensis, Peruphasma schultei, and Phryganistria guangxiensis. The lengths of the three mitochondrial genomes were 15,896 bp, 16,869 bp, and 17,005 bp, respectively, and the gene composition and structure of the three stick insects were identical to those of the most recent common ancestor of insects. The phylogenetic relationships among stick insects have been chaotic for a long time. In order to discuss the intra- and inter-ordinal relationship of Phasmatodea, we used the 13 protein-coding genes (PCGs) of 85 species for maximum likelihood (ML) and Bayesian inference (BI) analyses. Results showed that the internal topological structure of Phasmatodea had a few differences in both ML and BI trees and long-branch attraction (LBA) appeared between Embioptera and Zoraptera, which led to a non-monophyletic Phasmatodea. Consequently, after removal of the Embioptera and Zoraptera species, we re-performed ML and BI analyses with the remaining 81 species, which showed identical topology except for the position of Tectarchus ovobessus (Phasmatodea). We recovered the monophyly of Phasmatodea and the sister-group relationship between Phasmatodea and Mantophasmatodea. Our analyses also recovered the monophyly of Heteropterygidae and the paraphyly of Diapheromeridae, Phasmatidae, Lonchodidae, Lonchodinae, and Clitumninae. In this study, Peruphasma schultei (Pseudophasmatidae), Phraortes sp. YW-2014 (Lonchodidae), and species of Diapheromeridae clustered into the clade of Phasmatidae. Within Heteropterygidae, O. guangxiensis was the sister clade to O. mouhotii belonging to Dataminae, and the relationship of (Heteropteryginae + (Dataminae + Obriminae)) was recovered.


Introduction
Stick and leaf insects (Phasmatodea) belong to an order of polyneopteran insects, which includes over 3000 recognized species subdivided into approximately 500 genera, distributed across major landmasses [1]. Phasmatodea are the longest insects among extant species, mainly distributed in tropical and subtropical regions, with a few species in temperate regions [2]. Phasmatodea are herbivorous insects and have a strong ability to masquerade as bark, leaves, and twigs, which provide them with camouflage [3][4][5].
Consequently, it is difficult to explore the phylogenetic relationship of Phasmatodea by the degree of convergence of morphology [6].
The high-level phylogenetic relationships of Phasmatodea are currently not sufficiently clear [23,24]. Several authors have suggested that Phasmatodea should be divided into two suborders: Timematodea and Euphasmatodea (Verophasmatodea) [1,5,25,26]. The Euphasmatodea includes thirteen families: Aschiphasmatidae, Damasippoididae, Prisopodidae, Anisacanthidae, Bacillidae, Heteropterygidae, Phylliidae, Agathemeridae, Heteronemiidae, Pseudophasmatidae, Diapheromeridae, Lonchodidae, and Phasmatidae [2]. However, the wingless Nearctic walking-stick, Timema of family Timematidae, is the only genus within the one family in Timematodea [27]. Neverthless, Simon et al. demonstrated a basal dichotomy of Aschiphasmatodea and the Neophasmatodea in Euphasmatodea [28]. This result was also supported by Tihelka et al. [29]. The monophyly of Phasmatodea remains in dispute because of the phylogenetic position of genus Timema. Most studies based on morphology, transcriptome, and mitochondrial genome data have shown that the monophyly of Phasmatodea can be confirmed, and Timema is recognized as the sister group to all remaining phasmids [28,[30][31][32][33]. Nevertheless, data on mitochondrial genomes considered that Timema did not belong to the Euphasmatodea, but grouped with other orders (e.g., Orthoptera and Embioptera) [6,34]. This result coincided with the study about the morphology of Timema species in egg that indicated that Timema was a separate lineage [35]. Based on 16S, 12S RNA genes, and protein codon genes, Song et al. [22] constructed nine phylogenetic trees, and eight of these failed to recover the monophyly of Phasmatodea because the clade of Embioptera and Zoraptera that existed in long-branch attraction (LBA) was a sister group to Euphasmatodea, whereas only one phylogenetic tree supported the monophyly of Phasmatodea, as found in Song et al. [36]. In many phylogenetic analyses, especially those based on mitochondrial genomes, biases associated with LBA have been found [37].
Insect mitochondrial genomes normally have 37 genes (thirteen protein-coding genes, two ribosomal RNAs, and 22 transfer RNA genes) and a control region (CR) and are usually a 14-20 kb double-stranded circular molecule [38]. The mitochondrial genome has been widely used for phylogenetic analyses due to its simple and stable gene organization, lack of genetic recombination, and fast evolution rate, etc. [21,[39][40][41]. Many researchers have discussed the phylogenetic relationships among insect orders using mitochondrial genomes, such as Diptera [42], Orthoptera [43], Hymenoptera [44], and Coleoptera [45]. This included the first acquired complete mitochondrial genome from the Pseudophasmatidae and the first reported mitochondrial genome of Phryganistria (Phasmatidae). Moreover, we analyzed the gene composition and the structure of the three mitochondrial genomes.

Sampling Collection and DNA Extraction
Orestes guangxiensis and Ph. guangxiensis were collected from Jinxiu, Guangxi province, China, whereas Pe. schultei was retrieved from an insect pet market in China, the source area being Northern Peru. According to their morphological characters, these specimens were identified by JY Zhang and stored at −40°C in the Zhang laboratory, College of Life Sciences and Chemistry, Zhejiang Normal University, China. Total DNA was extracted from a piece of foreleg muscle using a Universal Genomic DNA Kit (Co Win Biosciences Company, Beijing, China).

PCR Amplification and Sequencing
The DNA from each of the three species was amplified using eight pairs of universal primers, as described in Zhang et al. [46], but there were still some vacancies. We then used Primer Premier 5.0 to design species-specific primers based on known sequences from universal primers [47] (Table S1) and used normal PCR (product length <3000 bp) as well as long PCR (product length >3000 bp) methods for amplification [46]. Takara Taq polymerase and Takara LATaq DNA polymerase were used, respectively (Takara, Dalian, China), in a 50 µL reaction volume. Reaction systems and cycling conditions for normal PCR and long PCR were as described in Zhang et al. [46]. All PCR products were sequenced in both directions using the primer-walking method and ABI3730XL by Sangon Biotech Company (Shanghai, China).

Mitochondrial Genome Annotation and Sequence Analyses
The fragments obtained by Sanger dideoxy sequencing were assembled with DNAS-TAR Package v.7.1 [48]. The tRNA genes were identified using the MITOS web server (http://mitos.bioinf.uni-leipzig.de/index.py (accessed on 15 July 2021)) [49]. Based on the homologous sequences of mitochondrial genomes from other stick insect species, we used Clustal X [50] to determine the two rRNA genes (12S and 16S rRNA). The remaining 13 protein-coding genes were analyzed using Mega 7.0 [51] to translate amino acids using the invertebrate mitochondrial genetic code and find open reading frames [52]. The AT content, codon usage, and relative synonymous codon usage (RSCU) of protein-coding genes were calculated by PhyloSuite 1.2.2 [53]. GC and AT skews were calculated according to the formula: AT skew = (A − T)/(A + T), GC skew = (G − C)/(G + C) [54].

Phylogenetic Analyses
To illuminate the phylogenetic relationships of Phasmatodea, we first performed maximum likelihood (ML) and Bayesian inference (BI) analyses based on data from 85 species, including the three newly determined sequences, sixty-five previously sequenced mitochondrial genomes, and 13 PCGs of seventeen species of Phasmatodea assembled from transcriptome data [14,33,34,36, (Table 1, Table 2 and Table S2). Long-branch attraction between Embioptera and Zoraptera appeared in both ML and BI phylogenetic trees and affected the stability of the topology. This phenomenon probably occurred because of phylogenetic artifacts generated by the mitochondrial genome [78]. For LBA, the following measures are recommended: increasing samples [79]; removing long-branch groups [80]; or removing sites with rapid evolutionary rates [81]. Adding sequences was not feasible because all current mitochondrial genome sequences of Embioptera and Zoraptera were used in this study. Then, we attempted to remove sites with fast evolutionary rates, but nevertheless, the results were not very reliable. Therefore, we reconstructed BI and ML phylogenetic trees with 13 concatenated PCG sequences of 81 species after removing representative species of Embioptera and Zoraptera. Three species of Archaeognatha, Nesomachilis australica, Pedetontus silvestrii, and Trigoniophthalmus alternatus were used as outgroups ( Table 2). We aligned each of the 13 protein-coding genes using Clustal W in the program Mega 7.0 [51] and used the program Gblock 0.91b to identify conserved regions [82]. The resulting alignments were concatenated with Geneious 8.1.6 [83]. Because the third codon positions had saturated by DAMBE 4.2.13 [84], we used the Bayesian inference (BI) and maximum likelihood (ML) methods with the first and second codon datasets to analyze the phylogenetic relationships. ML analysis was implemented by IQ-TREE v.2.1.2 with the best model GTR + I + G that was acquired by ModelFinder [85,86]. BI analysis was carried out by MrBayes 3.2 with GTR + I + G [87] and was set for 10 million generations with sampling every 1000 generations, and the first 25% of generations were discarded as burn-in.

Mitochondrial Genome Organization and Composition
The lengths of the three complete mitochondrial genomes of O. guangxiensis, Pe. schultei, and Ph. guangxiensis were 16,869 bp, 15,896 bp, and 17,005 bp, respectively ( Figure 1). All three genomes were deposited in GenBank, with accession numbers MW450873, MW450874, and MW450875, respectively. Mitochondrial genomes of the three species had the same genes and gene order as those of other published stick insects, which have 37 genes, including 13 PCGs, 22 tRNA genes, and two rRNA genes. Currently, the gene arrangement of published stick insects is similar to the assumed common ancestor of insects, except for Ramulus hainanense (CR-trnM-trnQ-trnI-CR-trnI-trnQ-trnM) and Megalophasma granulatum (trnR-trnA) [6,14,15,22,30,[89][90][91]. According to the previously published complete mitochondrial genomes of stick and leaf insects, we found that the differing lengths of the Phasmatodea genomes (15,590-18,248 bp) were caused mainly by the size of the A + T-rich region, gene overlaps, and different intergenic nucleotides (IGNs). The sequence length of Pe. schultei (15,896 bp), with a short A + T-rich region (<1500 bp), was the shortest after that of R. hainanense (15,590 bp). The three species have short intergenic regions Insects 2021, 12, 779 6 of 18 ranging from 1 to 18 bp. The whole mitochondrial genome of Ph. guangxiensis, which contained additional IGNs (136 bp), was longer than that of O. guangxiensis (Tables S3-S5). The nucleotide composition of the O. guangxiensis, Pe. schultei, and Ph. guangxiensis mitochondrial genomes had a high A + T bias of 75.6%, 76.6%, and 76.8%, respectively. All three species showed a positive AT skew and negative GC skew ( Table 3). The content of A was more than T, and the content of C was higher than G, which also occurred in the sequences of previously studied stick insects (Table S6) [6,14,15,22,30,[89][90][91].

Protein-Coding Genes and Codon Usages
The locations of the 13 PCGs within the three mitochondrial genomes were identical to those of most stick insects (Tables S3-S5). Four PCGs (ND1, ND4, ND4L, and ND5) were located in the minority strand (N-strand), whereas the remaining PCGs were coded on the majority strand (J-strand). The total lengths of the 13 protein-coding genes (PCGs) in O. guangxiensis, Pe. schultei, and Ph. guangxiensis were 11,112 bp, 11,100 bp, and 11,121 bp, respectively (Table 3). Among the three mitochondrial genomes, all PCG initiation codons used ATN (N represents A, G, C, or T), except for ND4L of Pe. schultei, which started with GTG. GTG as a start codon has also been found in ND4 of Megalophasma granulatum [91]. ATN is an accepted canonical initiation codon for insect mitochondrial genomes [92]. Of the stick insects that used ATN as an initiation codon, most used ATA, ATG, and ATT, with only a few using ATC [6,14,15,22,30,[89][90][91]; only ATP8 (Pe. schultei) used the ATC start codon in this study. The typical termination codons (TAA/TAG) were found in most PCGs, except for some incomplete terminal codons, such as T used for COX2 (O. guangxiensis, Pe. schultei, and Ph. guangxiensis), ND1 (Ph. guangxiensis), ND3 (O. guangxiensis and Pe. schultei), ND4L (Pe. schultei), and ND5 (O. guangxiensis and Pe. schultei). Incomplete termination codons have also been found in other insects [65,[93][94][95][96]. Incomplete stop codons play a significant role in polycistronic transcription cleavage and polyadenylation processes [97]. High A + T bias was also found in the PCGs of O. guangxiensis, Pe. schultei, and Ph. guangxiensis, which were 74.2%, 75.5%, and 75.8%, respectively. The PCGs of the majority strand displayed positive AT skews and negative GC skews, whereas the minority strand displayed negative AT skews and positive GC skews. The A skew (the content of A > T) and C skew (the content of C > G) of the minority strand was greater than on the majority strand (Table 3).
We calculated the relative synonymous codon usage (RSCU) of the three mitochondrial genomes ( Figure 2, Table S7). The results showed that A or T nucleotides were used in high frequency in the third codon position compared to other nucleotides, and that A was used more often than T. The most frequent codons used were UUU (Phe), UUA (Leu), AUU (Ile), and AUA (Met) and were used >280 times in the PCGs of O. guangxiensis, Pe. schultei, and Ph. guangxiensis mitochondrial genomes. In contrast, codons with a third codon G or C were used very rarely (≤10), such as CUC (Leu), UGC (Cys), CGC (Arg), GCG (Ala) (≤5), etc. This may be a kind of AT mutation bias that has an obvious influence on codon usage [98,99].

Ribosomal RNAs and Transfer RNAs
The   The mitochondrial genomes of these three stick insects, similar to other species in Phasmatodea, contained two rRNAs genes [6,15]. The 16S rRNA gene in O. guangxiensis, Pe. schultei, and Ph. guangxiensis was located between trnL1 and trnV, with a length of 1291 bp, 1281 bp, and 1328 bp, respectively, whereas the 12S rRNA was located between trnV and the CR, with sizes of 788 bp, 774 bp, and 796 bp, respectively. The AT content of the two rRNAs in O. guangxiensis (78.3%), Pe. schultei (77.7%), and Ph. guangxiensis (77.9%) were each higher than the average AT content of the 13 PCGs (Table 3). We found that the AT skew values of the two rRNAs in O. guangxiensis, Pe. schultei, and Ph. guangxiensis were 0.25, 0.22, and 0.20, respectively. Meanwhile, the GC skew was highly negative, with values around 0.3 (Table 3).

A + T-Rich Region
The large non-coding region of O. guangxiensis, Pe. schultei, and Ph. guangxiensis between trnI and 12S rRNA was an A + T-rich region with lengths of 2238 bp, 1294 bp, and 2286 bp, respectively (Table 3). Compared with mitogenomes from other Phasmatodea, the length of the A + T-rich region in Pe. schultei was the shortest, except for Ramulus hainanense (774 bp) (FJ156750). According to the published complete Phasmatodea mitochondrial genomes, the longest A + T-rich region was found in Cryptophyllium tibetense (3701 bp) (KX091862). In the mitochondrial genomes of O. guangxiensis, Pe. schultei, and Ph. guangxiensis, the content of A + T in the control regions was 79.6%, 82.5%, and 79.1%, respectively, which was higher than other partitions of mitochondrial genomes. The A + T-rich regions of the three species each showed positive AT skew values and negative GC skew values ( Table 3). The A + T region embodied the origin sites and essential regulatory elements needed for transcription and replication [105][106][107].
Repeat regions were observed in O. guangxiensis, Pe. schultei, and Ph. guangxiensis (Figure 3). The A + T-rich region of Ph. guangxiensis possessed seven copies of tandem repeats regions with a length of 106 bp, whereas the A + T-rich region of Pe. schultei contained six tandem repeat regions of a 32 bp sequence. However, two repeats (172 bp) in O. guangxiensis were not tandem ( Figure 3). Tandem repeats in the A + T rich region have also been observed in many Phasmatodea species. In the research of Kômoto et al., the presence of tandem repeats in the A + T region was also detected in eight Phasmatodea species, such as Entoria nuda (ten tandem repeats of a 128-129 bp sequence) and Ramulus mikado (two tandem repeats of a 149 bp sequence and seven tandem repeats of 125 bp) [6]. Two identical copies of a 64 bp tandem repeat were discovered in Ramulus hainanense (FJ15676), and Eurycantha calcarata included twenty-two tandem repeats of a 32 bp fragment (MW915467). Twenty-two tandem repeats of a 21 bp sequence were found in Cryptophyllium tibetense (KX091862), and Extatosoma tiaratum possessed three tandem repeats of a 128-129 bp fragment [30]. Phraortes sp. 1 NS-2020 contained two tandem repeats with lengths of 20 bp [22]. The secondary structure of these repeat units was predicted by RNAalifold in 15 species [108][109][110]. We found that most of these tandem repeats in the A + T region could form the stem-loop structure ( Figure S4). Repeat regions of the mitochondrial A + T region also existed in other insects. The A + T-rich region of Theopompa sp.YN (Mantodea: Mantidae) contained three tandem repeats of a 200 bp sequence [96]. Two tandem repeats of a 90 bp sequence, three tandem repeats of a 100 bp sequence, and six tandem repeats of a 50 bp sequence were observed in Serratella zapekinae (Ephemeroptera: Ephemerellidae) [63]. The cause of tandem repeats may be slipped-strand mispairing in mitochondrial genome replication [111].

Intergenic and Overlap Regions
The mitochondrial genomes of Phasmatodea, including the three stick insects in this study, were compact, with intergenic regions usually not exceeding 20 bp [6,14,15,22,30,[89][90][91]. The mitochondrial genomes of O. guangxiensis, Pe. schultei, and Ph. guangxiensis contained 7, 5, and 12 intergenic regions with total lengths of 17 bp, 41 bp, and 40 bp, respectively. We observed the longest intergenic spacer between ND1 and trnL1 and the second-longest between trnI and trnQ within the mitochondrial genomes of Pe. schultei, with lengths of 18 bp and 13 bp, respectively. Overall, insertions between genes ranged from 1 to 8 residues in the three stick insects (Tables S3-S5).

Intergenic and Overlap Regions
The mitochondrial genomes of Phasmatodea, including the three stick insects in this study, were compact, with intergenic regions usually not exceeding 20 bp [6,14,15,22,30,[89][90][91]. The mitochondrial genomes of O. guangxiensis, Pe. schultei, and Ph. guangxiensis contained 7, 5, and 12 intergenic regions with total lengths of 17 bp, 41 bp, and 40 bp, respectively. We observed the longest intergenic spacer between ND1 and trnL1 and the second-longest between trnI and trnQ within the mitochondrial genomes of Pe. schultei, with lengths of 18 bp and 13 bp, respectively. Overall, insertions between genes ranged from 1 to 8 residues in the three stick insects (Tables S3-S5).

Phylogenetic Analyses
When we analyzed the phylogenetic relationship using the 13PCGs of 85 species, including Embioptera and Zoraptera, the Bayesian tree was different from the maximum likelihood inference tree, mainly in the internal topological structure of Phasmatodea (Figures S5 and S6). We found that Phasmatodea was paraphyletic because the clade of Zoraptera and Embioptera clustered into the Phasmatodea, as also reported by Song et al. based on mitochondrial genome sequence data [22,36]. Zoraptera was the sister clade to Embioptera, caused by long-branch attraction, as found in Ma et al. [20].

Phylogenetic Analyses
When we analyzed the phylogenetic relationship using the 13PCGs of 85 species, including Embioptera and Zoraptera, the Bayesian tree was different from the maximum likelihood inference tree, mainly in the internal topological structure of Phasmatodea ( Figures S5 and S6). We found that Phasmatodea was paraphyletic because the clade of Zoraptera and Embioptera clustered into the Phasmatodea, as also reported by Song et al. based on mitochondrial genome sequence data [22,36]. Zoraptera was the sister clade to Embioptera, caused by long-branch attraction, as found in Ma et al. [20].
After removal of the Embioptera and Zoraptera species, we re-performed ML and BI analyses with the remaining 81 species, which showed identical topology except for the position of Tectarchus ovobessus (Phasmatodea). Figure 4 shows that Odonata was the basal group of Pterygota, and Ephemeroptera was a sister clade to the Polyneoptera, as also reported by some other molecular and morphological studies [62,112,113]. The monophyly of Polyneoptera also was supported. We recovered the monophyly of Phasmatodea, and the sister-group relationship between Phasmatodea and Mantophasmatodea was supported by current phylogenetic analyses after removing Zoraptera and Embioptera. Phasmatodea was divided into two branches: Timematodea and Euphasmatodea (Figure 4). group of Pterygota, and Ephemeroptera was a sister clade to the Polyneoptera, as also reported by some other molecular and morphological studies [62,112,113]. The monophyly of Polyneoptera also was supported. We recovered the monophyly of Phasmatodea, and the sister-group relationship between Phasmatodea and Mantophasmatodea was supported by current phylogenetic analyses after removing Zoraptera and Embioptera. Phasmatodea was divided into two branches: Timematodea and Euphasmatodea ( Figure 4). At the family level, our data supported the monophyly of Heteropterygidae but Diapheromeridae, Phasmatidae, and Lonchodidae were not recovered. Lonchodidae consists of two subfamilies, Lonchodinae and Necrosciinae, but did not form a clade, which was also reported by Forni et al., Kômoto et al.,and Song et al. [6,22,33]. However, the phylogenetic relationship among Dataminae, Heteropteryginae, and Obriminae of Heteropteridae is still controversial. In this study, O. guangxiensis was the sister clade to O. mouhotii belonging to Dataminae, and Heteropteryginae was the sister clade to (Dataminae + Obriminae). Meanwhile, the same conclusion was also obtained by some research that utilized morphological data [35,114] and molecular data [115]. By contrast, Bank et al., based on three nuclear data sets (18S, 28S and H3) and four mitochondrial data genes At the family level, our data supported the monophyly of Heteropterygidae but Diapheromeridae, Phasmatidae, and Lonchodidae were not recovered. Lonchodidae consists of two subfamilies, Lonchodinae and Necrosciinae, but did not form a clade, which was also reported by Forni et al., Kômoto et al.,and Song et al. [6,22,33]. However, the phylogenetic relationship among Dataminae, Heteropteryginae, and Obriminae of Heteropteridae is still controversial. In this study, O. guangxiensis was the sister clade to O. mouhotii belonging to Dataminae, and Heteropteryginae was the sister clade to (Dataminae + Obriminae). Meanwhile, the same conclusion was also obtained by some research that utilized morphological data [35,114] and molecular data [115]. By contrast, Bank et al., based on three nuclear data sets (18S, 28S and H3) and four mitochondrial data genes (COX1, COX2, 12S, and 16S), were in favor of the clade of Dataminae + (Heteropteryginae + Obriminae) [116], consistent with previous results [117][118][119]. By contrast, other studies hypothesized that Obriminae and the clade of (Heteropteryginae + Dataminae) had a close phylogenetic relationship [5,28,120,121]. The Pseudophasmatidae (Pe. schultei) should be a separate clade, but our molecular phylogenetic trees showed that it was classified into Phasmatidae, and, at the same time, Phraortes sp. YW-2014 (Lonchodidae) and species of Diapheromeridae clustered into the clade of Phasmatidae. However, some studies support a sister clade of Agathemera and Pseudophasmatidae [5,28].
At the subfamily level, both Lonchodinae and Clitumninae were recovered as a polyphyletic group. In our analysis, we showed that Lonchodinae was not monophyletic because Phraortes sp. YW-2014 (Lonchodinae) formed a clade with Medauroidea extradentata (Clitumninae) instead of the main clade of Lonchodinae, consistent with previous analyses using molecular markers [115,122] and mitochondrial genomes [6,33,91]. Phraortes sp. YW-2014 did not cluster with other Phraortes species, probably because of a species misidentification [33]. Meanwhile, some studies, even including transcriptomes, clearly recover the monophyly of Lonchodinae [5,26,28,114]. Therefore, the problem of the monophyly of Lonchodinae (or not) needs further study. Ph. guangxiensis (Phasmatidae: Clitumninae) formed a sister group to Pharnaciini sp. NS-2020, but the placement of this group and Phobaeticus serratipes was distant from the main clade of Clitumninae. In our work, we failed to recover the monophyly of Clitumninae, similar to the results presented in Bradler et al., Robertson et al.,and Song et al. [5,22,26].
Analyzing the phylogenetic relationships using ML and BI using the 85 species, Embioptera and Zoraptera were clustered into Phasmatodea, and the monophyly of Phasmatodea was not recovered, which was caused by long-branch attraction. After removal of the Embioptera and Zoraptera species, the phylogenetic relationships of ML and BI using the 81 species showed the monophyly of Phasmatodea and the relationships within subfamilies of Phasmatodea were supported.

Conclusions
In this study, we successfully determined the complete mitochondrial genomes of O. guangxiensis, Pe. schultei, and Ph. guangxiensis. The three stick insects shared a similar gene arrangement that has been previously reported for other stick insect species. In this study, after removing representatives of Embioptera and Zoraptera that showed long-branch attraction, BI tree and ML trees showed identical topology, except for the position of Tectarchus ovobessus (Phasmatodea). We recovered the monophyly of Phasmatodea and showed the sister-group relationship between Phasmatodea and Mantophasmatodea. We recovered the monophyly of Heteropterygidae and the paraphyly of Diapheromeridae, Phasmatidae, Lonchodidae, Lonchodinae, and Clitumninae. In this study, Peruphasma schultei (Pseudophasmatidae), Phraortes sp. YW-2014 (Lonchodidae), and species of Diapheromeridae clustered into the clade of Phasmatidae. Within Heteropterygidae, O. guangxiensis was the sister clade to O. mouhotii belonging to Dataminae, and Heteropteryginae was supported as the sister clade to (Dataminae + Obriminae). Future work may need to further explore the mitochondrial genomes of Embioptera and Zoraptera to evaluate the long-branch attraction and explore the phylogenetic relationships between Embioptera and Phasmatodea.  Figure  S5. Phylogenetic relationships of Phasmatodea inferred from ML analysis including Embioptera and Zoraptera. Figure S6. Phylogenetic relationships of Phasmatodea inferred from BI analysis including Embioptera and Zoraptera. Table S1. Specific primers used to amplify the mitochondrial genomes of O. guangxiensis, Pe. schultei and Ph. guangxiensis. Table S2. Species of Phasmatodea used to construct the phylogenetic relationships along with SRA numbers. Table S3. Location of features in the mtDNA of O. guangxiensis. Table S4. Location of features in the mtDNA of Pe. schultei. Table S5. Location of features in the mtDNA of Ph. guangxiensis. Table S6. Mitochondrial genome comparisons of Phasmatodea species. Table S7. The codon numbers and relative synonymous codon usage in mitochondrial protein-coding genes.