Insight into higher-level phylogeny of Neuropterida: Evidence from secondary structures of mitochondrial rRNA genes and mitogenomic data (original) (raw)

• Journal List
• PLoS One
• PMC5790268

PLoS One. 2018; 13(1): e0191826.

Nan Song, Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing,* Aili Lin, Formal analysis, and Xincheng Zhao, Formal analysis, Methodology

Nan Song

College of Plant Protection, Henan Agricultural University, Zhengzhou, China

Aili Lin

College of Plant Protection, Henan Agricultural University, Zhengzhou, China

Xincheng Zhao

College of Plant Protection, Henan Agricultural University, Zhengzhou, China

Bi-Song Yue, Editor

College of Plant Protection, Henan Agricultural University, Zhengzhou, China

Sichuan University, CHINA

Competing Interests: The authors have declared that no competing interests exist.

Received 2017 May 11; Accepted 2018 Jan 11.

This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Supplementary Materials

S1 Fig: Secondary structure drawing of the complete mitochondrial rrnS gene for ten neuropterid species. (A) Acanthacorydalis orientalis, (B) Ascalohybris subjacens, (C) Corydalus cornutus, (D) Dysmicohermes ingens, (E) Micromus angulatus, (F) Mongoloraphidia harmandi, (G) Neochauliodes fraternus, (H) Thaumatosmylus sp., (I) Rapisma sp., and (J) Sialis hamate.

(ZIP)

GUID: 7FA3B7A0-B569-45A9-9454-8B58D866D203

S2 Fig: Secondary structure drawing of the complete mitochondrial rrnL gene for ten neuropterid species. (A) Acanthacorydalis orientalis, (B) Ascalohybris subjacens, (C) Corydalus cornutus, (D) Dysmicohermes ingens, (E) Micromus angulatus, (F) Mongoloraphidia harmandi, (G) Neochauliodes fraternus, (H) Thaumatosmylus sp., (I) Rapisma sp., and (J) Sialis hamate.

(RAR)

GUID: 1A083371-83CB-4304-AF23-A32FAAB283D2

S3 Fig: Phylogenetic reconstruction from the data set of _rrnS-_SS using software IQ-TREE under GTR+G+I model. Node values represent bootstrap values.

(TIF)

GUID: 719DF29C-B52F-4FD2-BEE6-A5D227E92264

S1 Table: Taxa used in this study. (XLSX)

GUID: 2A8178E3-764B-414C-AF72-714A1A042CD9

S2 Table: The partition schemes and best-fitting models selected by PartitionFinder for each dataset. (XLSX)

GUID: 40120264-F74C-4F50-BDAB-4C1457F38E60

S3 Table: The organization of four newly sequenced mitogenomes. (XLSX)

GUID: 87C26D35-C5B6-41CA-95B8-3386A78CDCB4

S4 Table: The Chi-square tests. (XLSX)

GUID: 592BC868-BBF1-4D96-BD04-F6B4709F7C0D

S5 Table: The evolutionary rate analyses conducted for each species and for each gene partition. (XLSX)

GUID: 5B89A7A3-8ACF-49C2-BDD8-C1C442FA21EC

S6 Table: The nucleotide composition in each helix of the secondary structures of rrnS in 11 Neuropterida mitogenomes. (XLSX)

GUID: 7EECD126-8F0E-4D2B-A785-6355F4197135

S7 Table: The nucleotide composition in each helix of the secondary structures of rrnL in 11 Neuropterida mitogenomes. (XLSX)

GUID: B44D5BEC-DFE8-4E36-B93C-127A9EBC1C27

S1 File: All sequence alignment files in PHYLIP format compiled in this study. (ZIP)

GUID: 3D5B9C47-756F-4E6A-97F1-8372A00E6974

S2 File: All ML trees in NEXUS format reconstructed in this study. (ZIP)

GUID: CD1FF5D1-1AF7-4BCB-800E-C7287B9F6EFD

S3 File: All Bayesian trees in NEXUS format reconstructed in this study. (ZIP)

GUID: ADFDF172-1EC5-481D-8696-0C0BFD17C9C8

Data Availability Statement

All relevant data are within the paper and its Supporting Information files.

Abstract

It is well known that the rRNA structure information is important to assist phylogenetic analysis through identifying homologous positions to improve alignment accuracy. In addition, the secondary structure of some conserved motifs is highly stable among distantly related taxa, which can provide potentially informative characters for estimating phylogeny. In this paper, we applied the high-throughput pooled sequencing approach to the determination of neuropteran mitogenomes. Four complete mitogenome sequences were obtained: Micromus angulatus (Hemerobiidae), Chrysoperla nipponensis (Chrysopidae), Rapisma sp. (Ithonidae), and Thaumatosmylus sp. (Osmylidae). This allowed us to sample more complete mitochondrial RNA gene sequences. Secondary structure diagrams for the complete mitochondrial small and large ribosomal subunit RNA genes of eleven neuropterid species were predicted. Comparative analysis of the secondary structures indicated a closer relationship of Megaloptera and Neuroptera. This result was congruent with the resulting phylogeny inferred from sequence alignments of all 37 mitochondrial genes, namely the hypothesis of (Raphidioptera + (Megaloptera + Neuroptera)).

Introduction

Neuropterida is a name sometimes applied to an insect superorder within Holometabola. This superorder comprises three orders: Raphidioptera (snakeflies), Megaloptera (dobsonflies) and Neuroptera (lacewings). Some insect species of Neuropterida are economically important, owing to their significant role played by adults and/or larvae in the bio-control of insect pest species on agricultural crops [1,2]. Moreover, some species of Neuropterida exhibit an exceptionally wide range of morphological and biological diversity. Approximately 200 described species have been placed in two families of Raphidioptera, about 300 species in two families of Megaloptera, and about 5,800 species in 15 or 16 families of Neuroptera [3]. Furthermore, it is estimated that there are 10,000 extant neuropteridan insect species in the world [4–6]. Therefore, the Neuropterida are recognized as the fifth largest assemblage among the holometabolan groups [7].

Although with the economic importance and relatively rich fossil evidence, the phylogenetic position of the entire Neuropterida and interrelationships among the three orders of Neuropterida have been the most contentious questions in the higher-level phylogenetics of Holometabola. The Neuropterida is commonly considered to be one of the most basal lineage within Holometabola [6]. The Coleoptera is typically identified to be the sister group of Neuropterida [7–10]. However, Boudreaux (1979) [11] suggested the Mecopterida (including Mecoptera, Siphonaptera, Diptera, Trichoptera and Lepidoptera) to be the sister group of Neuropterida. Alternatively, the groupings of (Mecoptera + Siphonaptera) [12], ((Lepidoptera + Trichoptera) + (Mecoptera + Siphonaptera)) [13], (Strepsiptera + (Coleoptera + Diptera)) [14], and (Strepsiptera + Coleoptera) [15] were respectively proposed as the sister group of Neuropterida. More recently, a genome-scale study recovered the clade (Strepsiptera + Coleoptera) as the sister group of Neuropterida [16]. In our prior analyses of the entire Holometabola phylogeny based on mitogenomic data, the sister group relation between the Neuropterida and the clade (Strepsiptera + Coleoptera) was the mostly recovered relationship [17]. Another recent mitogenomic study also supported the same sister group relationship of Neuropterida and the clade (Strepsiptera + Coleoptera) [18]. In short, the majority of recent studies tended to reveal a close affinity of Neuropterida with Coleoptera.

In contrast to the unstable phylogenetic position of Neuropterida, the monophylies of the entire Neuropterida and of both Raphidioptera and Neuroptera have been generally accepted. Morphological evidence supporting the monophyly of Neuropterida have been summarized by Aspöck et al. (2012) [19]. Furthermore, the monophyletic Neuropterida was confirmed by several molecular studies [4,16–18]. Recently, the monophyly of Neuropterida was supported by a phylogenomic study based on anchored hybrid enrichment data [20]. Some earlier studies indicated Megaloptera as a non-monophyletic assemblage [15,21,22]. However, much more researches supported Megaloptera to be monophyletic [4,16–18,20,23–25]. For the inter-ordinal relationships within Neuropterida, three possible arrangements have been proposed. Traditionally, the sister group of (Megaloptera + Raphidioptera) was hypothesized in the context of morphological analyses [8,9,21]. Song et al. (2016) recovered a clade of (Neuroptera + (Megaloptera + Raphidioptera)) based on the mitogenomic data [18]. This relationship was also recovered by 18S rRNA data [26]. Wheeler et al. (2001) [27] provided an unresolved trichotomy of the three orders within Neuropterida based on the combined analysis of nuclear 18S rRNA and 28S rRNA gene sequences. In contrast, several authors provided good support for the sister-group relationship between Megaloptera and Neuroptera based on complete mitogenome data [7,28] and genome-wide data [16,20]. This relation was corroborated by morphological data of genital sclerites [23] and wing base structures [29].

