Characteristics of the complete mitochondrial genome of Suhpalacsa longialata (Neuroptera, Ascalaphidae) and its phylogenetic implications

The owlflies (Family Ascalaphidae) belong to the Neuroptera but are often mistaken as dragonflies because of morphological characters. To date, only three mitochondrial genomes of Ascalaphidae, namely Libelloides macaronius; Ascaloptynx appendiculatus; Ascalohybris subjacens, are published in GenBank, meaning that they are greatly under-represented in comparison with the 430 described species reported in this family. In this study, we sequenced and described the complete mitochondrial genome of Suhpalacsa longialata (Neuroptera, Ascalaphidae). The total length of the S. longialata mitogenome was 15,911 bp, which is the longest known to date among the available family members of Ascalaphidae. However, the size of each gene was similar to the other three Ascalaphidae species. The S. longialata mitogenome included a transposition of tRNACys and tRNATrp genes and formed an unusual gene arrangement tRNACys-tRNATrp-tRNATyr (CWY). It is likely that the transposition occurred by a duplication of both genes followed by random loss of partial duplicated genes. The nucleotide composition of the S. longialata mitogenome was as follows: A = 41.0%, T = 33.8%, C = 15.5%, G = 9.7%. Both Bayesian inference and ML analyses strongly supported S. longialata as a sister clade to (Ascalohybris subjacens + L. macaronius), and indicated that Ascalaphidae is not monophyletic.


INTRODUCTION
The study of mitochondrial genomes (mitogenomes) is of great interest to many scientific fields, including molecular evolution and evolutionary genomics (Avise et al., 1987;Salvato et al., 2008). Insect mitochondrial genomes are usually a double-stranded circular molecule with a length of 14-20 kbp, including 13 protein-coding genes (PCGs), 22 transfer RNAs (tRNAs), two ribosomal RNAs (rRNAs), and a control region (CR; AT-rich region) (Boore, 1999). The most widespread gene arrangement in insect mtDNAs is the number of sequenced species within the Neuroptera will be very helpful for phylogenetic reconstructions of Neuroptera relationships. Hence, in the present study we sequenced the complete mitogenome of Suhpalacsa longialata Yang 1992 (Neuroptera, Ascalaphidae) and analyzed its genomic structure and composition in comparison with the other three Ascalaphidae species including determining nucleotide composition, gene order, codon usage, and secondary structure of tRNAs. Additionally, we also analyzed evolutionary relationships within Neuroptera using Megaloptera as outgroups to discuss the relationship between Ascalaphidae and Myrmeleontidae, and the relationships of inter-families of Neuroptera.

Sample origin and DNA extraction
The sample of an adult S. longialata used for sequencing was collected from Hangzhou, Zhejiang province, China in July 2017 by LP Zhang. The specimen was identified by JY Zhang and preserved in 100% ethanol at -40 C in the lab of JY Zhang. Total DNA was isolated from one foreleg of S. longialata using Ezup Column Animal Genomic DNA Purification Kit (Sangon Biotech Company, Shanghai, China) according to the manufacturer's protocol.
PCR amplification and sequencing of S. longialata mtDNA A total of 12 universal primers for polymerase chain reaction (PCR) amplification were modified according to Simon et al. (2006), , and Zhang et al. (2018) (Table S1; Fig. 1) based on the mitogenome sequences of the three-known species of Ascalaphidae (L. macaronius, Ascaloptynx appendiculatus, and Ascalohybris subjacens). Then five specific primers (Table S1; Fig. 1) were designed based on the sequence information from universal primers using Primer Premier 5.0 (PREMIER Biosoft International, CA, USA). All PCR was performed with a BioRADMJMini Personal Thermal Cycler (made in Singapore) using Takara Taq DNA polymerase (TaKaRa Biotechnology Co., Ltd., Dalian, China) with the following cycling steps: denaturation at 94 C for 5 min, followed by 35 cycles of 94 C (50 s for denaturation), 48-60 C (30-50 s for annealing), and 72 C (1-3 min elongation), followed by a final elongation at 72 C for 10 min. PCR reactions were carried out in a 50 mL reaction volume consisting of 32.75 mL sterile deionized water, 5.0 mL 10ÂPCR buffer (Mg 2+ Free), 5.0 mL MgCl 2 (25 mM), 4.0 mL dNTP Mixture (2.5 mM each), 1.0 mL DNA template, 1.0 mL each primer (10 ppm), 0.25 mL Takara Taq DNA polymerase (5 U/mL). All PCR products were visualized by electrophoresis in a 1% agarose gel and sent to Sangon Biotech Company (Shanghai, China) for sequencing of both strands.

