Phylogenetic relationships among tribes of the green lacewing subfamily Chrysopinae recovered based on mitochondrial phylogenomics

Chrysopidae (green lacewings) is the second largest family in Neuroptera, and it includes medium-size lacewings largely recognized by the presence of golden-colored eyes, bright green bodies and delicate wings with dense venation patterns. The subfamily Chrysopinae includes 97% of the species diversity in the family and it is currently divided into four tribes: Ankylopterygini, Belonopterygini, Chrysopini and Leucochrysini. Here we sequenced and annotated the nearly complete mitochondrial genomes of four species of each these tribes: Abachrysa eureka, Italochrysa insignis, Leucochrysa pretiosa, Parankyloteryx sp. We then reconstructed the phylogenetic relationships with estimated divergence times among tribes of Chrysopinae based on the mt genomic data. Our results suggest that Chrysopinae sans Nothancyla verreauxi evolved as two reciprocally monophyletic lineages formed by stem members of the tribes Leucochrysini plus Belonopterygini on one hand, and the stem members of Ankylopterygini plus Chrysopini on the other. Our estimations of divergence times place the diversification of stem Chrysopinae into the extant tribes during the Middle Jurassic to Late Cretaceous. The relatively young ages previously estimated for the green lacewing divergences were probably underestimated due to false inferences of homology between non-sister taxa that are later correctly identified as homoplasy after more taxa are added.

In the case of Neuroptera, sequencing and annotation of mitogenomes began almost 10 years ago 9 and today there are 40 published mitogenomes; this achievement recently allowed the formulation of the first hypothesis for Neuropterida based exclusively on the phylogenetic signal provided by mitogenomes 8 . Likewise within Chrysopidae, mitogenomes have proven to be informative, as they helped place at last the enigmatic monotypic genus Nothancyla verreauxi, whose phylogenetic affinities to the other subfamilies of Chrysopidae remained ambiguous for many years 7 .
So far the phylogenetic hypotheses proposed for Chrysopidae vary greatly among authors; this is particularly true for the relationships among the tribes of Chrysopinae, which remain largely equivocal. Persistent questions regarding tribal relationships within Chrysopinae include the possible paraphyly of Chrysopini 5, 6 , the generic membership and delimitation of Leucochrysini and Belonopterygini 10,11 and the position of Nothancyla Navás within Chrysopinae 4,7,12 .
In this study, we provide the results of sequencing and annotating the mitochondrial genomes of four species of the subfamily Chrysopinae. These species each belong to Ankylopterygini, Belonopterygini and Leucochrysini, chrysopinae tribes whose mt genomes have not been sequenced previously. We use these and the mitochondrial genomes previously sequenced for five species of Chrysopidae 7, 13, 14 to elucidate the higher-level phylogeny and estimation of divergence times among chrysopine tribes of Chrysopinae.