As mentioned above, complete mitogenomes have been used to infer Neuropterida relationships at different phylogenetic levels [7,17,18,25,28,30], with various degree of success. However, these studies were based on the traditional Sanger sequencing, by which reconstructing mitogenomes were time and cost consuming due to the difficulty of long PCR and the necessity of designing a large set of species-specific PCR primers for primer-walking [31]. The advent of next-generation sequencing (NGS) technologies is revolutionizing the current biology. Specifically, NGS approaches have made insect mitogenome reconstruction to be less labor-intensive and more cost-effective. Several recent studies have successfully employed the NGS methods to assemble mitogenome sequences, which allow for multiplexing samples from diverse taxa and maximizing the number of taxa involved in the research [32–36].

The relative ease of acquiring complete mitogenomes by NGS methods further contributes to the mitogenomic phylogeny researches. Of the 37 mitochondrial genes, the mitochondrial small (rrnS) and large (rrnL) ribosomal subunit RNA sequences have been the most widely used phylogenetic markers [37–45]. The total length of both rrnS and rrnL gene sequences of neuropterid species is approximately 2.1 kb, which makes up ca. 14% phylogenetic information of the whole mitogenome. However, explicitly aligning these mitochondrial ribosomal RNA sequences is generally difficult, because of the variable sequence length and the indel-rich regions contained in both gene sequences. Inexactly aligned nucleotide sequences do not correctly reflect homology between sequences [46–48], which might lead to artificial phylogenetic hypotheses. Several previous studies have shown that the secondary reconstruction of ribosomal RNA molecules can aid in improving recognition of primary homology [46,51,52], and refine the alignment process [46–50]. Some structural motifs are highly stable among distantly related taxa, which can provide potentially informative characters for estimating phylogeny [53].

In the present study, we applied high-throughput pooled sequencing approach to reconstruct four complete mitogenomes of Neuroptera, and used them to investigate the interrelationships among Raphidioptera, Megaloptera and Neuroptera. In addition, secondary structures of the complete mitochondrial rrnS and rrnL genes were presented for eleven neuropterid species to provide new evidence for the phylogeny of Neuropterida.

Materials and methods

Ethics statement

No specific permits were required for the insect specimens collected for this study in China. The specimens of Micromus angulatus, Chrysoperla nipponensis EMHAU-15090613, and Thaumatosmylus sp. were collected by authors in Xinyang, Henan, China, while the Rapisma sp. was collected Linzhi, Tibet, China. The field studies did not involve endangered or protected species. The species sequenced in this study is not included in the ‘‘List of Protected Animals in China”.

Taxon sampling

Tissue samples from four species of Neuroptera, including M. angulatus (Hemerobiidae), C. nipponensis EMHAU-15090613 (Chrysopidae), Rapisma sp. (Ithonidae) and Thaumatosmylus sp. (Osmylidae) were collected for DNA extraction and subsequent sequencing.

In the phylogenetic analysis, a total of 35 taxa were included for tree reconstruction. Two species of Coleoptera were selected as the outgroup taxa. Ingroups included 33 species representing three orders of Neuropterida. S1 Table lists the taxa analyzed in this study.

De novo assembly of mitogenomes

Genomic DNA was extracted from 95% ethanol-preserved specimen individually from the thoracic leg muscle tissue using the TIANamp Micro DNA Kit (TIANGEN BIOTECH CO., LTD) as per the manufacturer’s protocol. DNA concentration was measured by Nucleic acid protein analyzer (QUAWELL TECHNOLOGY INC.).

Uniform quantities of genomic DNA from each of the neuropteran insects were pooled, and DNA was concentrated to 1.5 μg. The mixed DNA sample was utilized for library construction, using Illumina TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, CA, USA) and with an average insert size of 350 bp. The subsequent de novo genome sequencing were conducted on the Illumina HiSeq2500 platform in Shanghai OE Biotech CO., LTD. For the sequenced sample, 10 Gb paired-end reads of 125 bp length were generated. FastQC [54] was used for the quality control checks on raw sequence data. NGS QC toolkit [55] was applied to filter the data (cutoff read length for HQ = 70%, cutoff quality score = 20). In this step, reads containing adapters and ploy-N, and low quality reads were removed from raw data. At the same time, Q20, Q30, GC-content and sequence duplication level of the cleaned data were calculated. All the downstream analyses were based on clean data with high quality (avg. Q20 > 90%, and avg. Q30 > 85%). Finally, no less than 8 Gb high-quality reads were used in de novo assembly with IDBA-UD v. 1.1.1 [56]. The assemblies were constructed using 200 for the setting of minimum size of contig, and an initial k-mer size of 40, an iteration size of 10, and a maximum k-mer size of 90.

To identify the mitogenome assemblies from the pooled sequencing files, three different fragments of mtDNA (cox1 5’ region, cytb 3’ region, and rrnS 5’ region) were amplified as “Bait” sequences by standard PCR reactions using universal primers designed by the study of Song et al. (2016) [57]. Local BLAST searches were conducted with BioEdit version 7.0.5.3 [58] for each bait sequence reference against all corresponding assemblies. Only hit with 100% pairwise identity was considered as a successful identification. The mitochondrial contigs identified were inputted into MITOS [59] for initial mitogenome annotation. The resultant gene boundaries were checked and corrected by alignment with published Neuroptera mitogenome sequences (see the detailed species names in S1 Table). Additionally, we conducted mappings for the identified mitochondrial contigs using BWA v 0.7.13 [60], under default parameters. Mapping statistics were obtained to check the quality of the assembling with Qualimap [61] and Tablet [62], respectively.

Sequence alignment