Phylogenetic analyses
For Megaloptera as a sister clade to Neuroptera proposed by Engel, Winterton & Breitkreuz (2018) and Peters et al. (2014), four species of Megaloptera (Corydalus cornutus; Dysmicohermes ingens; Neochauliodes bowringi; Sialis hamata) (Beckenbach & Stewart, 2008;Cameron et al., 2009;Li et al., 2015;Wang, Liu & Yang, 2016) were used to as  Figure 1 Mitogenome map of S. longialata. The outermost circle shows the gene map of S. longialata and the genes outside the map are coded on the major strand (J-strand), whereas the genes on the inside of the map are coded on the minor strand (N-strand). The middle circle (black) displays the GC content and the paracentral circle (purple & green) displays the GC skew. Both GC content and GC skew are plotted as the deviation from the average value of the total sequence. A total of 17 arcs display the PCR amplification methods. All primers are shown in Table S1.
Full-size  DOI: 10.7717/peerj.5914/ fig-1 outgroups in phylogenetic analyses. We downloaded the data from previously sequenced species of Neuroptera as ingroups including S. longialata (Beckenbach & Stewart, 2008;Cameron et al., 2009;Cheng et al., 2014Cheng et al., , 2015Haruyama et al., 2011;He et al., 2012;Jiang et al., 2017;Lan et al., 2016;Negrisolo, Babbucci & Patarnello, 2011;Wang et al., 2012Wang et al., , 2013Wang et al., , 2017Yan et al., 2014;Zhao et al., 2013Zhao et al., , 2016Zhang & Wang, 2016;Zhang & Yang, 2017) to discuss family-level phylogenetic relationships of Neuroptera. Accession numbers of all mitochondrial genomes are listed in Table S2. Nucleotide sequences of the 13 PCGs were employed for construction of Bayesian inference (BI) and maximum likelihood (ML) phylogenetic trees according to Cheng et al. (2016) and Zhang et al. (2018). DNA alignment was acquired from the amino acid alignment of the 13 PCGs using Clustal W in Mega 7.0 (Kumar, Stecher & Tamura, 2016), and the conserved regions were found by Gblock 0.91b (Castresana, 2000). We estimated the best partitioning scheme and model by the program PartionFinder 1.1.1 (Lanfear et al., 2012) on the basis of Bayesian information criterion. The ML tree was constructed in RAxML 8.2.0 with the best model of GTRGAMMA and the branch support inferred from 1,000 bootstrap replications (Stamatakis, 2014). BI analysis was carried out in MrBayes 3.2 (Ronquist et al., 2012) with the model of GTR + I + G; the analysis was set for 10 million generations with sampling every 1,000 generations; the initial 25% of generations was discarded as burn-in. Because long branch attraction (LBA) can cause a wrong relationship (Bergsten, 2005;Philippe et al., 2005), we obtained a second data set using 40 species of Neuroptera (40SN) as the ingroup by excluding Semidalis aleyrodiformis, Coniopteryx sp., and Dilar sp. that showed LBA. The ML and the BI analyses of data 40SN were then performed as above.
The duplication-random loss model may be a possible explanation for the transposition of contiguous genes. Similar to the report by Beckenbach & Stewart (2008), it is likely that the tRNA Trp -tRNA Cys (WC) genes were duplicated in tandem to form a tRNA cluster WCWC, which was then followed by random loss of partial duplicated genes to produce the final CW gene order. The mitogenome of S. longialata (15,911 bp) is the longest as compared with those of other Ascalaphidae species, whose mitogenomes range from 15,873 to 15,890 bp. The greater length of the S. longialata mitogenome is due largely to 16 intergenic regions ranging from 1 to 54 bp and a long typical A+T-rich region (1,088 bp) as compared to 1,049 bp for L. macaronius (Negrisolo, Babbucci & Patarnello, 2011), 1,066 bp for Ascaloptynx appendiculatus (Beckenbach & Stewart, 2008), and 1,051 bp for Ascalohybris subjacens (Cheng et al., 2014). The nucleotide composition of the S. longialata mitogenome is as follows: A = 41.0%, T = 33.8%, C = 15.5%, G = 9.7%. It is obvious that the S. longialata had a strong A+T bias of 74.8%, which is similar to other species of the Ascalaphidae: 74.5% for L. macaronius; 75.5% for Ascaloptynx appendiculatus; 75.7% for Ascalohybris subjacens (Beckenbach & Stewart, 2008;Negrisolo, Babbucci & Patarnello, 2011;Cheng et al., 2014) (Table 1). The high A+T bias was found in PCGs, rRNA genes, tRNA genes, and the CR. Previous studies pointed out that the strand bias in nucleotide composition may be attributed to mutational damage primarily affecting the lagging strand during asymmetric replication (Francino & Ochman, 1997;Hassanin, Leger & Deutsch, 2005). The skew statistics indicated that S. longialata had a positive AT-skew and negative GC-skew (Table 1).