Results and Discussion
Gene structure. The sequence of Abachrysa eureka (Banks, 1931) (Belonopterygini) is 14500 bp in size, containing genes, with 13 PCGs, 19 tRNAs and 2 rRNAs (Fig. 1), while that of Italochrysa insignis (Walker, 1853) (Belonopterygini) is 14433 bp in length consisting the same 33 complete genes and partial sequence of nad2. The arrows indicated the orientation of gene transcription. The tRNAs are denoted by the color blocks and are labelled according to the IUPACIUB single-letter amino acid codes (L1: UUR; L2: CNU; S1: AGN; S2: UCN). The GC content was plotted using a black sliding window, as the deviation from the average GC content of the entire sequence. GC-skew was plotted as the deviation from the average GC-skew of the entire sequence. The inner cycle indicated the location of the genes in the mt genome.
SCientifiC REPORTS | 7: 7218 | DOI:10.1038/s41598-017-07431-1 The sequences of Leucochrysa (Leucochrysa) pretiosa (Banks, 1910) (Leucochrysini) and Parankylopteryx sp. (Ankylopterygini) were 14185 bp and 14364 bp respectively, containing the same 32 complete genes with A. eureka and partial sequences of nad2 and rrnS. Because of the high A + T composition and complicated secondary structure, we were unable to amplify the control regions as well as the genes approaching it in the four species. Gene orders of all species sequenced are identical and as a group, they all differ from the ancestral insect mt genome arrangement 15 , at the translated positions of trnC and trnW, which are translocated. As shown previously by multiple studies, this translocation is universal and exclusive to Neuroptera 7,8,13,[16][17][18] , the exceptions are Coniopterygidae, Nevrorthidae, Sisyridae and Osmylidae, whose mt genomes exhibit the ancestral arrangement with no translocation of trnC and trnW 8,19 .
There are overlaps between every two genes in the mt genomes we sequenced. Two gene pairs including atp8-atp6 and nad4-nad4L overlap 7 nucleotides, also sharing the sequence ATGNTAA which has been reported in many other insect mt genomes with difference appeared as ATGATAA and ATGTTAA respectively 20,21 . In addition, there are small non-coding intergenic spacers upstream and downstream of nad1 about 20-30 bp in length.
Comparison of tRNA structure. There are 19 tRNAs found in the partial mt genome sequences of the four species sequenced here. The length of the tRNAs ranged from 62 to 72 bp. Twelve genes are located on the J-strand and the other seven are located on the N-strand. Most of the tRNAs can be folded into the typical clover-leaf structure (Fig. 2) except trnS1 (AGN), with its dihydorouridine (DHU) arm forming a simple loop, which is a common phenomenon in sequenced Neuropterida mt genomes 7,22 . The aminoacyl (AA) stem (7 bp) and the anticodon loop (7 bp) are conservative in length. The length of the anticodon (AC) stems is constantly 5 bp, with exception of trnK, where a mismatch is present in the first pair of their anticodons. The DHU and TΨC (T) stems are really variable in the loops and stems (Fig. 2). Based on the secondary structure models, mismatched base pairs (i.e. U-U and A-C pairs) were found in the tRNAs. Nucleotide substitutions rarely occur on acceptor and anticodon stems. The changes are fully compensatory base changes (cbcs) (e.g. G-C vs A-T on the acceptor stem of trnR and trnP) or hemi-cbc (e.g. G-T vs A-T on the acceptor stem of trnT and trnW) 23 , which are more restricted than that on TΨC (T) and DHU loops. In addition, substitutions involving two full cbcs are found in the trnT anticodon stem and the trnR acceptor stem.
In previous studies of insect phylogeny, tRNAs have rarely been used. However, the presence of full cbcs and hemi-cbcs may characterize taxa at different taxonomic levels, as well as the secondary structure of tRNAs 24 . The base substitution in stems and loops may provide some evidence in reference to the phylogeny. Sequences of the variable tRNAs of Parankylopteryx sp. (Ankylopterygini) are quite different from that in the other three species we sequenced (Fig. 2), e.g. trnR and trnS2 (UCN). In addition, among these variable tRNAs, the acceptor stems of trnL2 (CUN) and trnP possess T-A pair in Parankylopteryx sp., while in the other taxa the C-G pair is present instead. The sequences of tRNAs in Abachrysa eureka (Belonopterygini) and Italochrysa insignis (Belonopterygini) appear to be invariable, and even identical in trnG and trnV. The sequences of tRNAs in Leucochrysa pretiosa (Leucochrysini) are more similar to that of Abachrysa eureka than Parankylopteryx sp. (see Fig. 2). This supports the result of phylogenetic analysis that Belonopterygini is the sister group of Leucochrysini, while Ankylopterygini is distantly related to them.
Comparison of rrnL secondary structures. The lengths of rrnL of the four green lacewings we sequenced were determined from 1301 bp to 1318 bp (see bioinformatics analysis). The inferred secondary structure model for the rrnL of A. eureka is provided in the Fig. 3. The structure of rrnL largely resembles previously published structures for Libelloides macaronius 17 and Nothancyla verreauxi Navás 7 . The secondary structure consists of five canonical domains (I-II, IV-V) with domain III absent, which is a typical trait in arthropods 25 (Fig. 3), including 50 helices. The predicted secondary structure of domain I includes eight helices and is consistent among all nine chrysopid taxa. There are conserved sites unevenly distributed throughout the rrnL secondary structure. For example, the sequence of rrnL of A. eureka is in a visible uniformity to 69.01% with the multiple alignments of rrnL in Chrysopidae and 71.84% in Chrysopinae. The region with highest level of invariable positions is located in domain IV, especially the triangle region formed by helices H30 to H34, while that with lowest level of invariable positions is observed respectively in domains I-II.
Phylogenetic relationships among tribes of Chrysopinae. Bayesian and maximum likelihood analyses generated phylogenetic trees with the same topology as the chronogram in Fig. 4. The sampling list of the present phylogenomic analysis is shown in Table 1. Accordingly, Chrysopinae was recovered sister to Apochrysinae plus Nothochrysinae, which is consistent with results of previous studies using nuclear genes only 5,6 . The sister relationship between Nothochrysinae and Apochrysinae recovered here had a high bootstrap value and strong posterior probability support resulted under partition-specific models, but with rather low posterior probability support resulted under CAT model, which indicates instability concerning the relationships between these subfamilies. A previous mitogenomic analysis for the family 7 found Apochrysinae sister to the remaining Chrysopidae, a result that agrees partially with hypotheses based on morphology 1,26 . Taken together these results suggest that the signal provided by mitogenomes seems equivocal and unable to resolve relationships in the deepest part of the Chrysopidae tree.
As a genus that exhibits characteristics of both Apochrysinae and Chrysopinae, Nothancyla was recovered here as the sister group to all other Chrysopinae, in agreement with the result from a previous study using fewer mt genomes 7 . The remaining Chrysopinae are grouped into two clades, Leucochrysini + Belonopterygini and Ankylopterygini + Chrysopini, which is in accordance with the relationships proposed by Brooks 3 . Leucochrysini and Belonopterygini share certain synapomorphic genitalic characters suggesting a closer relationship between the two tribes. For example, the gonarcus is broad, transverse and hardly arcuate with gonocornua, and the arcessus is short and broad with a strong hook and lateral membranous lobes at the apex in both tribes. Nevertheless, Belonopterygini is characterized by the long and broad intramedian cell, which is usually quadrangular. The synapomorphy supporting the sister group relationship between Ankylopterygini and Chrysopini is the narrow arcessus. Ankylopterygini genera can be distinguished by the narrow hind wings (length:breadth > 3.3:1) 1 , densely setose forewings with broad costal region, scythe-like mandibles and apically constricted maxillary palps.
Remarkably, the relationships among tribes of Chrysopinae in this study are different from other results from DNA data with weak support in various topologies. Chrysopini was recovered as sister to the remaining Chrysopinae the other three tribes in the study by Winterton and de Freitas 4 , with Ankylopteryigini sister to Leucochrysini. Moreover, subsequent studies using different molecular datasets did not support the monophyly of Chrysopini 5, 6 . A recurring issue with molecular phylogenetic studies of Chrysopidae is either limited character sampling 13 or limited taxon sampling 4,7 . This study again presents large amount of DNA sequence data, but also begins to build upon the limited taxon sampling of previous studies to include representatives of all major chrysopid lineages. The next step is to expand on this taxon sampling considerably, to begin testing the monophyly of these lineages in a phylogenetic context. The larval stages of chrysopids exhibit the unusual behavior is entangling debris in elongate, recurved setae on their dorsum to form a trash packet. This is used for camouflage and as a physical shield against predators and parasites 27 . However, not all green lacewing larva carry a trash packet and those that do not (so called 'naked' larvae) display different behaviors (i.e. generally nocturnal, faster moving). Trash carrying is not found in larvae of the closely related Hemerobiidae. While Nothochrysa larvae do carry a trash packet, most Nothochrysinae do not. It is found in most representative sampled in Chrysopinae, except Chrysoperla and Chrysopa. These results suggest two evolutionary scenarios that (1) larval trash carrying evolved once in stem chrysopids and then has  Letters on nodes correspond to the age estimates featured in Table 2. Scale units are in millions of years. Stars indicate the node with fossil calibration and letters are used as arbitrary labels.
been lost at least twice, in both Nothochrysinae and in at least some Chrysopini; and (2) it evolved twice in stem Chrysopinae and within Nothochrysinae and then has been lost once in at least some Chrysopini.