All thirty-seven genes of insect mitogenome, including 13 protein-coding genes, 22 tRNA genes and two rRNA genes, were respectively aligned and used to further analyses. For protein-coding genes, firstly stop codons were excluded. Subsequently, each was aligned based on the invertebrate mitochondrial genetic code with Perl script TransAlign [63]. Each of tRNA and rRNA genes was aligned using MAFFT (version 7) under iterative refinement method incorporating the most accurate local (E-INS-i) pairwise alignment information [64]. The resultant alignments were checked in MEGA 6 [65]. Gaps were automatically striped by Gap Strip/Squeeze v2.1.0 with 40% Gap tolerance (http://www.hiv.lanl.gov/content/sequence/GAPSTREEZE/gap.html). All alignments were concatenated into the final data matrices with FASconCAT_v1.0 [66]. Two classes of data were compiled, namely with or without RNA genes.

To reduce the random similarity in the sequence alignments, each alignment of protein-coding genes, tRNA genes and rRNA genes was individually masked with Aliscore version 2.0 [67,68]. The sequence alignments were screened separately with the following options: -r 2000, -w 6, -N, -o outgroups. Positions identified as randomly similar by Aliscore were removed using Alicut version 2.0 [67,68]. After masking, the single gene alignments were concatenated with FASconCAT.

Potential saturation in the combined data sets was assessed using the index of substitution saturation (Iss) implemented in the DAMBE 5 [69,70]. Nucleotide homogeneity across taxa was assessed using the Chi-square test [71] implemented in PAUP*4.0b10 [72]. In addition, the results of IQ-TREE [73,74] also contain the Chi-square test, by which we can compare the composition of each sequence to the average composition across all data. The settings of IQ-TREE were the same as the following section of phylogenetic analysis.

Estimates of nonsynonymous (dN) and synonymous (dS) substitution rates of concatenated protein-coding genes were obtained by Yang and Nielsen (2000) [75] method using the program yn00 as implemented in PAML 4.9 [76]. The One-way Analysis of Variance (ANOVA) is performed in Excel 2016. In addition, the evolutionary rates of individual protein-coding genes were calculated using the same method mentioned above. Divergences of alignments of each protein-coding gene, rrnL, _rrnL_-SS, rrnS, _rrnS_-SS and tRNA were calculated using MEGA with the Maximum Composite Likelihood model.

Secondary structure prediction

The secondary structures of mitochondrial rrnS and rrnL were predicted mainly with reference to the method of “Comparative sequence analysis” by Ouvrard et al. (2000) [52]. Slight modifications were implemented in this study. Firstly, we used the MAFFT under the method of E-INS-i to produce the preliminary alignments, on which we could identify the conserved regions shared by the analyzed taxa and distinguish the hyper-variable portions present in each individual species. In the second step, the secondary structures of each conserved region and of each hyper-variable region were inferred, with reference to the models predicted for Apis mellifera [77] and/or Drosophila virilis [78]. Thirdly, putative helices were located through searching for uninterrupted base pairing in regions similar to the model of A. mellifera [77] or D. virilis [78]. Finally, Helix numbering followed the convention established at the CRW site [78]. The secondary structure diagrams of the complete mitochondrial rrnS and rrnL for eleven neuropterid insects are presented in Figs ​1 and ​2, Figures A-J in S1 Fig and Figures A-J in S2 Fig.

Secondary structure drawing of the complete mitochondrial rrnS gene for Chrysoperla nipponensis EMHAU-15090613.

The grey boxes indicated the helices aligned across all eleven sampled neuropterids.

Secondary structure drawing of the complete mitochondrial rrnL gene for Chrysoperla nipponensis EMHAU-15090613.

The grey boxes indicated the helices aligned across all eleven sampled neuropterids.

Based on secondary structures, sequences of rrnS and rrnL were realigned. Furthermore, both alignments were concatenated with subsets of protein-coding genes and tRNA genes, in order to investigate the effect of different sequence refinement algorithms on the phylogenetic reconstructions.

Phylogenetic analysis

Data partitioning

Prior to phylogenetic analyses, PartitionFinder [79] was employed to infer the optimal partitioning strategy, under a greedy search with RAxML [80]. For the data sets of PCG and PCGRNA, the data blocks were defined by gene types and by codon positions. For the remaining data sets, the data blocks were defined only by gene types. The partition schemes for each data set are provided in Tables A-N of S2 Table.

Phylogenetic reconstruction

The phylogenetic analyses were conducted using the following data sets: 1) PCG: 13 protein-coding genes; 2) PCG12: 13 protein-coding genes without the third codon positions; 3) PCG_AA: deduced amino acids of 13 protein-coding genes; 4) PCGRNA: combination of PCG, 2 rRNA and 22 tRNA genes; 5) PCG12RNA: combination of PCG12, 2 rRNA and 22 tRNA genes. Additionally, five corresponding data sets with masking were compiled: 1) Aliscore_PCG; 2) Aliscore_PCG12; 3) Aliscore_PCG_AA; 4) Aliscore_PCGRNA; 5) Aliscore_PCG12RNA.

To investigate the effect of mitochondrial rRNA realignments derived from secondary structures, eight corresponding data sets were utilized in phylogenetic analyses: 1) rrnS: rrnS aligned automatically by MAFFT; 2) rrnL: rrnL aligned automatically by MAFFT; 3) _rrnS-_SS: rrnS aligned manually based on secondary structures; 4) _rrnL-_SS: rrnL aligned manually based on secondary structures; 5) PCGRNA-SS: combining PCG, tRNA, and rRNA sequences realigned on secondary structures; 6) PCG12RNA-SS: combining PCG12, tRNA, and rRNA sequences realigned on secondary structures; 7) Aliscore_PCGRNA-SS: PCGRNA-SS with masking; 8) Aliscore_PCG12RNA-SS: PCG12RNA-SS with masking. In total, eighteen data sets were utilized in tree searches.

Maximum likelihood (ML) analyses were conducted using IQ-TREE [73,74], with various data partition schemes and the corresponding substitution models determined by PartitionFinder (Tables A-N in S2 Table). Branch support was estimated using Ultrafast bootstrap analysis with 1000 replicates.

Bayesian analyses were performed using a parallel version of PhyloBayes (pb_mpi1.5a) [81,82] as implemented on a HP server with twenty-four CPU and 320 G memory. The model CAT-GTR was used for nucleotide analyses, while the model CAT for amino acids. Two chains were run in parallel, and started from a random topology. The bpcomp program contained in the package of PhyloBayes was used to calculate the largest “maxdiff” and mean “meandiff” discrepancy observed across all bipartitions, and to yield the final consensus tree with the default options. The program tracecomp was also used to summarize the discrepancies and the effective sizes estimated for each column of the trace file. When the maximum “maxdiff” value was lower than 0.1 and minimum effective size was higher than 100, the Bayesian runs were recognized to be reached good convergence.

For the trees obtained, the bootstrap supports (BS) ≥ 75 and the posterior probabilities (PP) ≥ 0.90 were considered to be credible support values for the internal nodes. All sequence alignment files and tree files built in this article are available in Supporting Information S1–S3 Files.

Results

Mitogenomic assembly

Four complete neuropteran mitogenome sequences were obtained from each of single contigs. Mapping statistics showed that every site in each contig corresponded to the same base composition. Table 1 summarizes the statistics associated to the sequencing of the mitogenome assemblies. The total number of 125 bp reads ranged from 288,827,144 (M. angulatus) to 326,450,911 (C. nipponensis EMHAU-15090613). From these sequences, 0.22% to 0.27% corresponded to mitochondrial reads. The sequencing coverage ranged from 2,759× to 3,283×. Four mitogenomes had sequence lengths typical for insects, which ranged from 15,958 bp (Rapisma sp.) to 16,186 bp (Thaumatosmylus sp.). Each of them consisted of the full 37 mitochondrial genes and a putative control region. The mitogenomes showed the gene order proposed as ancestral for insects [31], except for the arrangement of tRNA cluster of trnC(gca)-trnW(tca)-trnY(gta). The variability of genome length occurred mainly in the control regions which ranged in size from 1,166 bp (Rapisma sp.) to 1,310 bp (M. angulatus). The organization of each mitogenome is presented in Tables A-D of S3 Table. All newly determined sequences have been deposited in GenBank (accession numbers: KX670539-KX670542).

Table 1

Statistics associated to the sequencing of mitogenomes using NGS-Illumina technology in four neuropteran species.

Sequence characteristics of data matrices

The saturation tests showed the transitions at the third codon positions to be saturated when assuming an asymmetrical topology (Iss > Iss.cAsym, Table 2). Nevertheless, all the remaining data partitions passed the test (Iss < Iss.cSym, and Iss < Iss.cAsym). Chi-square tests by PAUP based on the nucleotide data sets revealed significant heterogeneity among taxa (p < 0.05). To explore the source of base composition heterogeneity within each combined data set, the Chi-square tests were also performed by IQ-TREE. The results showed that several species contributed mainly to the sequence composition heterogeneity, for example, the Mongoloraphidia harmandi (Raphidioptera) and the three representatives of Ithonidae (Neuroptera) (Tables A and B in S4 Table). Comparing the number of taxa passing the Chi-square test across data sets indicated that the addition of RNA genes exacerbated the base composition difference among taxa. Removing the third codon positions, translating nucleotides into amino acids and applying alignment masking effectively reduced the composition heterogeneity.

Table 2

Saturation test.

For outgroup taxa, two coleopteran species had the higher values of dN (avg. 0.2290), dS (avg. 5.3490) and dN/dS (avg. 0.0428) (Table A in S5 Table). For the ingroup exemplars of Neuropterida, the single representative of Raphidioptera had a similar dN (0.2405), dS (5.2542) or dN/dS (0.0458) value as that of outgroup taxa. Whereas, the remaining ingroups displayed the lower rates of sequence evolution, with the average dN, dS and dN/dS values of 0.1640, 5.0798 and 0.0323 for Megaloptera, and 0.1560, 4.5043 and 0.0357 for Neuroptera. The one-way ANOVA analyses revealed incongruence in the dN or dN/dS values between groups (P < 0.05). However, there was no significant difference in the dS values (P = 0.4355).

Both the evolutionary rate analysis and the sequence divergence analysis showed the same trend for the protein-coding genes (Table B in S5 Table). The atp8 and nad6 genes have undergone accelerated evolution, as evidenced by the highest values of dN or dN/dS and the highest level of divergence. The cox1 and cytb genes are the relatively conserved protein-coding genes. Alignments based on the secondary structure information reduced the level of divergence of mitochondrial ribosomal RNA genes.

Secondary structures of mitochondrial rRNA genes

For the newly sequenced neuropteran species, the rrnS and rrnL were respectively identified between trnV (tac) and control region and between trnL1 (tag) and trnV (tac), with length ranging from 778 bp to 799 bp for rrnS, and from 1,314 bp to 1,324 bp for rrnL (Tables A-D in S3 Table). Both the location and the sequence lengths of rrnS and rrnL were consistent with other published neuropteran insects.