Protein-coding genes and codon usages
Nine PCGs (ND2, COX1, COX2, ATP8, ATP6, COX3, ND3, ND6, and CYTB) were located on the major strand (J-strand) with the remaining PCGs on the minor strand (N-strand). All PCGs genes used ATN (N represents A, G, C, or T) as initiation codons, which have been accepted as the canonical mitochondrial start codons for insect mitogenomes (Wolstenholme, 1992). Termination codons for S. longialata were mostly complete (TAA) with some incomplete (TA or T). Such incomplete stop codons have been found in various insect species (Ma et al., 2015;Nardi et al., 2001;Fenn, Cameron & Whiting, 2007), and it has been determined that incomplete stop codons can produce functional stop codons in polycistronic transcription cleavage and polyadenylation processes (Ojala, Montoya & Attardi, 1981). The only exception was detected in ND1, where S. longialata exhibited TAG as the stop codon. The infrequent use of TAG may be because of the high A+T composition of the PCGs, although TAG is the conservative stop codon in most insect mitogenomes (Liu et al., 2015). However, in the other three published Ascalaphidae mitogenomes, COX1 of L. macaronius (Negrisolo, Babbucci & Patarnello, 2011), Ascaloptynx appendiculatus (Beckenbach & Stewart, 2008), and Ascalohybris subjacens (Cheng et al., 2014) used ACG as the start codons, and ND1 of Ascalohybris subjacens used TTG. The other start/stop codons were identical to the S. longialata situation. The total length of the 13 PCGs in the S. longialata mitogenome was 11,169 bp, with an average AT content of 73.0%. The PCGs displayed A-skews (A > T) and C-skews (C > G) (Table 1). We calculated the RSCU of the S. longialata mitogenome, excluding stop codons (Fig. 2). The RSCU proved that codons with A or T in the third position are always overused when compared to the other synonymous codons. The codons of amino acids being NNW (NNA/NNU) were higher than 1.0 without exception in S. longialata. The most frequently encoded amino acids were Leu (UUR), Phe, Ile (>300), and the least frequently used amino acid was Cys (<45) (Table S3), which was similar to the other Ascalaphidae mitogenomes (Fig. 2).

