The first complete mitochondrial genome of marigold pest thrips, Neohydatothrips samayunkur (Sericothripinae) and comparative analysis

Complete mitogenomes from the order Thysanoptera are limited to representatives of the subfamily Thripinae. Therefore, in the present study, we sequenced the mitochondrial genome of Neohydatothrips samayunkur (15,295 bp), a member of subfamily Sericothripinae. The genome possesses the canonical 13 protein-coding genes (PCGs), 22 transfer RNA genes (tRNAs), and two ribosomal RNA genes (rRNAs) as well as two putative control regions (CRs). The majority strand was 77.42% A + T content, and 22.58% G + C with weakly positive AT skew (0.04) and negative GC skew (−0.03). The majority of PCGs start with ATN codons as observed in other insect mitochondrial genomes. The GCG codon (Alanine) was not used in N. samayunkur. Most tRNAs have the typical cloverleaf secondary structure, however the DHU stem and loop were absent in trnV and trnS1, while the TΨC loop was absent in trnR and trnT. The two putative control regions (CR1 and CR2) show 99% sequence similarity indicated a possible duplication, and shared 57 bp repeats were identified. N. samayunkur showed extensive gene rearrangements, with 11 PCGs, 22 tRNAs, and two rRNAs translocated when compared to the ancestral insect. The gene trnL2 was separated from the ‘trnL2-cox2’ gene block, which is a conserved, ancestral gene order found in all previously sequenced thrips mitogenomes. Both maximum likelihood (ML) and Bayesian inference (BI) phylogenetic trees resulted in similar topologies. The phylogenetic position of N. samayunkur indicates that subfamily Sericothripinae is sister to subfamily Thripinae. More molecular data from different taxonomic groups is needed to understand thrips phylogeny and evolution.

To date, the highly rearranged mitogenomes of five thrips species (Anaphothrips obscurus, Frankliniella intonsa, Frankliniella occidentalis, Scirtothrips dorsalis and Thrips imaginis) are available [17][18][19][20][21] . However, the availability of thrips mitogenomes is limited to the subfamily Thripinae. In this study, we sequenced the complete mitochondrial genome of N. samayunkur, a member of subfamily Sericothripinae using next-generation sequencing (NGS) technology and compared it to other thrips mitogenomes, analysing genome organization, codon usage patterns, tRNA secondary structure and strand asymmetry. Phylogenetic relationships were inferred by analysing the 13 PCGs from published thrips mitogenomes using maximum likelihood (ML) and Bayesian inference (BI).