Fig 1 depicts the inferred secondary structure for the mitochondrial rrnS gene of C. nipponensis EMHAU-15090613, while the remaining rrnS secondary structures for ten representatives of Neuropterida are available in Figures A-J in S1 Fig. The alignment derived from these secondary structures spanned 784 nucleotide positions, which contained 402 variable positions.

The predicted secondary structures of the mitochondrial rrnS genes were similar to those of other arthropods [77,83,84], in which there were four typical domains identified (Fig 1, Figures A-J in S1 Fig and S6 Table). In the Domain I, Megaloptera and Neuroptera shared a similar secondary structure of H47, which was comprised of two long stems, one short stem and three loops (Fig 3B). However, the secondary structure of H47 predicted for M. harmandi (representing Raphidioptera) contained only one long stem, one short stem and two loops (Fig 3B). An alignment of H47 based on secondary structure information was performed for eleven neuropterid species (Fig 4A). The secondary structure of this motif can serve as a potential synapomorphy supporting the sister-group relationship between Megaloptera and Neuroptera. H17 was another highly conserved motif of Domain I, which included six identical nucleotides across eleven neuropterid rrnS sequences. The H885 in the Domain II and the H921, H939 and H1074 in the Domain III were relatively conserved, whose heterogeneity were mainly caused by the M. harmandi. Domain IV included two helices: H1399 and H1506, which were highly conserved in the sequence composition (conserved positions were 83.33% and 70.00%, respectively) (S6 Table). In particular, all neuropterids had an identical nucleotide composition for the H1506, with the exception of M. harmandi. Domain III and Domain IV were more conserved than Domain I and Domain II, with regard to the secondary structure and the base composition.

Phylogenetic reconstruction and secondary structure drawings of the helices.

(A) Phylogenetic reconstruction from PCG12RNA using software IQ-TREE under the partitioned models automatically selected. Node values represent bootstrap values (left) and posterior probabilities (right). The “-” indicates the relationship not being retrieved by the data set of PCG12RNA using software PhyloBayes under CAT-GTR model. Asterisks designate the species newly sequenced in this study. (B) Secondary structure drawings of the helix 47 in rrnS for eleven representative neuropterid species. (C) Secondary structure drawings of the helix 837 in rrnL for eleven representative neuropterid species.

Alignments of helices for eleven representative neuropterid species.

(A) Alignments of the helix 47 in rrnS. (B) Alignments of the helix 837 in rrnL. Helices are highlighted by grey box and numbered as diagrams of the full secondary structures. No color areas in each helix indicate the unpaired regions.

Secondary structure diagrams for the mitochondrial rrnL gene of C. nipponensis EMHAU-15090613 is shown in Fig 2, while the remaining rrnL secondary structures for ten representatives of Neuropterida are available in Figures A-J in S2 Fig. All predicted secondary structures of rrnL comprised of five canonical structural domains (I-II, IV-VI) and 44 helices (Fig 2, Figures A-J in S2 Fig and S7 Table), which were also found in other arthropod taxa [77,83–86]. The inferred rrnL secondary structures did not contain domain III, due to a shorten sequence region between domain II and IV (Fig 2 and Figures A-J in S2 Fig). The alignment of rrnL genes for eleven neuropterids based on the secondary structures spanned 1,250 nucleotide positions, with 679 variable positions. Domains I-II were more variable than the rest of domains (S7 Table). The secondary structures of H837 predicted for Megaloptera and Neuroptera were similar (Fig 3C). Whereas, the M. harmandi had a distinct structure of this helix, which contained an obviously longer unpaired region in the middle region (Fig 3C). Fig 4B illustrated the alignment derived from the secondary structure of H837 for eleven neuropterid species. The structure information of this helix may also be potential for the phylogenetic inference. Domains IV-VI were highly conserved, in which H1755, H1835, H1906, H1925, H2064, H2359, H2507 and H2735 had more than 75% sequence similarity, respectively.

Phylogenetic analyses

Phylogeny of Neuropterida

For the inter-order relationships within Neuropterida, the hypothesis of (Raphidioptera + (Megaloptera + Neuroptera)) was consistently recovered by all combined analyses (e.g. Fig 3A), albeit with strong supports for main nodes (BP > 90, PP > 0.90).

Across ML analyses, the differences occurred in the resolution of internal relationships within the Megaloptera and those within the Neuroptera. In the Megaloptera, two families were represented by more than two exemplars, namely the Sialidae and the Corydalidae. The monophyly of both families were strongly supported (BP = 100). At the genus level, the Neochauliodes was retrieved as non-monophyletic. Four species of Neochauliodes scattered among three clades in all ML analyses. Our sampling within the most diverse Neuroptera included 20 species representing eight families. Four families with more than one exemplars were congruently recovered as monophyletic: Osmylidae, Ithonidae, Hemerobiidae, and Chrysopidae. However, the Myrmeleontidae was retrieved as non-monophyletic with respect to Ascalaphidae in the ML analyses based on Aliscore_PCG and Aliscore_PCG-SS. All other ML analyses supported the Myrmeleontidae as a monophyletic group. In addition, the Ascalaphidae was recovered as monophyletic by several data sets (i.e. PCG, PCGRNA, PCG12RNA, Aliscore_PCG, Aliscore_PCGRNA and Aliscore_PCG12RNA). The remaining ML analyses recovered the Ascalaphidae to be a paraphyletic group, with Ascalohybris subjacens splitting as an independent branch. For the inter-familial relationships within Neuroptera, the Osmylidae was recovered as the most basal lineage. All other neuropteran taxa could be divided into two large groups. One included three families: Nymphidae, Myrmeleontidae, and Ascalaphidae. The other contained the Ithonidae, Mantispidae, Hemerobiidae, and Chrysopidae. The interrelationships within each grouping varied with different data sets. However, two sister-group relationships were often recovered: (Myrmeleontidae + Ascalaphidae) and (Hemerobiidae + Chrysopidae).

Bayesian analyses under the site-heterogeneous model produced the similar tree topologies to ML analyses. The monophyly of Neuropterida was strongly supported, and the sister-group relationship (Megaloptera + Neuroptera) was consistently recovered. Furthermore, the branching patterns within the Megaloptera and the Neuroptera were basically identical to those in ML trees. Bayesian and ML analyses differed only in the resolution of deep nodes. In the Bayesian analyses of Aliscore_PCGRNA, Aliscore_PCG12RNA, Aliscore_PCGRNA-SS and Aliscore_PCG12RNA-SS, the relationship among three neuropteridan orders was almost unresolved owing to the extremely short internal branch. At the family level, Bayesian analysis from Aliscore-PCG recovered a sister-group relationship between Hemerobiidae and Chrysopidae, while ML analysis from the same data set did not.

Effect of data treatment methods

Comparison between topologies from various data treatment methods showed that removing the third codon positions and translating nucleotides into amino acids had a marked influence on the ingroup relationships. In the analyses of PCG12 and PCG_AA, the monophyly of Ascalaphidae was not retrieved, while the sister group (Hemerobiidae + Chrysopidae) was always recovered. In contrast, sequence alignment masking (e.g. Aliscore_PCG) did not significantly alter the ingroup relationships under the same tree reconstruction method.

Realignment of mitochondrial rrnS gene sequences based on the secondary structure information significantly improved the resolution of relationships within Neuropterida (S3 Fig). However, the alignment of rrnL derived from the secondary structure produced an unexpected topology, in which the M. harmandi nested within the Megaloptera and was recovered as sister to the Sialidae. This may be due to the more hyper-variable regions contained in rrnL gene sequences, for which they are difficult to be aligned unambiguously even with reference to secondary structure information.

Discussion

NGS approach to mitogenome assembly