Ribosomal and transfer RNAs
The mtDNA of S. longialata contained the entire content of 2 rRNAs and 22 tRNAs genes that were also found in other neuropterid mitogenomes (Boore, 1999;Song, Lin & Zhao, 2018;Wang et al., 2017). The 16S rRNA gene with a length of 1,314 bp was located between tRNA Leu (CUN) and tRNA Val whereas the 12S rRNA gene with a size of 739 bp was located between tRNA Val and the CR; these locations were also detected in the other ascalaphid owlfly species (Beckenbach & Stewart, 2008;Negrisolo, Babbucci & Patarnello, 2011;Cheng et al., 2014). The AT content of rRNAs in the S. longialata mitogenome was the highest (77.8%) except for the A+T-rich region (85.1%). We found that the AT-skew was strongly positive whereas the GC-skew was highly negative, which showed that the contents of A and C were higher than those of T and G, respectively.
The size of the tRNAs was 1,476 bp with an average A+T content of 76.2%. Among the 22 tRNAs, most tRNA genes displayed the common cloverleaf secondary structure, whereas the tRNA Ser(AGN) had lost the dihydrouridine (DHU) arm (Fig. 3). The absence of this arm in tRNA Ser(AGN) is a typical feature of many insect mtDNAs (Wolstenholme, 1992;Salvato et al., 2008;Sheffield et al., 2008;Negrisolo, Babbucci & Patarnello, 2011;Yan et al., 2014;Du et al., 2017;, and is usually demonstrated to be functional (Hanada et al., 2001;Stewart & Beckenbach, 2003). We also found that the tRNA Phe and tRNA Leu (CUN) lack the TcC loops. Furthermore, unmatched U-U base pairs were observed in tRNA Trp (Fig. 3).
In terms of the tRNA gene structures of the other three ascalaphid owlflies, the tRNA Phe in L. macaronius and Ascalohybris subjacens showed the loss of TcC loops, and the tRNA Ser (AGN) in Ascalohybris subjacens lost the DHU loop, whereas the tRNA genes of Ascaloptynx appendiculatus displayed the typical cloverleaf secondary structure.

A+T-rich region and intergenic regions
Generally speaking, the A+T-rich region was the largest non-coding region, which was located between 12S rRNA and tRNA Ile . The A+T-rich region of S. longialata mtDNA having a length of 1,088 bp was the longest when compared to the other three species of Ascalaphidae, for example, the L. macaronius (1,049 bp), Ascaloptynx appendiculatus (1,066), and Ascalohybris subjacens (1,051 bp). Additionally, the composition of A+T was 85.1% in S. longialata, which was higher than in L. macaronius (84.5%) and lower than Ascaloptynx appendiculatus (85.7%) and Ascalohybris subjacens (86.2%).
The mitochondrial genomes of most insects are compact (Boore, 1999), although large intergenic regions occur in some species. In the S. longialata mitogenome the longest intergenic region was a 54 bp insertion between tRNA Ile and tRNA Gln . This spacer was also present in L. macaronius, Ascaloptynx appendiculatus, and Ascalohybris subjacens and spanned 55, 42, 54 bp, respectively (Beckenbach & Stewart, 2008;Negrisolo, Babbucci & Patarnello, 2011;Cheng et al., 2014). This intergenic region of the four species also shared a 12 bp long congruent motif A(A/G)TTAA(A/C)TAAAT adjacent to tRNA Gln . It has previously been reported that this spacer may diverge quickly among different families of the same order (Negrisolo, Babbucci & Patarnello, 2011). Aside from this spacer, gaps between genes ranged from 1 to 18 residues in the S. longialata sequence.

CONCLUSION
We successfully sequenced the entire mitochondrial genome of S. longialata, which showed similar gene characteristics to the other three species of Ascalaphidae. Both BI and ML analyses supported S. longialata as a clade sister to (Ascalohybris subjacens + L. macaronius), but Ascalaphidae is not monophyletic. From the results obtained in the present study, we believe the different topologies of phylogenetic relationships were caused mainly by LBA of Coniopteryx sp., Dilar sp., and Semidalis aleyrodiformis.