Divergence time estimates. The chronogram with divergences times estimations obtained with
PhyloBayes using fossil calibrations is featured in Fig. 4. Mean age values and 95% high posterior density (HPD) intervals for each node are presented in Table 2. The estimated age at the node on the chronogram represents a mean of the probability distribution of ages for that node, with time intervals within 95% probability of the mean age. Estimations of divergence times based on the concatenation of the 13 PCGs and rRNAs placed the divergence between stem Nothochrysinae + Apochrysinae from stem Chrysopinae in the Late Triassic to the Early Jurassic (~201 MA). In addition, stem Apochrysinae diverged from stem Nothochrysinae at ~149 MA, in the end of Jurassic to onset of Cretaceous. While the first split within Chrysopinae was that stem Nothancyla verreauxi diverged from stem Chrysopinae around 193 MA. The diversification of Chrysopinae into two reciprocally monophyletic lineages occurred approximately 185 MA. The tribes Leucochrysini was estimated to be the oldest of the four tribes, diverging from Belonopterygini in the Middle Jurassic (~171 MA), in contrast to Ankylopterygini and Chrysopini which diverged from their common ancestor close to the end of Cretaceous at ~81 MA.
When comparing the divergence times estimations from this study with those from that of Dai et al. 7 (based also on mitogenomes but with a smaller taxonomic sampling), we found that the divergence time estimates obtained here are significantly older than those obtained previously. This surprising result (given the similarities in the sourced data and calibrations used by both studies) is presumably due to the "node density effect" 28 on branch lengths, which results from the detection of homoplasy, initially mistaken as homology, once more taxa is included. Homoplasy erroneously inferred as synapomorphies when less taxa are included, will make unrelated taxa look more similar and (under a clock) more recently diverged than they are.
As deep phylogenetic divergences were included in this study, theoretically given the known rapid nucleotide substitutions of mt DNA, a strict molecular clock model failed (or underestimate) at accounting for the whole extend of the sequence divergence because of saturation (aka. Multiple hits). The use of a strict molecular clock model does not leave room to express uncertainty around the estimations, as it provides direct measurements  Table 2. Mean divergence times and 95% high posterior density (HPD) intervals for each node of the topology presented in Fig. 4. Time-scale units are in millions of years.
resulted from the application of the rate. In addition, previous research suggested that there is noticeable rate variation among lineages so that the use of a strict molecular clock model could not account for such rate variation 8 . As a result of this, the divergence time estimated using a strict molecular clock model was too young to be believable (see Table 2 for mean age values and 95% HPD intervals for each node). The fact that branch lengths estimations are affected by taxon sampling is a problem for relaxed molecular clocks that has been reported multiple times [29][30][31] . Therefore, this means that as old as the ages obtained here with a larger taxonomic sample are, they represent a better hypothesis than the one obtained by Dai et al. 7 . Of course all things being equal (e.g. calibrations), we would expect future divergence times estimations for Chrysopidae to be even older if they are estimated for a more inclusive taxon sampling.