Mitogenome sequences have proved to be useful for studying Neuropterida phylogeny [6,18,25,28,30,87]. However, the availability of complete mitogenomes for neuropteridan insects is still limited. To date, there is no report on the NGS technology applied to the determination of neuropterid mitogenome. In this study, we employed high-throughput pooled sequencing approach to mitogenome assembling, and successfully reconstructed four complete mitogenomes of Neuroptera. Combined with the existing neuropteran mitogenome sequences, the newly sequenced mitogenomes were utilized for phylogenetic reconstruction of Neuropterida. The tree building results presented congruent phylogenetic hypotheses with recent studies on the higher-level phylogeny of Neuropteria [7,16,20,23,28,29,30], in which all four newly determined mitogenome sequences were consistently positioned in the clade including their close relatives. In addition, the gene arrangement, nucleotide composition, gene length, anti-codons of tRNA, start and stop codon usage of protein-coding genes and sequence evolutionary rate were similar to most of neuropteran mitogenomes published. All these results validated the applicability of acquiring large number of neuropterid mitogenomes from de novo assembly of high-throughput pooled sequencing. Furthermore, more insect species can be added into the pool to improve sequencing effectiveness from the point of view of deeper sequencing coverage applied. At present, sequencing cost has been decreased to < US$300 for constructing a single library and generating 10 G data by a standard NGS platform. Thus, the NGS approach for reconstructing mitogenomes will become popular for insect mitochondrial phylogenomics due to its relatively easy protocol and remarkably lower cost.

Genome organization

Just when we finished writing this paper, another research of mitogenomic phylogeny of Neuropterida had been published [30]. In the analysis of Wang et al. (2017) [30], a greater number of mitogenomes of represented Neuropterida were sequenced by PCR amplification and traditional Sanger sequencing method. Wang et al. (2017) [30] revealed the tRNA arrangement of trnW(tca)-trnC(gca)-trnY(gta) (i.e. WCY) as an ancestral gene organization harbored by the neuropteran families Coniopterygidae, Nevrorthidae, Osmylidae and Sisyridae. Furthermore, they suggested the tRNA arrangement of WCY as the potential evidence for the basal position of Osmylidae in Neuroptera. However, our sequencing result demonstrated that the tRNA arrangement of trnC(gca)-trnW(tca)-trnY(gta) (i.e. CWY in Thaumatosmylus sp.) also occurred in the Osmylidae. Therefore, whether the arrangement of this tRNA cluster can be regarded as an evidence for neuropterid phylogeny requires additional data.

Comparative analysis of mitochondrial rRNA secondary structures

By comparing mitochondrial rRNA secondary structures between species, conserved motifs and highly variable portions within each canonical structural domain were identified (Figs ​1 and ​2, Figures A-J in S1 Fig, Figures A-J in S2 Fig, and S6 and S7 Tables). Subsequently, the conserved sequence motifs and their associated secondary structure elements were utilized to identify the homologous positions in sequences. In this paper, we analyzed the complete mitochondrial rrnS and rrnL gene sequences for eleven neuropterid species. The results demonstrated that realigned rRNA data could improve the resolution of relationships within Neuropterida (S3 Fig). Two conserved motifs (i.e. H47 for rrnS and H837 for rrnL) were found to be the potential evidence for a close relation between Megaloptera and Neuroptera (Fig 3B and 3C). In the prior studies, only one neuropteran insect mitochondrial rrnS and rrnL secondary structures were predicted, namely the Libelloides macaronius [87]. Some regions inferred for the secondary structures of L. macaronius mitochondrial rrnS and rrnL were different from other insects published [77,78]. For example, the Domain I proposed for the secondary structure of L. macaronius mitochondrial rrnL was distinct from the commonly utilized model of A. mellifera [77]. The prediction for the structure of this region by Negrisolo et al. (2011) [87] was not corroborated by the current study, either. The eleven neuropterid rrnL sequences can be inferred as the A. mellifera model [77]. With a denser taxon sampling, the diagrams of mitochondrial rrnS and rrnL secondary structures presented in this study can be useful to improve the accuracy of the neuropterid secondary structures, and allow better visualization and comparison of mitochondrial rRNA structural features. These will be of fundamental importance to their application in phylogenetic studies and to understanding the function of these molecules for neuropterid insects.

Neuropterida phylogeny

With respect to the higher-level phylogeny of Neuropterida, the sister-group relationship between Raphidioptera and Megaloptera has long been recognized [8,9,88] by some of Neuropterida systematists. This view was also advocated in the study of Beutel et al. (2010) [21]. Synapomorphies supporting this relation were mainly from the morphological structures of larval and adult head. But these characters were re-considered as symplesiomorphic by Aspöck and Haring (2012) [19]. In the trees recovered in this study, only two analyses based on the single rRNA gene alignments (i.e. _rrnL_-SS and rrnS) under ML inference recovered a close relationship between Raphidioptera and Megaloptera, but with weak nodal support values (BP < 75). This result was likely to be an artifact of inaccurate alignments of mitochondrial rRNA genes. Because the rrnL gene sequences contained more highly variable regions, while the automatically aligned rrnS gene sequences had insufficient phylogenetic information on the resolution of deep nodes in the neuropterid tree. In contrast, utilizing secondary structure information resulted in a marked increase in the phylogenetic signal and the lower sequence divergence of alignment of _rrnS_-SS. The conserved sites, variable sites and parsimony-information sites relative to all sites for _rrnS_-SS and rrnS are 254/548/433/802 and 282/492/401/774, respectively. Thus, a congruent topology with the well-supported tree from the combined data can be recovered by _rrnS_-SS alignment.

The majority of trees from mitogenomic data provided strong support for the hypothesis of (Raphidioptera + (Megaloptera + Neuroptera)). In addition, we provided a potential evidence from the secondary structures of two conserved motifs (i.e. H47 in rrnS and H837 in rrnL) for this inter-ordinal relationship. This branching pattern was also supported by recent morphological [23,29] and molecular [7,16,28,30] studies. Supporting evidence mainly included the larval morphological characters [23], the wing base structures [29] and the mitogenomic data [30]. In the study of Wang et al. (2017) [30], significant compositional heterogeneity of mitogenome sequences across lineages had an adverse effect on the resolution of higher-level phylogeny of Neuropterida, especially under homogeneous model and parsimony analyses. More recently, Winterton et al. (2017) [20] provided the most comprehensive study of Neuropterida phylogeny using anchored phylogenomics to date. And they recovered the monophyletic Raphidioptera as sister to the rest of Neuropterida in all ML and Bayesian analyses, with strong statistical support [20]. However, they indicated that sequence heterogeneity had no apparent negative impact on phylogeny estimation. Our sequence characteristic analysis also showed a rate heterogeneity across miogenome sequences of three neuropterid orders. The monophyly of each order and the similar inter-ordinal relationships as recent studies [20,30] were retrieved in all analyses, regardless of models or algorithms utilized. The difference at the lower-level relationships between studies may be due to the various taxon sampling. Because the effect of sequence heterogeneity is likely to be exacerbated with expanding sequence data from limited species coverage. We should acknowledge that all current studies included still sparse taxon sampling for a highly diverse Neuropterida, which might have an influence on the resulting phylogeny. Addition of more mitogenome or the whole genome-wide data, especially those for species of Raphidioptera, to future neuropteridan phylogenetic studies can provide new insights into the evolutionary relationships among three orders. Thus, we considered the current phylogenetic result on the ordinal relationships of Neuropterida to be still tentative.

The monophyly of Megaloptera was questioned by several previous studies [19,20]. However, the results of these studies were considered as incredible by Aspöck and Haring (2012) [19], due to the various tree searching algorithms, character or taxon sampling limitation. In this paper, the sister-group relationship between Sialidae and Corydalidae were consistently recovered by the present mitogenomic data, with strong nodal support (BP = 100, PP = 1). This relationship was also supported by recent morphological [5,89] and molecular studies [4,30]. Therefore, our results confirmed a monophyletic Megaloptera, encompassing two families of Sialidae and Corydalidae.