Results and Discussion
Genome structure, organization and composition. The complete sequence of the mitochondrial genome of N. samayunkur (accession number MF991901) is 15,295 base pairs (bp) in length. This is longer than those of A. obscurus, F. intonsa, F. occidentalis, and S. dorsalis South Asia (SA1), but smaller than the genomes of T. imaginis and S. dorsalis East Asia (EA1) ( Table S1). As in other insect species, the mitochondrial genome of N. samayunkur included 37 genes: 13 PCGs, large and small rRNAs and 22 tRNAs but two putative CRs (Fig. 1). There are 204 bp intergenic nucleotides in total, across 23 locations, with individual spacer length of 1 to 41 bp. The longest intergenic spacer (41 bp) was located between the trnS2 and cox1 gene, with an extremely high AT content (85.37%). Three pairs of genes overlap with lengths ranging from 1 to 24 bp. Thirty genes are on the majority strand and seven on the minority (Table 1). Nucleotide composition was 77.42% A + T content and 22.58% G + C content ( Table 2); similar to other thrips mitogenomes. A + T content was highest at 80.87%, in tRNAs, followed by rRNAs (79.31%), PCGs (77.15%), and CRs (71.12%). The mitogenome showed weakly positive AT (0.04) and negative GC (−0.03) skews (Table 2).
Protein-coding genes. All 13 PCGs used ATN start codons (five with ATA, four with ATT, three with ATG and one with ATC) as is observed in most of the insect mitochondrial genomes 22,23 . The stop codon TAA was used by 10 PCGs, and TAG for atp6, while an incomplete stop codon is present in nad1 and nad2. Comparative analysis of start and stop codons among thrips showed the unique features of N. samayunkur: ATT start codon in cytb and nad6, ATC in nad2, and ATG in atp6. The complete stop codon TAA was used by atp8 in N. samayunkur, while it was terminated by an incomplete stop codon T(AA) in other thrips species (S2 Table). The detection of Figure 1. The mitochondrial genome of the marigold thrips, N. samayunkur. Direction of gene transcription is indicated by arrows in entire complete genome. PCGs are shown as purple arrows, rRNA genes as sea green arrows, tRNA genes as blue arrows and CR regions as red rectangles. The GC content is plotted using a black sliding window, as the deviation from the average GC content of the entire sequence. GC-skew is plotted using a colored sliding window (green and orchid color), as the deviation from the average GC-skew of the entire sequence. The figure was drawn using CGView online server (http://stothard.afns.ualberta.ca/cgview_server/) with default parameters. The species photograph was taken by the second author (KT) using Leica Microscope DM1000 with Leica software application suite (LAS EZ) and edited manually in Adobe Photoshop CS 8.0. an incomplete stop codon in atp8 gene may be due to misannotation, as atp8-atp6 is a conserved ancestral gene block with no tRNA between them 24 . The entire length of PCGs of N. samayunkur was 10,982 bp. Overall A + T content of 13 PCGs was 77.15% in N. samayunkur, while it ranges from 74.53% to 77.29% across thrips. Codon usage in N. samayunkur shows a significant bias towards A/T rich codons. Relative synonymous codon usage analysis of N. samayunkur revealed that the codon GCG (Alanine) was not present at all. The most frequently utilized amino acids were Lysine (K), Phenylalanine (F), Leucine (L), Isoleucine (I), Tyrosine (Y), and Serine (S) as in other insects (S3 Table).
Ribosomal and transfer RNA genes. N. samayunkur has two rRNAs as in other insects. The large ribosomal gene (16S) was 1085 bp long, and located between nad6 and trnS1; the small (12S) was 723 bp long, located between trnF and atp8 ( Table 2). A + T content of two rRNAs was 79.31%, while it ranges from 79.16% (S. dorsalis EA1) to 79.87% (F. occidentalis) observed in other thrips. N. samayunkur contained a complete set of 22 tRNAs (total length 1,401 bp) individually ranging from 58 to 72 bp in length. Collectively tRNAs have the highest A + T content 80.87% of any gene group (78.05% to 80.87% in thrips) ( Table 2). Most tRNAs have the typical cloverleaf secondary structure except trnV, trnS1, trnR, and trnT. The DHU stem and loop were absent in trnV and trnS1 while TΨ C loop absent in trnR and trnT (Fig. S1). The absence of DHU stem and loop in trnV is consistent across all thrips species sequenced to date.

Control regions.
Control regions (CRs) in mitogenomes play an important role in transcription and replication 25 . A CR was found with following conserved elements; a poly T stretch at the 5′ end, a TA(A)n-like stretch, a stem and loop structure, a TATA motif, and a GAT motif [26][27][28] . The N. samayunkur mitogenome contains two putative control regions, CR1 (628 bp) and CR2 (300 bp), located between trnV and nad6, and trnV and nad4L respectively. CR2 had 99% sequence similarity with CR1, indicating a possible duplication. Three near tandem repeats (57 bp) were identified in CR1, while one repeat sequence was present in CR2 (Fig. 2). Most thrips species have been documented to have multiple CRs except A. obscurus. Three CRs are present in F. intonsa, F. occidentalis, and S. dorsalis SA1, two in T. imaginis, S. dorsalis EA1, and one in A. obscurus [18][19][20][21] . A location of CR1 upstream of nad5 gene has been suggested to be ancestral condition of thrips 17 , however, the CR locations in N. samayunkur (subfamily Sericothripinae) differ from those of other thrips. Gene arrangement. The mitogenome gene arrangements have been characterized by following patterns, transpositions, inversions, and inverse transpositions 11,29,30 . Tandem duplication-random loss (TDRL) is the most widely accepted process to explain transpositions 11 . The gene arrangement of N. samayunkur was assessed by comparing the common intervals with the ancestral insect gene order as an outgroup 17,31 . CREx 32 analysis identified eight gene rearrangement events in N. samayunkur, including four inversions plus four TDRLs, assignable to two sets of alternative scenarios (Fig. S2). CREx detected inversions of atp8, trnF, trnC and gene block nad1-rrnS in both scenarios. N. samayunkur is a highly rearranged mitogenome with rearrangements of 11 PCGs, 22tRNAs, and two rRNAs as compared with the ancestral insect (Fig. 3). The majority of rearrangements were transpositions, while nine rearrangements (nad1, atp8, trnF, trnL1, trnQ, trnV, trnC, rrnS, and rrnL) were inverse transpositions. Further, when N. samayunkur was compared to other thrips species, the following derived gene blocks: trnG-cox3, trnN-trnE, trnY-nad1, trnF-atp6, and nad5-nad4L were conserved in all thrips species. Within the conserved gene block trnF-atp6, atp8 was subsequently inverted in N. samayunkur. The following tRNAs were inverted in thrips species as compared to the ancestral insect: trnY in S. dorsalis SA1, trnP in both S. dorsalis, trnS1 in T. imaginis, and trnF in all species except S. dorsalis SA1. The gene trnL2 was transposed in N. samayunkur away from the gene block trnL2-cox2, which is conserved in most insects including thrips (Fig. 3). The gene block trnD-cox3 is conserved in five thrips species including N. samayunkur, while trnD and trnR were translocated in T. imaginis and interrupted by the CR2 in S. dorsalis EA1. The gene block nad4L-nad5 is ancestral in insects and conserved in all thrips species. The conserved gene blocks trnY-nad1 and atp6-trnF are separated by trnM and trnA in most of the thrips species except N. samayunkur (trnM alone) and A. obscurus (trnA alone).
Strand asymmetry. AT and GC skews on the majority strand are used to measure strand asymmetry [33][34][35] .
Most insects have positive AT skew (A > T) and negative GC skew (C > G). A reversal of the strand asymmetry (T > A and G > C) has been observed in a few species, and is proposed to be caused by the inversion of the replication origin within the control region [26][27][28]35 . N. samayunkur showed weakly positive AT skew (0.04) and negative GC skew (−0.03), similar to most other thrips species (Table 2). AT skew in other thrips species, ranges from

Phylogenetic analysis. Both Maximum likelihood (ML) and Bayesian Inference (BI) phylogenetic trees
resulted in similar topologies (Fig. S3, Fig. 4). BI posterior probabilities (PP) were higher than ML bootstrap support (BS) values. It has been suggested that PP and BS values are not directly comparable and interchangeable, as PP gives higher nodal support than BS 36 . Species within the same genus, F. intonsa and F. occidentalis were grouped together and closely related to T. imaginis. The two cryptic species of S. dorsalis (EA1 and SA1) also clustered together and were closely related to the Frankliniella + Thrips clade. The four genera in the subfamily Thripinae: Anaphothrips, Frankliniella, Scirtothrips, and Thrips grouped together. N. samayunkur shows a sister relatioship to the Thripinae clade in the present phylogeny. Although, gene order is extensively rearranged among thrips mitogenomes, the branching pattern inferred by MLGO 37 is congruent with the PCGs based ML and BI phylogeny (Fig. 5). Previous studies showed a close relationship between T. imaginis with S. dorsalis 17 , however, our study found that T. imaginis was closer to Frankliniella than to Scirtothrips, congruent with morphological understanding of these taxa 38 . Previous studies described a close relationships between Scirtothrips (subfamily Thripinae) and Neohydatothrips (subfamily Sericothripinae) using both morphological 6 and molecular data 8,9 . These two genera share the following morphological characters: 6 presence of closely spaced rows of microtrichia on lateral thirds of abdominal tergites; median pair of tergal setae close together; campaniform sensilla absent on tergite IX, tergite X not split longitudinally. The suborder Terebrentia is typically classified into eight families with four subfamilies 39 ; however, an alternative view proposes 28 families based on highly conserved taxonomic characters and elevates the subfamily Sericothripinae to family rank 40 . According to Bhatti, the proposed family Sericothripidae can be separated from other families of suborder Terebrantia by the presence of sublateral callosities on antecostal line on tergites II to VII and sternites II to VI (II to VIII in male), prominent anteriorly directed sublateral apodeme on each side on female sternites VII, one pair of cervical sclerites, annular rows of microtrichia on femora and   tibiae, short and straight hind coxal apodeme, metathoracic furca not elongate and lyre-shaped, metasternellum somewhat to strongly enlarged, forming a transversely striate area on each side of the mid line, mesothorax with sternal coxal process small and inconspicuous; trochantin small and inconspicuous, absence of sclerites anterior to mesoacrotergite, one anteriorly directed apodeme on each side of sternite I. Moreover, previous studies clearly stated that the relationships between Thripinae (1779 species in 234 genera) and Sericothripinae (168 species in 3 genera) were unclear due to the absence of molecular data 8,41 . The present phylogenetic analysis contradicts a close relationship between Scirtothrips and Neohydatothrips.
The mitogenomes of two cryptic species of S. dorsalis (SA1 and EA1) vary considerably with respect to gene rearrangement and chromosome size 20 . S. dorsalis was described from "castor and chillies" at Coimbatore, India 42 . It is a polyphagous pest and a vector of tospoviruses with a global distribution. Earlier studies indicated that this species is a complex, consisting of many morphologically indistinguishable species 43,44 . Recently, nine cryptic species of S. dorsalis were delimited using multilocus molecular data 45 , however, these cryptic species has never been morphologically treated to validate and describe these species. Tagging specimens with the correct species name is a major problem as it is difficult to ascertain which of these cryptic species represent the true S. dorsalis.
To date, molecular phylogenetic studies of thrips is in its early stages due to lack of large scale data and taxonomic sampling.The generation of comprehensive molecular data on families/subfamilies is still needed.

Materials and Methods
Sample collection and DNA extraction. Adult specimens of N. samayunkur were collected from Odisha State, India. The studied species are common pests of crops, thus no prior permission was required for collection. Specimens were morphologically identified by the second author (K.T.) with available taxonomic keys 2,10 , and preserved in absolute ethyl alcohol at −30 °C in Centre for DNA Taxonomy, Molecular Systematics Division, Zoological Survey of India, Kolkata. Genomic DNA was extracted using DNeasy (QIAGEN) following the manufacturer's standard protocol. Concentration of DNA was determined using a Qubit fluorometer with a dsDNA high-sensitivity kit (Invitrogen), and by agarose gel (0.8%) electrophoresis.
Mitogenome sequencing and assembly. The whole genome library of genomic DNA was sequenced using the Illumina Hiseq2500 (2 × 150 base paired-end reads) (Illumina, USA) platform which yielded ~23 million reads. The paired-end library was constructed according to standard protocols for the TruSeq DNA Library Preparation kit (https://support.illumina.com/downloads/truseq). Raw sequencing reads were trimmed and quality filtered using the NGS-Toolkit 46 to removing adapter contamination and low-quality reads (N's or more than 70% of bases with a quality score < 20). High quality reads were filtered by using the Burrows-Wheeler Alignment (BWA) tool 47 and assembled with SPAdes 3.9.0 48 , using default parameters, and the S. dorsalis mitochondrial genome (NC_025241.1) as a reference. Aligned reads were used for de novo mitochondrial genome assembly.
Genome annotation, visualization, and comparative analysis. The assembled mitogenome was annotated using the MITOS web-server (http://mitos.bioinf.uni-leipzig.de/index.py) 49 . PCGs and rRNAs were confirmed manually by BLASTn, BLASTp and ORF Finder in NCBI 50,51 (https://www.ncbi.nlm.nih.gov/ orffinder/). Nucleotide sequences from protein coding genes (PCGs) were translated into putative proteins on the basis of the invertebrate mitochondrial genetic code. Initiation and termination codons were identified in ClustalX 52 using other thrips reference mitogenome sequences. MEGA6 53 was used for the alignment of homologous sequences across thrips species. The complete annotated mitogenome was submitted to NCBI GenBank using Sequin tool (http://www.ncbi.nlm.nih.gov/Sequin/). The circular map of the N. samayunkur mitogenome was illustrated by the CGView online server (http://stothard.afns.ualberta.ca/cgview_server/) with default parameters 54 . MEGA6 was used for estimation of nucleotide composition, codon usages, relative synonymous codon usage (RSCU) and composition of skewness with the following formula: AT skew = (A − T)/ (A + T) and GC skew = (G − C)/(G + C) 55 . Secondary structures of transfer RNA (tRNA) genes were predicted by MITOS and further confirmed using tRNAscan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/) 56 and ARWEN 1.2 57 . RNAstructure version 6.0.1 was used to predict possible secondary structure within CRs 58 . Homology between CR1 and CR2 in N. samayunkur was determined through the ClustalW sequence alignment tool implemented in MEGA6. Gene arrangements pathways in N. samayunkur were evaluated by CREx (Common Interval Rearrangement Explorer) 32 .
Phylogenetic analysis. Six complete mitogenomes of five thrips species were retrieved from GenBank on 1 st November 2017 for phylogenetic inference (S1 Table). The A. bakeri mitogenome was used as an out group 31 . Each PCG was aligned individually using the MAFFT algorithm in the TranslatorX 59 online platform under the L-INS-i strategy based on codon-based multiple alignment. Poorly aligned nucleotides (1652 bp) were removed from the protein alignment using GBlocks (within TranslatorX) with default settings. The resulting alignments were concatenated by using Sequence Matrix1.7.8 60 . Concatenated dataset (9330 bp) was used for Bayesian inference (BI) and maximum likelihood (ML) analysis. PartitionFinder version 2.1.1 61 , with the greedy algorithm was used to find the best substitution models and partition schemes. Partitions were predefined for the codon positions for each PCGs (13 genes X 3 codons = 39 partitions). The BI analysis was performed using Mr. Bayes 3.2 62 with HKY + I + G, TVM + G, TRN + G, GTR + I + G, HKY + G, GTR + I, TVM + I + G model estimated by PartitionFinder (S4 Table). Two runs each with four chains (three heated and one cold) for 500,000 generations, and trees were sampled every 100 generations. A consensus tree was acquired and visualized after excluding the first 25% trees as burn-in. The ML analysis was performed using the IQ-TREE 63 Web Server in W-IQ-TREE 64 (http://iqtree.cibiv.univie.ac.at/) with 1,000 replicates of ultrafast likelihood bootstrap 65 66 . Phylogenetic relationships of studied taxa were also estimated based on gene arrangement patterns in the MLGO web server 37 (S5 Table).