Materials and Methods
Specimens and DNA extraction. Specimens from four species belonging to the Chrysopinae tribes: Belonopterygini Navás, 1913, Leucochrysini Adams, 1978, andAnkylopterygini Navás, 1910 and Parankylopteryx sp. were deposited in GenBank with accession numbers KY587199, KY587200, KY587201 and KY587202, respectively. Sequences assembly was done by using ContigExpress. tRNAs were identified by tRNAscan-SE Search Server v.1.21 32 . The secondary structure of the tRNAs were generated by tRNAscan-SE as stated, while some of them were predicted by the relevant genes from green lacewing species Nothancyla verreauxi Navás and Nothochrysa californica Banks 7 which could not be generated by tRNAscan-SE. The PCGs were identified as open reading frames corresponding to the 13 PCGs in metazoan mt genomes, whose boundaries were identified by comparing the boundaries of the same PCGs in related neuropteran mt genomes with MEGA 5.0 33 . It is impossible to infer the boundaries of rRNAs accurately because there are no start or stop codons in the rRNA genes. Therefore, the rRNAs were assumed to extend to the boundaries of flanking genes 34,35 and were confirmed by alignment to relevant genes of other species of Neuroptera (i.e., Apochrysa matsumurae Okamoto 13 , Nothancyla verreauxi Navás and Nothochrysa californica Banks 7 using MEGA 5.0. The control regions were identified afterwards by the boundary of the rRNAs genes and compared with other insect mt genomes. Alignment of the sequences of the 13 PCGs was inferred from the amino acid sequences using ClustalW in MEGA 5.0 33 , while the rRNA alignments were conducted by G-blocks Server (http://molevol.cmima.csic.es/ castresana/Gblocks_server.html). Then, alignments of individual genes were concatenated without stop codons but including all three codon positions for phylogenetic analysis. Bayesian (BI) and maximum likelihood (ML) analyses were made based on the partitioned datasets with partition-specific models (Table S5-6) estimated by PartitionFinder 37 . The BI inference was performed in MrBayes Version 3.1.2 38 , in which two simultaneous runs of 2,000,000 generations were conducted, and the dataset was sampled every 200 generations with a burn-in of 25%. In the ML algorithms, nodal supporting values were assessed by bootstrap values (BP) 39 calculated with 10000 replicates from 1000 RAxML runs. Also phylogenies were inferred under the heterogeneous model CAT using PhyloBayes 3.3 40 . Four independent tree searches were run until the likelihood of the sampled trees was stabilized and every two runs had satisfactorily converged (maxdiff less than 0.3).

Divergence time estimation.
Divergence time among the green lacewing subfamilies was estimated in PhyloBayes 3.3 40 . The topology obtained in each of the phylogenetic analysis earlier was constrained using the following five minimum age calibrations, (1) the most recent common ancestor (MRCA) of Italochrysa and Abachrysa was calibrated with an undescribed specimen from the late Eocene in Baltic Amber to a minimum age of 38 Ma 41 ; (2) the split of Leucochrysa from Italochrysa + Abachrysa was calibrated with the mean age of SCientifiC REPORTS | 7: 7218 | DOI:10.1038/s41598-017-07431-1 Leucochrysa (Nodita) prisca to a minimum age of 21 Ma 42 ; (3) the split of Chrysoperla was constrained to a min age of 21 Ma based on Chrysopa glaesaria 42 ; (4) the MRCA of Apochrysinae and Nothochrysinae was calibrated with Adamsochrysa aspera to a min age of 53 Ma 43 ; (5) the MRCA of Hemerobiidae and Chrysopidae was calibrated with a mean age of 163 Ma based on Mesypochrysa 43 . Additionally, a secondary calibration of 260 Ma for the root was placed based on the divergence of Chrysopidae + Hemerobiidae from Mantispidae. In PhyloBayes two chains of autocorrelated lognormal relaxed clock model 44 were run, using a birth-death process as tree prior for 28 000 cycles with the first 5000 cycles 40 discarded as burn in.