Within Neuroptera, eight families were included in the present study. The Osmylidae were recovered with strong support (BP = 100, PP = 1) as the earliest offshoot to all other neuropteran families. The basal position of Osmylidae was largely congruent with previous molecular study based on the single nuclear and mitochondrial gene fragments [4] and the mitogenomic data [30]. The clade (Nymphidae + (Myrmeleontidae + Ascalaphidae)) was recovered by most of combined analyses. The same result was also retrieved by some prior studies [4,5,19,22,23,90,91,92,93]. However, Michel et al (2017) [94] showed that Nemopteridae had a closer relationship to the clade (Myrmeleontidae + Ascalaphidae) than Nymphidae. Winterton et al. (2017) [20] also placed the Nemopteridae as sister to the assemblage including Myrmeleontidae and Ascalaphidae. Therefore, the clade comprising (Nymphidae + (Myrmeleontidae + Ascalaphidae)) recovered by the present data may be due to the absence of Nemopteridae and Psychopsidae. More extensive sampling is needed to definitively elucidate relationships between these families. In the study by Winterton et al. (2017) [20], the nested position of Ascalaphidae led to a paraphyletic Myrmeleontidae. Similarly, the paraphyly of Myrmeleontidae was also proposed by Wang et al. (2017) [30]. However, most of our analyses retrieved the Ascalaphidae and Myrmeleontidae as diverging clades. The cluster of the remaining four families recovered in our analyses had the same constitution as those in the study by Haring and Aspöck (2004) [4], who recognized the group (Ithonidae + Polystoechotidae + Hemerobiidae + Chrysopidae + Mantispidae) according to A + T content of nuclear 18S rRNA sequences. Of which, the Polystoechotidae is more often considered as a member of the family Ithonidae [20,30]. Winterton and Makarkin (2010) [95] actually demonstrated the clade comprising Ithonidae and Polystoechotidae as a monophyletic group (i.e. an unique family) based on molecular, morphological and fossil evidence. In this study, we argued for merging the Polystoechotidae into the family Ithonidae. The limited taxon sampling led to an aberrant hierarchy arrangement of four families in several analyses (e.g. the topology depicted in Fig 3A). This branch pattern is different from the most recent analyses. In particular, the Chrysopidae was retrieved in a terminal position, which contrasted with a deeper placement within Neuroptera in two recent studies [20,30]. Nevertheless, the sister-group (Hemerobiidae + Chrysopidae) was frequently recovered in our analyses (e.g. the ML and Bayesian analyses of PCG, PCG12 and PCG_AA). This result agreed with the previous molecular studies by Haring and Aspöck (2004) [4] and Wang et al. (2017) [30], but disagreed with that of Winterton et al. (2017) [20].

Supporting information

S1 Fig

Secondary structure drawing of the complete mitochondrial rrnS gene for ten neuropterid species.

(A) Acanthacorydalis orientalis, (B) Ascalohybris subjacens, (C) Corydalus cornutus, (D) Dysmicohermes ingens, (E) Micromus angulatus, (F) Mongoloraphidia harmandi, (G) Neochauliodes fraternus, (H) Thaumatosmylus sp., (I) Rapisma sp., and (J) Sialis hamate.

(ZIP)

S2 Fig

Secondary structure drawing of the complete mitochondrial rrnL gene for ten neuropterid species.

(A) Acanthacorydalis orientalis, (B) Ascalohybris subjacens, (C) Corydalus cornutus, (D) Dysmicohermes ingens, (E) Micromus angulatus, (F) Mongoloraphidia harmandi, (G) Neochauliodes fraternus, (H) Thaumatosmylus sp., (I) Rapisma sp., and (J) Sialis hamate.

(RAR)

S3 Fig

Phylogenetic reconstruction from the data set of _rrnS-_SS using software IQ-TREE under GTR+G+I model.

Node values represent bootstrap values.

(TIF)

S1 Table

Taxa used in this study.

(XLSX)

S2 Table

The partition schemes and best-fitting models selected by PartitionFinder for each dataset.

(XLSX)

S3 Table

The organization of four newly sequenced mitogenomes.

(XLSX)

S4 Table

The Chi-square tests.

(XLSX)

S5 Table

The evolutionary rate analyses conducted for each species and for each gene partition.

(XLSX)

S6 Table

The nucleotide composition in each helix of the secondary structures of rrnS in 11 Neuropterida mitogenomes.

(XLSX)

S7 Table

The nucleotide composition in each helix of the secondary structures of rrnL in 11 Neuropterida mitogenomes.

(XLSX)

S1 File

All sequence alignment files in PHYLIP format compiled in this study.

(ZIP)

S2 File

All ML trees in NEXUS format reconstructed in this study.

(ZIP)

S3 File

All Bayesian trees in NEXUS format reconstructed in this study.

(ZIP)

Funding Statement

This research is supported by grants from the National Natural Science Foundation of China (no. 31402002), Key Scientific Research Projects of Henan Province (grant no. 16A210029), and Henan Academician Workstation of Pest Green Prevention and Control for Plants in Southern Henan (YZ201601). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability

All relevant data are within the paper and its Supporting Information files.

References

1. Senior LJ, McEwen P. The use of lacewings in biological control In: McEwen P, New TR & Whittington AE (Eds.). Lacewings in the Crop Environment; 2001. pp. 296–302. [Google Scholar]

2. Daane KM, Hagen KS. An evalution of lacewing release in North America In: McEwen P, New TR & Whittington AE (Eds.). Lacewings in the Crop Environment; 2001. pp. 398–407. [Google Scholar]

4. Haring E, Aspöck U. Phylogeny of the Neuropterida: a first molecular approach. Syst Entomol. 2004;29: 415–430. [Google Scholar]

5. Aspöck U, Plant JD, Nemeschkal HL. Cladistic analysis of Neuroptera and their systematic position within Neuropterida (Insecta: Holometabola: Neuropterida: Neuroptera). Syst Entomol. 2001;26: 73–86. [Google Scholar]

6. Aspöck U. Phylogeny of the Neuropterida (Insecta: Holometabola). Zool Scr. 2002;31: 51–55. [Google Scholar]

7. Cameron SL, Sullivan J, Song H, Miller KB, Whiting MF. A mitochondrial genome phylogeny of the Neuropterida (lace-wings, alderflies and snakeflies) and their relationship to the other holometabolous insect orders. Zool Scr. 2009;38: 575–590. [Google Scholar]

8. Kristensen NP. Phylogeny of extant hexapodsVolume I and IISecond Edition Naumann I. D., Carne P. B., Lawrence J. F., Nielsen E. S., Spradberry J. P., Taylor R. W., Whitten M. J. and Littlejohn M. J. eds. Carlton, Victoria: Melbourne University Press; Insects of Australia: A Textbook for Students and Research Workers; 1991. pp. 125–140. [Google Scholar]

9. Kristensen NP. Phylogeny of endopterygote insects, the most successful lineage of living organisms. Eur J Entomol. 1999;96: 237–253. [Google Scholar]

10. Kjer KM, Carle FL, Litman J, Ware J. A molecular phylogeny of Hexapoda. Arthropod Syst Phylo. 2006;64: 35–44. [Google Scholar]

11. Boudreaux HB. Arthropod phylogeny with special reference to insects. John Wiley & Sons, New York, Chichester, Brisbane, Toronto: 1979. [Google Scholar]

12. Whiting MF. Phylogeny of the holometabolous insect orders: molecular evidence. Zool Scr. 2002;31: 3–15. [Google Scholar]

13. Caterino MS, Shull VL, Hammond PM, Vogler AP. Basal relationships of Coleoptera inferred from 18S rDNA sequences. Zool Scr. 2002;31: 41–49. [Google Scholar]

14. Misof B, Niehuis O, Bischoff I, Rickert A, Erpenbeck D, Staniczek A. Towards and 18S phylogeny of hexapods: Accounting for group-specific character covariance in optimized mixed nucleotide/doublet models. Zoology. 2007;110: 409–429. doi: 10.1016/j.zool.2007.08.003 [PubMed] [Google Scholar]

15. Wiegmann BM, Trautwein MD, Kim JW, Cassel BK, Bertone MA, Winterton SL, et al. Single copy nuclear genes resolve the phylogeny of the holometabolous insects. BMC Biol. 2009;7: 34 doi: 10.1186/1741-7007-7-34 [PMC free article] [PubMed] [Google Scholar]

16. Misof B, Liu S, Meusemann K, Peters RS, Donath A, Mayer C, et al. Phylogenomics resolves the timing and pattern of insect evolution. Science. 2014;346: 763–767. doi: 10.1126/science.1257570 [PubMed] [Google Scholar]

17. Song N, An SH, Yin XM, Zhao T, Wang XY. Insufficient resolving power of mitogenome data in deciphering deep phylogeny of Holometabola. J Syst Evol. 2016;54: 545–559. [Google Scholar]

18. Song F, Li H, Jiang P, Zhou X, Liu J, Sun C, et al. Capturing the phylogeny of Holometabola with mitochondrial genome data and Bayesian site-heterogeneous mixture models. Genome Biol Evol. 2016;8: 1411–1426. doi: 10.1093/gbe/evw086 [PMC free article] [PubMed] [Google Scholar]

19. Aspöck U, Haring E. The phylogeny of the Neuropterida: long lasting and current controversies and challenges (Insecta: Endopterygota). Arthropod Syst Phy. 2012;70: 119–129. [Google Scholar]

20. Winterton SL, Lemmon AR, Gillung JP, Garzon IJ, Badano D, et al. Evolution of lacewings and allied orders using anchored phylogenomics (Neuroptera, Megaloptera, Raphidioptera). Syst. Entomol. 2017. doi: 10.1111/syen.12278 [Google Scholar]

21. Beutel RG, Friedrich F, Hörnschemeyer T, Pohl H, Hünefeld F, Beckmann F, et al. Morphological and molecular evidence converge upon a robust phylogeny of the megadiverse Holometabola. Cladistics. 2010;26: 1–15. [PubMed] [Google Scholar]

22. Winterton SL, Hardy NB, Wiegmann BM. On wings of lace: phylogeny and bayesian divergence time estimates of neuropterida (insecta) based on morphological and molecular data. Syst Entomol. 2010;35: 349–378. [Google Scholar]

23. Aspöck U, Aspöck H. Phylogenetic relevance of the genital sclerites of Neuropterida (Insecta: Holometabola). Syst Entomol. 2008;33: 97–127. [Google Scholar]

24. Yang D, Liu XY. Fauna Sinica InsectaVol 51 Megaloptera. Beijing in Chinese: Science Press; 2010. [Google Scholar]

25. Zhao J, Li H, Winterton SL, Liu Z. Ancestral gene organization in the mitochondrial genome of Thyridosmylus langii (McLachlan, 1870) (Neuroptera: Osmylidae) and implications for lacewing evolution. PLoS One. 2013;8: e62943 doi: 10.1371/journal.pone.0062943 [PMC free article] [PubMed] [Google Scholar]

26. Whiting MF. Phylogeny of the holometabolous insect orders: molecular evidence. Zool Scr. 2002;31: 3–15. [Google Scholar]

27. Wheeler WC, Whiting MF, Wheeler QD, Carpenter JM. The phylogeny of the extant hexapod orders. Cladistics. 2001;17: 113–169. [PubMed] [Google Scholar]

28. Wang YY, Liu XY, Winterton SL, Yang D. The first mitochondrial genome for the fishfly subfamily Chauliodinae and implications for the higher phylogeny of Megaloptera. PLoS One. 2012;7: e47302 doi: 10.1371/journal.pone.0047302 [PMC free article] [PubMed] [Google Scholar]

30. Wang Y, Liu X, Garzón-Orduña IJ, Winterton SL, Yan Y, Aspöck U, et al. Mitochondrial phylogenomics illuminates the evolutionary history of Neuropterida. Cladistics. 2017;33: 617–636. [PubMed] [Google Scholar]

31. Cameron SL. How to sequence and annotate insect mitochondrial genomes for systematic and comparative genomics research. Syst Entomol. 2014;39: 400–411. [Google Scholar]

32. Timmermans MJ, Dodsworth S, Culverwell CL, Bocak L, Ahrens D, Littlewood DT, et al. Why barcode? High-throughput multiplex sequencing of mitochondrial genomes for molecular systematics. Nucleic Acids Res. 2010;38: e197 doi: 10.1093/nar/gkq807 [PMC free article] [PubMed] [Google Scholar]

33. Coates BS. Assembly and annotation of full mitochondrial genomes for the corn rootworm species, Diabrotica virgifera and Diabrotica barberi (Insecta: Coleoptera: Chrysomelidae), using Next Generation Sequence data. Gene. 2014;542: 190–197. doi: 10.1016/j.gene.2014.03.035 [PubMed] [Google Scholar]

34. Gillett CP, Crampton-Platt A, Timmermans MJ, Jordal BH, Emerson BC, Vogler AP, et al. Bulk de novo mitogenome assembly from pooled total DNA elucidates the phylogeny of weevils (Coleoptera: Curculionoidea). Mol Biol Evol. 2014;31: 2223–2237. doi: 10.1093/molbev/msu154 [PMC free article] [PubMed] [Google Scholar]

35. Crampton-Platt A, Timmermans MJ, Gimmel ML, Kutty SN, Cockerill TD, Vun Khen C, et al. Soup to tree: the phylogeny of beetles inferred by mitochondrial metagenomics of a bornean rainforest sample. Mol Biol Evol. 2015;32: 2302–2316. doi: 10.1093/molbev/msv111 [PMC free article] [PubMed] [Google Scholar]

36. Ramakodi MP, Singh B, Wells JD, Guerrero F, Ray DA. A 454 sequencing approach to dipteran mitochondrial genome research. Genomics. 2015;105: 53–60. doi: 10.1016/j.ygeno.2014.10.014 [PubMed] [Google Scholar]

37. Dalevi D, Hugenholtz P, Blackall LL. A multiple-outgroup approach to resolving division-level phylogenetic relationships using 16S rDNA data. Int J Syst Evol Micr. 2001;51: 385–391. [PubMed] [Google Scholar]

38. Lee IM, Hammond RW, Davis RE, Gundersen DE. Universal amplification and analysis of pathogen 16S rDNA for classification and identification of mycoplasmalike organisms. Phytopathology. 1993;83: 834–842. [Google Scholar]

39. Black WC, Piesman J. Phylogeny of hard- and soft-tick taxa (Acari: Ixodida) based on mitochondrial 16S rDNA sequences. P Natl Acad Sci Usa. 1994;91: 10034–10038. [PMC free article] [PubMed] [Google Scholar]

40. Hasegawa E, Kasuya E. Phylogenetic analysis of the insect order odonata using 28S and 16S rDNA sequences: a comparison between data sets with different evolutionary rates. Entomol Sci. 2006;9: 55–66. [Google Scholar]

41. Yang Q, Zhao S, Kucerová Z, Stejskal V, Opit G, Qin M, et al. Validation of the 16S rDNA and COI DNA barcoding technique for rapid molecular identification of stored product psocids (Insecta: Psocodea: Liposcelididae). J Econ Entomol. 2013;106: 419–425. [PubMed] [Google Scholar]

42. Halanych KM, Robinson TJ. Multiple substitutions affect the phylogenetic utility of Cytochrome b, and 12S rDNA data: examining a rapid radiation in Leporid (Lagomorpha) Evolution. J Mol Evol. 1999;48: 369–379. [PubMed] [Google Scholar]

43. Colgan DJ, Costa PD. Megachiropteran evolution studied with 12S rDNA and c-mos dna sequences. J Mamm Evol. 2002;9: 3–22. [Google Scholar]

44. Murrell A, Campbell NJ, Barker SC. Mitochondrial 12S rDNA indicates that the rhipicephalinae (acari: ixodida) is paraphyletic. Mol Phylogenet Evo. 1999;12: 83–86. [PubMed] [Google Scholar]

45. Møller PR, Gravlund P. Phylogeny of the eelpout genus Lycodes (Pisces, Zoarcidae) as inferred from mitochondrial cytochrome b and 12S rDNA. Mol Phylogenet Evo. 2003;26: 369–388. [PubMed] [Google Scholar]

46. Kjer KM. Use of rrna secondary structure in phylogenetic studies to identify homologous positions: an example of alignment and data presentation from the frogs. Mol Phylogenet Evo. 1995;4: 314–330. [PubMed] [Google Scholar]

47. Titus TA, Frost DR. Molecular homology assessment and phylogeny in the lizard family Opluridae (Squamata: Iguania). Mol Phylogenet Evol. 1996;6: 49–62. doi: 10.1006/mpev.1996.0057 [PubMed] [Google Scholar]

48. Buckley TR, Simon C, Flook PK, Misof B. Secondary structure and conserved motifs of the frequently sequenced domains iv and v of the insect mitochondrial large subunit rrna gene. Insect Mol Biol. 2000;9: 565–580. [PubMed] [Google Scholar]

49. Xie Q, Tian Y, Zheng L, Bu W. 18S rRNA hyper-elongation and the phylogeny of euhemiptera (Insecta: Hemiptera). Mol Phylogenet Evol. 2008;47: 463–471. doi: 10.1016/j.ympev.2008.01.024 [PubMed] [Google Scholar]

50. Tsagkogeorga G, Turon X, Hopcroft RR, Tilak MK, Feldstein T, Shenkar N, et al. An updated 18S rRNA phylogeny of tunicates based on mixture and secondary structure models. BMC Evol Biol. 2009;9: 187 doi: 10.1186/1471-2148-9-187 [PMC free article] [PubMed] [Google Scholar]

51. De RP, De WR. DCSE, an interactive tool for sequence alignment and secondary structure research. Bioinformatics. 1993;9: 735–740. [PubMed] [Google Scholar]

52. Ouvrard D, Campbell BC, Bourgoin T, Chan KL. 18s rRNA secondary structure and phylogenetic position of Peloridiidae (Insecta: Hemiptera). Mol Phylogenet Evol. 2000;16: 403–417. doi: 10.1006/mpev.2000.0797 [PubMed] [Google Scholar]

53. Caetano-Anolle´s G. Tracing the evolution of RNA structure in ribosomes. Nucleic Acids Res. 2002;30: 2575–2587. [PMC free article] [PubMed] [Google Scholar]

54. Andrews S. FastQC, A quality control tool for high throughput sequence data. Bioinformatics Group, Babraham Institute, Cambridge, United Kingdom. 2010. (www.bioinformatics.babraham.ac.uk/projects/fastqc).

56. Peng Y, Leung HC, Yiu SM, Chin FY. IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics. 2012;28: 1420–1428. doi: 10.1093/bioinformatics/bts174 [PubMed] [Google Scholar]

58. Hall TA. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp Ser. 1999;41: 95–98. [Google Scholar]

59. Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013;69: 313–319. doi: 10.1016/j.ympev.2012.08.023 [PubMed] [Google Scholar]

61. García-Alcalde F, Okonechnikov K, Carbonell J, Cruz LM, Götz S, Tarazona S, et al. Qualimap: evaluating next-generation sequencing alignment data. Bioinformatics. 2012;28: 2678–2679. doi: 10.1093/bioinformatics/bts503 [PubMed] [Google Scholar]

62. Milne I, Stephen G, Bayer M, Cock PJ, Pritchard L, Cardle L, et al. Using Tablet for visual exploration of second-generation sequencing data. Briefings in Bioinformatics. 2013;14: 193–202. doi: 10.1093/bib/bbs012 [PubMed] [Google Scholar]

67. Misof B, Misof K. A Monte Carlo approach successfully identifies randomness in multiple sequence alignments: a more objective means of data exclusion. Syst Biol. 2009;58: 21–34. doi: 10.1093/sysbio/syp006 [PubMed] [Google Scholar]

68. Kück P, Meusemann K, Dambach J, Thormann B, von Reumont BM, Wägele JW, et al. Parametric and non-parametric masking of randomness in sequence alignments can be improved and leads to better resolved trees. Front Zool. 2010;7: 1–12. [PMC free article] [PubMed] [Google Scholar]

69. Xia X, Xie Z, Salemi M, Chen L, Wang Y. An index of substitution saturation and its application. Mol Phylogenet Evol. 2003;26: 1–7. [PubMed] [Google Scholar]

71. Farris JS, Källersjö M, Kluge AG, Bult C. Testing significance of incongruence. Cladistics1994;10: 315–319. [Google Scholar]

72. Swofford DL. PAUP: phylogenetic analysis using parsimony and other methods. Version 4 beta 10 Sinauer Associates, Sunderland, Massachusetts: 2002. [Google Scholar]

73. Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Mol Biol Evol. 2015;32: 268–274. doi: 10.1093/molbev/msu300 [PMC free article] [PubMed] [Google Scholar]

74. Trifinopoulos J, Nguyen LT, Von HA, Minh BQ. W-iq-tree: a fast online phylogenetic tool for maximum likelihood analysis. Nucleic Acids Res. 2016;44: W232–W235. doi: 10.1093/nar/gkw256 [PMC free article] [PubMed] [Google Scholar]

77. Gillespie JJ, Johnston JS, Cannone JJ, Gutell RR. Characteristics of the nuclear (18S, 5.8S, 28S and 5S) and mitochondrial (12S and 16S) rRNA genes of Apis mellifera (Insecta: Hymenoptera): structure, organization, and retrotransposable elements. Insect Mol Biol. 2006;15: 657–686. doi: 10.1111/j.1365-2583.2006.00689.x [PMC free article] [PubMed] [Google Scholar]

78. Cannone JJ, Subramanian S, Schnare MN, Collett JR, D’Souza LM, Du Y, et al. The comparative RNA web (CRW) site: An online database of comparative sequence and structure information for ribosomal, intron, and other RNAs. BMC Bioinformatics. 2002;3: 15. [PMC free article] [PubMed] [Google Scholar]

79. Lanfear R, Calcott B, Ho SY, Guindon S. Partitionfinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29: 1695–1701. doi: 10.1093/molbev/mss020 [PubMed] [Google Scholar]

80. Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22: 2688–2690. doi: 10.1093/bioinformatics/btl446 [PubMed] [Google Scholar]

81. Lartillot N, Philippe HABayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Mol Biol Evol. 2004;21: 1095–1109. doi: 10.1093/molbev/msh112 [PubMed] [Google Scholar]

82. Lartillot N, Lepage T, Blanquart S. PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating. Bioinformatics. 2009;25: 2286–2288. doi: 10.1093/bioinformatics/btp368 [PubMed] [Google Scholar]

83. Lavrov DV, Boore JL, Brown WM. The complete mitochondrial DNA sequence of the horseshoe crab Limulus polyphemus. Mol Evol. 2000;17: 813e24. [PubMed] [Google Scholar]

84. Jeyaprakash A, Hoy MA. The mitochondrial genome of the predatory mite Metaseiulus occidentalis (Arthropoda: Chelicerata: Acari: Phytoseiidae) is unexpectedly large and contains several novel features. Gene. 2007;391: 264e74. [PubMed] [Google Scholar]

85. Gutell RR, Fox GE. A compilation of large subunit RNA sequences presented in structural format. Nucleic Acids Res. 1988;16: R175e313. [PMC free article] [PubMed] [Google Scholar]

86. Misof B, Fleck G. Comparative analysis of mt LSU rRNA secondary structures of Odonates: structural variability and phylogenetic signal. Insect Mol Biol. 2003;12: 535e47. [PubMed] [Google Scholar]

87. Negrisolo E, Babbucci M, Patarnello T. The mitochondrial genome of the ascalaphid owlfly Libelloides macaronius and comparative evolutionary mitochondriomics of neuropterid insects. BMC Genomics. 2011;12: 221 doi: 10.1186/1471-2164-12-221 [PMC free article] [PubMed] [Google Scholar]

88. Aspöck H, Aspöck U, Rausch H. Die Raphidiopteren der Erde. Eine monographische Darstellung der Systematik, Taxonomie, Biologie, Ökologie und Chorologie der rezenten Raphidiopteren der Erde, mit einer Zusammenfassenden UÈ bersicht der Fossilen Raphidiopteren (Insecta: Neuropteroidea). Vol 2. Goecke & Evers, Krefeld, Germany. 1991.

89. Liu X, Yanan LÜ, Aspöck H, Yang D, Aspöck U. Homology of the genital sclerites of megaloptera (insecta: neuropterida) and their phylogenetic relevance. Syst Entomol. 2016;41: 256–286. [Google Scholar]

90. Stange LA. Reclassification of the New World antlion genera formerly included in the tribe Brachynemurini. Insecta Mundi. 1994;8: 67–119. [Google Scholar]

91. Kuznetsova VG, Khabiev GN, Krivokhatsky VA. Chromosome numbers in antlions (Myrmeleontidae) and owlflies (Ascalaphidae) (Insecta, Neuroptera). Zookeys. 2014;538: 47–61. [PMC free article] [PubMed] [Google Scholar]

92. Badano D, Aspöck U, Aspöck H, Cerretti P. Phylogeny of Myrmeleontiformia based on larval morphology (Neuropterida: Neuroptera). Syst Entomol. 2017;42: 94–117. [Google Scholar]

93. Zhang L, Yang J. The mitochondrial genome of Gatzara jezoensis (Neuroptera: Myrmeleontidae) and phylogenetic analysis of Neuroptera. Biochem Syst Ecol. 2017;71: 230–235. [Google Scholar]

94. Michel B, Clamens AL, Béthoux O, Kergoat GJ, Condamine FL. A first higher-level time-calibrated phylogeny of antlions (Neuroptera: Myrmeleontidae). Mol Phylogenet Evol. 2017;107: 103–116. doi: 10.1016/j.ympev.2016.10.014 [PubMed] [Google Scholar]

95. Winterton SL, Makarkin VN. Phylogeny of moth lacewings and giant lacewings (Neuroptera: Ithonidae, Polystoechotidae) using DNA sequence data, morphology, and fossils. Ann Entomol Soc Am. 2010;103: 511–522. [Google Scholar]

Articles from PLOS ONE are provided here courtesy of PLOS