Complete Mitochondrial Genome of Helicoverpa zea (Lepidoptera: Noctuidae) and Expression Profiles of Mitochondrial-Encoded Genes in Early and Late Embryos.

The mitochondrial genome (mitogenome) of the bollworm, Helicoverpa zea (Boddie), was assembled using paired-end nucleotide sequence reads generated with a next-generation sequencing platform. Assembly resulted in a mitogenome of 15,348 bp with greater than 17,000-fold average coverage. Organization of the H. zea mitogenome (gene order and orientation) was identical to other known lepidopteran mitogenome sequences. Compared with Helicoverpa armigera (Hübner) mitogenome, there were a few differences in the lengths of gaps between genes, but the lengths of nucleotide overlaps were essentially conserved between the two species. Nucleotide composition of the H. zea mitochondrial genome was very similar to those of the related species H. armigera and Helicoverpa punctigera Wallengren. Mapping of RNA-Seq reads obtained from 2-h eggs and 48-h embryos to protein coding genes (PCG) revealed that all H. zea PCGs were processed as single mature gene transcripts except for the bicistronic atp8 + atp6 transcript. A tRNA-like sequence predicted to form a hammer-head-like secondary structure that may play a role in transcription start and mitogenome replication was identified within the control region of the H. zea mitogenome. Similar structures were also found within the control regions of several other lepidopteran species. Expression analysis revealed significant differences in levels of expression of PCGs within each developmental stage, but the pattern of variation was similar in both developmental stages analyzed in this study. Mapping of RNA-Seq reads to PCG transcripts also identified transcription termination and polyadenylation sites that differed from the sites described in other lepidopteran species. Published by Oxford University Press on behalf of the Entomological Society of America 2016. This work is written by US Government employees and is in the public domain in the United States.

Helicoverpa zea (Boddie), commonly known as the bollworm or corn earworm, is a polyphagous insect in the family Noctuidae and the only species of Helicoverpa in the subfamily Heliothentinae (Hardwick 1996) in North America. Its geographic range includes most temperate and tropical regions of the Americas with the exception of Alaska and northern Canada (Metcalf and Flint 1962, Hardwick 1965) and overlaps with the range of Helicoverpa gelotopoeon Dyar in South America. Both have very similar host ranges, and morphology of the larval stages is similar. Corn earworm damages a large number of agricultural crops and estimated annual losses in the United States alone exceed one billion dollars (Light et al. 1993). Control costs for this pest are a significant component of crop production costs in the United States and insecticides are routinely used for suppression of H. zea in cotton (Williams 2014) and soybean (Musser et al. 2013). Helicoverpa zea is one of the primary targets of transgenic corn hybrids and cotton varieties that express insecticidal proteins derived from Bacillus thuringiensis Berliner (Luttrell and Jackson 2012). Of the 10 other species of Heliothentinae (Hardwick 1996) in North America, Heliothis virescens (Fabricius) is the only other significant pest of agricultural crops.

Both H. zea and Helio. virescens are polyphagous and often attack the same crops at the same time. Eggs and larvae are almost indistinguishable in the field. Morphological traits can be used to identify both species, but the process is labor intensive and requires inspection under magnification. Diagnostic tools for larvae based on monoclonal antibodies and genetic differences are available and have been used on a limited scale to characterize mixed populations in field situations (Zeng et al. 1999, 1998). An expanding base of genetic knowledge and genetic markers (Grasela and McIntosh 2005, Behere et al. 2007, Perera et al. 2007) has provided additional diagnostic tools for species identification and population genetics.

The recent introduction of Helicoverpa armigera (Hübner) into Brazil (Czepak et al. 2013, Tay et al. 2013), further spread into Argentina, Paraguay, and Uruguay (Murúa et al. 2014), and potential future movement into North America further complicates identification of the Heliothentinae pest species. The establishment of H. armigera in the United States represents a significant risk of increased economic costs to agriculture, both in control costs and in direct crop losses. Inter-specific hybridizations between heliothine species that yield fertile progeny have been documented. Partial sterility was observed between Helio. virescens and Heliothis subflexa Guenée (Laster et al. 1988), H. armigera and Helicoverpa assulta (Guenée) (Wang 2007), and H. zea and H. assulta (Guenée) (Degrugillier and Newman 1993). However, no evidence of sterility or incompatibility was identified between H. armigera and H. zea (Laster and Hardee 1995, Laster and Sheng 1995). Any hybridization between invasive H. armigera and the two Helicoverpa species native to the Americas resulting in nuclear gene introgression would further complicate taxonomic identification based on morphology as well as nuclear genes. Mitogenome analysis facilitates the positive identification of the maternal lineage of any suspect hybrid and can be compared with nuclear markers.

Insect mitochondrial (mt) genomes tend to range from 14 to 18 kb in size and usually possess 37 genes; 13 protein coding genes (PCGs), 22 tRNAs, 2 rRNA genes, and several noncoding regions (Boore 1999). Mitochondrial DNA (mtDNA) sequences have been the tool of choice in a large number of studies on population genetics and phylogenetics of insects (Castro et al. 1998, Caterino et al. 2000). This is largely due to the highly conserved nature of certain regions that have allowed the development of universal primers. However, these highly conserved regions of the mitochondrial genomes are not ideal for intraspecific genetic studies as they may underestimate actual genetic diversity. The use of entire mitogenome sequences is increasing in intraspecies genetic studies as they can provide greater resolution (Yu et al. 2007, Teacher et al. 2012, Pabijan et al. 2013).

A large number of lepidopteran mitogenomes have been sequenced to date (>100 NCBI) from diverse taxonomic groups. Missing from many of the sequenced mitogenomes is information on expression of the genes and confirmation of predicted annotations using expressed gene transcripts. In most vertebrate and invertebrate mitochondrial genomes there is little or no intergenic space and PCGs are interspersed by tRNA genes, the usual transcript cleavage sites (Ojala et al. 1980). However, the mitochondrial genomes can also produce polycistronic protein coding RNA transcripts where there is overlap within the genome sequence. This has been identified in a number of insect species including Maruca vitrata (Fabricius) and Drosophila melanogaster Meigen but with some differences in the genes identified as polycistronic (Stewart and Beckenbach 2009, Margam et al. 2011). Analysis of transcriptome data can help accurately annotate these regions and identify genes that are polycistronic and those that are not. Mapping of mitochondrial transcripts, expression mechanisms, and gene expression profiles have been studied in various animal species (Costanzo and Fox 1990, Stewart and Beckenbach 2009, Torres et al. 2009, Wang et al. 2013, Marková et al. 2015), but so far only one study has examined the transcript processing in a moth species (Margam et al. 2011). A reliable means of identifying samples to species and performing intragenetic studies would be useful to properly understand the distribution, ecology, agricultural impacts, and the insecticide resistance status in these heliothine pest species. The full-length mitogenome of H. zea is a useful tool in developing markers and tests suitable for species- specific and intraspecific studies. In this work, we determined the nucleotide sequences of the complete mitogenome of the moth, H. zea and used RNA-Seq data to study processing of polycistronic mitochondrial transcripts.

Materials and Methods

DNA Extraction and Nucleotide Sequencing

A thorax of a single adult female H. zea from a laboratory colony maintained at the USDA-ARS Southern Insect Management Research Unit, Stoneville, MS was used in mtDNA extraction. Tissue was homogenized in a buffer containing 90 mM KCl, 55 mM CaCl2, 30 mM NaCl, 15 mM MgSO4, and 250 mM sucrose and filtered using a nylon mesh cloth. The filtrate was spun at 200 × g for 25 min. to pellet nuclei and large tissue debris. Supernatant enriched with mitochondria was carefully aspirated and transferred to a new 1.5-ml microcentrifuge tube. The supernatant was centrifuged at 5,000 × g to pellet mitochondria. DNA was extracted from mitochondria using MasterPure DNA extraction kit (Epicentre, Charlotte, NC) following manufacturer’s instructions. Nextera XT (Illumina, San Diego, CA) was used to prepare a small insert sequencing library from the mtDNA enriched nucleic acid and 100 bp paired end reads were generated on a MiSeq (Illumina) at the USDA-ARS Genomics and Bioinformatics Research Unit.

A de novo assembly of the paired-end miSeq reads was performed using SeqMan NGen 12.0 (DNAStar, Madison, WI). Assembly parameters minimum match percentage, match size, match spacing, mismatch penalty, gap penalty, and maximum gap length were set to 93, 50, 10, 20, 30, and 6%, respectively. Consensus sequence was exported and ends were manually edited to remove duplicated nucleotides.

Annotation of the Mitogenome

tRNA genes were identified using the tRNAscan-SE Search Server (Lowe and Eddy 1997, Schattner et al. 2005). PGCs, rRNA genes, and the tRNAs that were not identified by the tRNAscan-SE program were annotated by comparing with published mt genomes of H. armigera and Helicoverpa punctigera Wallengren (Yin et al. 2010, Walsh 2014). Nucleotide compositional skew for the whole genome was calculated using the formula GC-skew = (G − C)/(G + C) and AT-skew = (A − T)/(A + T), where A, C, G, and T are the total number of nucleotides represented by the respective IUPAC letter code in the nucleotide sequence (Perna and Kocher 1995).

RNA-Seq Library Preparation

Freshly emerged females and males (n = 50) from a laboratory colony were placed in a 4-liter plastic bucket and allowed to mate for 3 d. Eggs were collected on the fourth day by placing a sheet of wax paper over the open end of the container for 1 h. Twenty-five eggs were dislodged from the wax paper using a small brush at the end of the second h (2-h eggs) and total RNA was extracted using TriZol reagent (Invitrogen, Carlsbad, CA) following manufacturer’s instructions. The remainder of the eggs was incubated at 27°C for a total of 48 h. At the end of incubation period, 25 eggs containing developed embryos with dark head capsules (48-h embryos) were selected using a microscope and dislodged from the wax paper. TriZol reagent was used to extract total RNA from these eggs. Freshly extracted total RNA was treated with 1 unit of RNase-free DNase per 2 µg of total RNA. Double-stranded cDNA was synthesized from 2 µg of total RNA using SuperScript III cDNA synthesis reagents (Invitrogen) following manufacturer’s instructions. First strand cDNA synthesis was primed using a custom-designed oligo d(T) primer (5′-GGT AA T ACGACTCACTATAGGGAGAAGAGGCGAGC ACAGAATTA AT ACGACTTTTTTTTTTTTTTTTTTTTV-3′). After the second strand synthesis, cDNA was purified and concentrated and Illumina sequencing libraries were prepared with 50 ng of double-stranded cDNA using Nextera library construction reagents. Two replicate libraries each were prepared for 2-h eggs and 48 embryos. Library quality was evaluated by running 1 µl of each library on an Agilent 2200 Tape Station (Agilent Technologies, Santa Clara, CA). The libraries were size fractionated to select 400–500 bp fragments using the Pipin Prep instrument (Sage Science, Inc., Beverly, MA). One-hundred base pair single-end reads were generated using the Illumina HiSeq 2500 instrument at the USDA-ARS Genomics and Bioinformatics Research Unit.

Transcript and Expression Analysis

Complete nucleotide sequences of 13 PCGs were extracted from H. zea mitogenome and used in calculating synonymous and nonsynonymous substitution rates compared with orthologous mitochondrial PCGs of H. armigera. The number of nonsynonymous substitutions per nonsynonymous site (KA) and the number of synonymous substitutions per synonymous site (KS) between the two Helicoverpa species were estimated using the software KaKs Calculator (Zhang et al. 2006) with the maximum-likelihood model averaging method (Posada 2003).

Mitochondrial PCG sequences were also used for mapping RNA-Seq reads generated from 2-h eggs and 48-h embryos. SeqMan NGen 4.0 assembly parameters used in this reference assembly were identical to those used in the reference assembly of the mitogenome. Reads mapping to each gene were manually examined to verify ends of the genes and the presence or absence of polyadenylation. Coding sequences of the 13 PCG were also used in expression profiling using ArrayStar 5.0 software (DNAStar). Reads per kilobase per million was used to normalize the expression levels of the PCG. Normalized expression values and the standard deviations of the PCG were exported from ArrayStar 5.0 software and each expression value was divided by the value of the gene with the lowest expression level (nad4L). Covariance (CV) of the quotient was calculated by the formula , where SD1 and SD2 are the standard deviations of the two expression values used to calculate the relative expression level. Standard deviations of the relative expression values were calculated by multiplying the relative expression value by the CV of each relative expression value.

In addition, all Illumina sequence reads from 2-h eggs and 48-h embryos were mapped separately to full-length sequence using CLC Genome Workbench 8.5 (Qiagen, Redwood City, CA) to obtain sequence read density profiles across the complete mitogenome. Mapped reads were analyzed for sequence polymorphisms using a basic variant detection tool implemented in CLC Genome Workbench.

Results

Nucleotide Composition and Genome Organization

The mtDNA enriched H. zea sequencing library yielded 2,226,753 sequence read pairs with Phred quality values ≥35. Nucleotide sequence reads were deposited in the NCBI SRA database (Submission no. SRX690295). De novo assembly of paired-end miSeq reads yielded a contig of 15,984 bp with 17,173-fold average coverage from a total of 1,447,494 read pairs consistent within the contig. Manual editing of duplicated sequences at the ends resulted in a final H. zea mitogenome of 15,348 bp (Accession no. KJ930516). The mitogenome of H. zea is 81.00% A + T, similar to 80.97% in H. armigera and 81.35% in H. punctigera (Yin et al. 2010, Walsh 2014). The GC-skew and the AT-skew in the complete mitogenome of H. zea were −0.2032 and 0.0024, respectively, indicating that this mitogenome contains more C than G and more A than T. This nucleotide skew pattern is similar to that of closely related species H. armigera and H. punctigera. Overall nucleotide identity between H. zea mitogenome and those from H. armigera and H. punctigera were 97.4 and 94.6%, respectively.

Annotation of the mt genome of H. zea identified 13 PCG, 2 rRNA genes, 22 tRNA genes, and a control region with high A + T content and repetitive sequences. The gene order and the orientation of the genes in H. zea mitogenome were identical to those of H. armigera and other lepidopteran species. Organization of the genes in H. zea mitogenome is shown in Fig. 1 and the coordinates of the genes and control region are given in Table 1. Although there are a few differences in the lengths of space (i.e., number of nucleotides) between some of the mitochondrial genes, the lengths of overlaps between trnI and trnQ, trnW and trnC, atp8 and atp6, atp6 and cox3, trnE and trnF are conserved between H. armigera and H. zea (Table 1). One exception noted in this study is the processing of nad3 and trnA genes. Although there is a native in-frame TAG stop codon (nucleotides 5923–5925) in the nad3 gene, the stop codon in all processed transcripts of (i.e., all polyadenylated transcripts) found in this study contained a TAA stop codon created by polyadenylation. The trnA sequence starts at the nucleotide position 5924 and the excision of the trnA at this nucleotide would essentially remove the TAG stop codon, and generation of the TAA stop codon when the polyA+ tail was added. There were one to two nucleotide differences in length between trnY, trnN, trnE, trnT, trnS2 genes of H. zea and H. armigera and the rrnL gene in H. zea is 5 bp shorter than that of H. armigera due to one to two nucleotide indels distributed across the rrnL gene.

Fig. 1.

A graphical illustration of gene organization in the mitogenome of Helicoverpa zea. The outer circles represent the forward or heavy strand and the inner circle represents the reverse or light strand of the mitogenome. The protein coding sequences and rRNA are marked in between the circles and the genes coding for tRNA are marked outside the circles. Underlined text indicates gene coded by the light strand of DNA. Putative palindromic repeat containing trnI-like sequence predicted inside the control region is denoted by I*.

Table 1.

Organization of Helicoverpa zea mitogenome

Gene name	Start	Stop	Length	Start	Stop	tRNA anticodon	Overlap/gap (bp)
trnM	1	69	69			CAT	0
trnI	70	133	64			GAT	−3
trnQ	199	131	69			TTG	45
nad2	245	1,255	1,011	ATT	TAA		2
trnW	1,258	1,325	68			TCA	−8
trnC	1,386	1,318	69			GCA	10
trnY	1,462	1,397	66			GTA	9
cox1	1,465	2,995	1,531	CGA	T+polyA+		0
trnL2	2,996	3,062	67			TAA	0
cox2	3,063	3,744	682	ATG	T+PolyA+		0
trnK	3,745	3,815	71			CTT	4
trnD	3,820	3,885	66			GTC	0
atp8	3,886	4,047	162	ATC	TAA		−7
atp6	4,041	4,718	678	ATG	TAA		−1
cox3	4,718	5,503	786	ATG	TAA		2
trnG	5,506	5,571	66			TCC	0
nad3	5,572	5,925	354	ATT	T+polyA+		−2
trnA	5,924	5,992	69			TGC	0
trnR	5,993	6,059	67			TCG	27
trnN	6,087	6,152	66			GTT	1
trnS1	6,154	6,221	68			GCT	0
trnE	6,222	6,288	67			TTC	−2
trnF	6,353	6,287	67			GAA	0
nad5	8,099	6,354	1,746	ATA	TAA		0
trnH	8,166	8,100	67			GTG	0
nad4	9,505	8,167	1,339	ATG	T+PolyA+		0–2
nad4l	9,840	9,550	291	ATG	TAG		3
trnT	9,844	9,909	66			TGT	0
trnP	9,974	9,910	65			TGG	7
nad6	9,982	10,515	534	ATT	TAA		0
cob	10,516	11,667	1,152	ATG	TAA		6
trnS2	11,673	11,738	66			TGA	23
nad1	12,700	11,762	939	ATG	TAA		1
trnL1	12,769	12,702	68			TAG	0
rrnL	14,159	12,770	1,390				0
trnV	14,225	14,160	66			TAC	0
rrnS	15,019	14,226	794				0
Control	15,020	15,348	329
trnI-likea	15,173	15,268				AAT

Start and stop nucleotide positions are given for all genes. Translation start and stop codons for open reading frames of PCGs and the anti-codons for tRNA genes are given in respective columns. The number of nucleotide in a gap or an overlap between two successive genes is indicated by a positive or negative number, respectively.

Putative tRNA-like sequence with palindromic repeats.

Results of the comparative analysis H. zea and H. armigera of PCGs are given in Table 2. Alignments of the transcripts and translated peptide sequences of the PCGs of two species are given in Fig. S1. The ML models tested whether the KA/KS ratio for a given sequence pair significantly deviate from 1 (i.e., neutral selection). Significantly low KA/KS ratios (P < 0.00001) were observed for all PCG comparisons except for atp8 and nad4l, indicating that stabilizing or purifying selection pressure on most mitochondrial PCGs. The highest number of nonsynonymous substitutions was found in nad5 (6), nad6 (6), nad1 (4), nad4l (4), and cob (3). Nonsignificant test results indicated that the KA/KS ratios of atp8 and nad4l genes did not deviate from KA/KS ratio of 1 and may be selection neutral. The power of KA/KS ratio test was dependent on the sequence length and the degree of divergence between the sequences being compared (Nekrutenko et al. 2002). Therefore, short sequence length of atp8 and nad4l genes may have contributed to the nonsignificant test results.

Table 2.

Total nucleotide substitutions, maximum-likelihood (ML) estimates of synonymous (S) and nonsynonymous (N) substitutions, proportion of synonymous substitutions per synonymous site (KS), proportion of nonsynonymous substitutions per nonsynonymous site (KA), the KA/KS ratio, and the P-value of the Fisher exact test for the 13 mitochondrial PCGs of Helicoverpa zea compared with the orthologous genes of H. armigera

PCG	Actual substitutions	ML estimated S-substitutions	ML estimated N-substitutions	KA	KS	KA/KS	P-value (Fisher)
nad2	24	21.124	2.876	0.00536	0.64407	0.00833	1.8014E−27
cox1	39	38.529	0.472	0.00064	0.47556	0.00134	1.6880E−42
cox2	21	20.728	0.272	0.00069	0.69224	0.00100	2.0111E−27
atp8	3	0.659	2.341	0.02107	0.04452	0.47329	5.6586E−01
atp6	16	15.908	0.092	0.00022	0.21872	0.00100	5.6774E−15
cox3	21	20.861	0.139	0.00031	0.31421	0.00100	1.3740E−20
nad3	13	12.928	0.072	0.00035	0.34691	0.00100	0.0000E+00
nad5	35	31.137	3.863	0.00412	0.39477	0.01044	4.0051E−34
nad4	27	26.787	0.213	0.00197	1.77085	0.00111	1.3641E−26
nad4l	5	1.554	3.446	0.01416	0.08010	0.17675	4.5459E−02
nad6	18	11.706	6.294	0.01829	0.31909	0.05732	2.1480E−09
cob	39	36.792	2.208	0.00385	0.46802	0.00823	5.7313E−34
nad1	31	29.368	1.633	0.00490	1.05782	0.00464	1.5389E−35

Mitochondrial Transcripts

In the present study, a cDNA library produced by reverse transcription using an oligo d(T) primer was used to generate sequence reads from the polyadenylated 3′-ends, facilitating identification of polyadenylation sites of the PCG transcripts. Illumina sequence reads obtained from both developmental stages were submitted to NCBI SRA database (Project no. PRJNA259919). Reference assembly of RNA-Seq reads from 2-h eggs and 48-h embryos yielded 411,145 and 168,606 reads mapped to the reference mitochondrial gene set, respectively. Examination of the RNA-Seq reads mapped to the PCG revealed that in addition to polycistronic primary transcripts, both 2-h eggs and 48-h embryos produced mature polyadenylated transcripts for all genes except for atp8. All the reads mapped to atp8 extended beyond the putative translation end site into the atp6 gene, indicating that polypeptides of the gene pair atp8 and atp6 are translated from bicistronic mature transcripts. In H. zea, we did not find any evidence for tricistronic atp8-atp6-cox3 transcripts reported in M. vitrata (Margam et al. 2011).

Bicistronic nad4L-nad4 transcripts may be an evolutionarily conserved feature of mitogenome transcription in other species studied so far (Berthier et al. 1986, Scheibye-Alsing et al. 2007, Stewart and Beckenbach 2009, Margam et al. 2011, Wang et al. 2013). However, we found polyadenylated transcripts of nad4L and nad4 in both 2-h eggs and 48-h embryos, indicating that the transcripts are processed into single transcripts in H. zea. Moreover, the polyadenylation most likely started 42–44 bp downstream of the nad4L translation stop codon (at the nucleotide position 9,506–9,508), but the precise location could not be predicted because we could not determine whether the two adenosines present at the nucleotide positions 9,506–9,507 were part of the polyA+ tail or not. Unlike overlapping nad4L-nad4 genes in M. vitrata and D. melanogaster, these two genes in H. zea were separated by a 44-bp gap and processing them into individual mature transcripts would not affect the translation into polypeptides.

The mitochondrial genes nad6 and cob were considered as a bicistronic transcript in other insect species. We detected polyadenylated nad6 transcripts in both 2-h eggs and 48-h embryos, clearly indicating that nad6 and cob are processed as independent mature transcripts in all developmental stages of H. zea. Differences in processing of the cob transcript were also observed between 2-h eggs and 48-h embryos. In 48-h embryos, 40.5% of 121 reads mapped to the 3′-end of the cob transcript were polyadenylated. However, we did not detect any evidence for polyadenylation in the 396 reads from 2-h eggs mapped to the 3′-end of the cob transcript.

Although the transcription termination and polyadenylation sites of mitochondrial PCG are commonly inferred by comparison with a published data from a closely related species, our study found that some of the H. zea PCG have entirely different or multiple alternative polyadenylation sites compared with those established for H. armigera and other related species. In this study, we found polyadenylation of cob transcripts started at a site different from predicted by comparison with other closely related species. Transcripts of cox3, nad1, nad2, and nad5 had alternative polyadenylation sites in addition to the site predicted by comparison with other mitogenomes (Table 3). These changes in polyadenylation sites also affect the start sites of the gene immediately downstream of the affected gene.

Table 3.

Protein coding sequences in Helicoverpa zea mitogenome with different polyadenylation start sites

Gene	Positions of alternative Ply A+ site(s)
nad2	1255; 1257
cox3	5005; 5007; 5011
nad5	6351; 6354
cob	11677
nad1	11757; 11759; 11760

Other Features

Analysis of the control region of H. zea mitogenome with tRNAscan-SE Search program indicated the presence of a putative trnI-like sequence with a 23-bp intron. Control regions of H. punctigera and Bombyx mandarina mitogenomes also contained putative trnN-like structures with possible intron sequences (Table 4). Mitogenome control regions of Spodoptera exigua and Spodoptera litura contained trnI-like (AAU) and an Ochre suppressor tRNA-like (UAA) sequences with possible intron sequences, respectively. A putative trnY-like structure without an intron was predicted in the Bombyx mori control region. Control regions of H. armigera and a few other insects did not have tRNA-like structures detectable by the tRNAscan-SE program (Table 4). Vertebrate mitogenome control regions contain nucleotide sequences forming tRNA-like secondary structures that may play a role in DNA replication (Shadel and Clayton 1997). Evaluation of the 96 bp putative H. zea trnL-like sequence with mFold (http://mfold.rna.albany.edu/?q=DINAMelt/Quickfold) indicated the presence of a hammerhead-like structure (Fig. 2). Examination of the putative tRNA-like sequences from H. punctigera, B. mandarina, B. mori, S. exigua, and S. litura with mFold revealed secondary structures similar to that of predicted trnI-like sequence in H. zea. Mapping of RNA-Seq reads from 2-h eggs and 48-h embryos to the 329 bp control region resulted in a total of 269 mapped reads yielding an average coverage of 47.6. However, no transcripts were mapped to the nucleotides between 15,186 and 15,257 of the H. zea mitogenome. This stretch of control region overlaps with the predicted tRNA-like sequence. The presence of transcripts flanking this region may indicate that the polycistronic transcripts of mitochondrial genes start in this region and the secondary structure possibly plays a role in transcription initiation as well as replication initiation.

Table 4.

Analysis of control regions of selected lepidopteran mitogenomes with tRNAScan-SE for putative tRNA-like sequences

Species	Control region length	tRNA	Intron	Accession number
Agrotis ipsilon	332	–	−	KF163965.1
Athyma sulpitia	349	–	−	JQ347260
H. punctigera	327	trnN	+	KF977797.1
Bombyx mandarina	747	trnN	+	AB070263.1
Bombyx mori	499	trnY	−	AF149768.1
Glyphodes quadrimaculalis	327	–	−	KF234079
Helicoverpa armigera	328	–	−	GU188273.1
Helicoverpa zea	329	trnI	+	KJ930516
Manduca sexta	324	–	−	NC010266.1
Mythimna separata	374	−	–	KF730242
Paracymoriza distinctali	352	–	−	KF859965.1
Spodoptera exigua	331	trnI	+	JX316220.1
Spodoptera litura	326	sup(UUA)	+	KF701043.1

A “+”or “−”in the intron column indicates presence or absence of a possible intron within the predicted tRNA-like structure.

Fig. 2.

Secondary structure predicted for putative trnI-like region identified in the control region of the Helicoverpa zea mitogenome. mFold web server was used to predict the secondary structure.

Expression Profiles in Early Eggs and Late Embryos

In total, 54,570,728 and 39,320,735 reads from 2-h eggs and 48-h embryos, respectively, were used for expression profiling and mapping reads to the complete mitogenome. The number of reads mapped to the mitochondrial PCGs and the two rRNA genes from 2-h eggs and 48-h embryos were 6,361,316 (11.6%) and 1,506,558 (3.8%), respectively. However, relative expression levels of the transcripts, except cox1, did not deviate significantly between two developmental time points. Expression profiles of PCGs indicated only 1.46- to 2.41-fold differences in expression levels of two transcripts (nad2 and cox1) between 2-h eggs and 48-h embryos (Fig. 3), but these differences did not represent physiologically significant differences in mitochondrial respiratory chain complexes between two developmental stages. Expression levels of different PCG transcripts varied broadly, but in a similar pattern (Fig. 3), in both 2-h eggs and 48-h embryos. Normalized expression levels of cox complex genes varied from 9.49 ± 1.16 for cox2 in 48-h embryos to 90.58 ± 1.28 for cox1 in 2-h eggs. High level of cox gene transcripts in both developmental stages may represent physiologically significant gene expression patterns in developing embryos or an artifact resulting from mapping short nucleotide reads to cox PCGs, especially if nucleotide regions conserved between nuclear cox genes (including pseudo cox genes expressed in the nucleus) were mapped to the reference transcripts. Further studies are needed to fully understand whether these differences represent actual gene expression patterns or experimental artifacts. Sequence read coverage of the large subunit of rRNA (rrnL) was 33- and 65-fold higher in 2-h eggs and 48-h embryos, respectively, compared with the small subunit of rRNA (rrnS). Although the oligo d(T) primer used for reverse transcription was expected to convert mRNA with a poly A-tail, rRNA contamination is common in even in cDNA synthesized with poly(A)+ selected mRNA (Sooknanan et al. 2010). Therefore, unequal representation of the two rRNA subunit transcripts was most likely due to uneven conversion of contaminating rRNA transcripts to cDNA.

Fig. 3.

Expression levels of the protein coding sequences of the mitogenome of Helicoverpa zea in 2-h-old eggs and 48-h-old embryos, relative to the expression level of nad4L. Expression levels of the PCGs were normalized to that of nad4L which had the lowest expression level.

Mapping the nucleotide sequence reads from two developmental times to the complete mitogenome also resulted in a coverage maps similar to those observed in expression analysis (Fig. 4). Peaks in the coverage maps represent high read coverage, usually in the middle of the transcript, and valleys represent the ends of transcripts with low coverage. This pattern was observed for all individual PCGs except for atp8 and atp6 transcripts where a read mapping pattern consistent with a single bicistronic atp8-atp6 transcript was observed, indicating that atp8-atp6 is the only mature multicistronic transcript present in mitochondrial of H. zea. Evaluation of the mapped reads for nucleotide sequence variants did not identify any polymorphisms in mitochondrial transcripts other than random variations due to sequencing errors. This result is not surprising because the eggs used for generating RNA-Seq libraries were obtained from a highly inbred laboratory colony.

Fig. 4.

Coverage of RNA-Seq reads from 2-h eggs (orange) and 48-h embryos (blue) mapped to complete mitogenome of Helicoverpa zea. All PCGs and rRNA subunit genes are shown below the coverage map. Peaks represent high read coverage, usually in the middle of the transcript, and valleys represent the ends of transcripts with low coverage.

Discussion

In this study, we determined the complete nucleotide sequence of the mitogenome of H. zea, identified PCG transcripts that were processed differently compared with other insect species, and evaluated expression levels of PCGs in early eggs and late embryos. Conserved lepidopteran gene organization was preserved in H. zea mitogenome with some species specific characteristics such as differences in gaps or overlaps between mitochondrial genes.

Genetic studies using mitogenomes sequences are abundant in literature (reviewed in: Moritz 1994, Cameron 2014). Segments of mitochondrial genes such as cox1 and cox2 have been used in evolutionary and population genetic studies as well as for identification of Helicoverpa species (Behere et al. 2007, 2008; Tay et al. 2013; Leite et al. 2014; Mastrangelo et al. 2014). Although these highly conserved genes provide sufficient level of divergence between species, they may not provide sufficient resolution for intraspecific genetic studies and multilocus or complete mitogenome-based analyses are more suitable for intraspecific genetic studies (Anstead et al. 2002, Zhang et al. 2012, Cordon-Obras et al. 2014, De Jager et al. 2014, Morales et al. 2015). Alignment of H. zea and H. armigera mitogenomes indicated 327 nucleotide substitutions between two species (Supplementary data S1 and S2) of which 292 (89.3%) were in the 13 PCGs (Table 2). There were 47 alignment gaps resulting from one to four nucleotide differences (indels) in nonprotein coding regions. Coding regions of cox1 and cox2, two genes commonly used in genetic studies, had 39 and 21, respectively, nucleotide substitutions between two Helicoverpa species. A recent study of H. armigera and H. zea in Brazil using cox1 gene reported low nucleotide diversity within both species (Leite et al. 2014). Multilocus analysis that included less conserved genes such as nad2, nad6, and cob would have been more suitable for studying these insects due to inherent low genetic diversity (Mallet et al. 1993, Perera and Blanco 2011).

Mitogenomes are tRNA punctuated and polycistronic mitochondrial primary gene transcripts cleaved off at the 5′- and 3′-end of each tRNA gene by mitochondrial RNase P and RNase Z, respectively, to produce mature transcripts (Rackham and Filipovska 2013). It has been suggested that when PCG are not punctuated by tRNA genes, tRNA-like secondary structures at the boundaries of the PCG are recognized by the processing machinery (Ojala et al. 1981, Montoya et al. 1983, Fernandez-Silva et al. 2003, Sanchez et al. 2011, Rackham and Filipovska 2013). Insect mitogenomes have seven genes that are not punctuated by tRNA genes. These seven genes are arranged in three contiguous clusters containing atp8-atp6-cox3, nad4L-nad4, and cob-nad6. Mitogenome transcript analyses in other insects have detected mature bicistronic transcripts containing atp8-atp6, nad4L-nad4, and cob-nad6 (Berthier et al. 1986, Falkenberg et al. 2007, Stewart and Beckenbach 2009) and a tricistronic contig containing atp8-atp6-cox3 assembled using Roche-454 sequence reads in the legume pod borer, M. vitrata (Margam et al. 2011). The only gene identified in reference assemblies without mapped PolyA+ containing reads in H. zea was atp8, which could form the bicistronic transcript for atp8-atp6 genes. We could exclude mature tricistronic atp8-atp6-cox3 transcripts in H. zea because 3′-end of the atp6 gene was polyadenylated. We manually examined raw reads mapped to reference mitochondrial transcripts to identify polyadenylated reads rather than using contigs assembled from raw reads. At any given time point, active mitochondria contain polycistronic primary transcripts as well as processed mature transcripts. Models proposed for mitochondrial transcript processing in D. melanogaster indicate immediate polyadenylation of the 3′-end of primary transcripts followed by processing of individual genes starting at the 3′-end of the primary transcript (Stewart and Beckenbach 2009). Therefore, in theory, it is possible to generate cDNA molecules containing multiple genes from unprocessed primary transcripts. Sequence reads generated from these partially processed primary transcripts would resemble polycistronic transcripts. Software programs designed to assemble high-throughput sequence data can also generate chimeric contigs from independent transcript sequences with matching nucleotide sequences at the ends (Mundry et al. 2012). Therefore, the tricistronic atp8-atp6-cox3 transcript detected in M. vitrata could possibly be an assembly artifact, but further studies are needed to verify this phenomenon. The detection of polyadenylated transcripts for all mitochondrial PCGs except for atp8 in H. zea may indicate either differences in mitochondrial transcript processing mechanisms compared with other lepidopteran species or erroneous application of models predicted for other insects (e.g., Drosophila) without actual experimental evidence. Observation of coverage map patterns representing single transcripts for all PCGs except atp8 and atp6 may indicate that atp8-atp6 transcript was the only mature multicistronic mitochondrial PCG in H. zea. This observation was consistent with the identification of polyadenylated transcripts for all H. zea mitochondrial PCGs except for atp8 in this study.

Differences in processing of the cob transcript observed between 2-h old eggs and 48-h embryos (Mitochondrial Transcripts) may be due to the transcriptional activity differences between the two developmental stages. In H. zea eggs incubated at 25°C, the egg nucleus complete the first and second maturation divisions within 90 min. after oviposition and cleavage nuclei in 16-cell stage were observed at 3 h (Presser and Rutschky 1957). Therefore, 2-h-old fertilized eggs of H. zea represent pre-embryonic stages that were either transcriptionally inactive or beginning to be active. Visual differentiation of unfertilized and fertilized eggs is not possible at very early stages and the sequence reads obtained from 2-h-old eggs could have come from a mixture of both fertilized and unfertilized eggs. All nuclear gene transcripts in eggs are maternally deposited during oocyte maturation. Although mitochondrial gene expression in vertebrate oocytes starts before replication of mitochondria (Ammini and Hauswirth 1999) no information is available on the mitochondrial gene transcription during insect embryonic development. In reference assembly of RNA-Seq reads, we detected mapped sequence reads without polyA+ tails extending past the 3′-ends all PCGs (i.e., reads originated from partially processed primary transcripts) and polyA+ containing reads at the 3′-ends of all PCGs except atp8 (both 2-h eggs and 48-h embryos) and cob (2-h eggs) indicating that both early and late embryos contain a mixture of mature and primary transcripts at various stages of processing. However, determination of developmental stage at which mitochondrial gene transcription activated was not possible with the data obtained in this study.

Transcript mapping analysis also revealed differences in transcription termination in some H. zea PCGs compared with published data for other lepidopteran species. Annotation of most mitigenomes has been carried out not by analyzing actual transcripts, but by using annotations from a mitogenome of a closely related species. For example, stop codon of the nad3 gene was considered the native TAG of all other lepidopteran species. Overlap of two-nucleotides between nad3 and trnA would require that the mitochondrial transcript processing machinery excise the trnA between the T and AG of the native stop codon, leaving the T at the 3′-end of the nad3. Because PCGs and tRNAs do not remain as bicistronic transcripts, it is only logical to assume that all species with overlapping nad3 and trnA would be processed similarly and the annotations in those mitogenomes were simply an overlook of the transcript processing requirements. Similarly, novel transcription end sites we described here for some PCGs may not be novel after all; detailed transcript analysis may reveal similar transcription termination or polyadenylation sites in other closely related lepidopteran species. It is also possible that the pattern of transcript processing observed developing embryos may be different from other life stages. Therefore, further studies will be necessary to determine the patterns of mitochondrial transcripts in other life stages. RNA-Seq libraries used in this analysis were not suitable to determine actual transcription start sites or 5′-ends of mature transcripts. Sanger dideoxy sequencing of a 5′-RACE cDNA library would be the most suitable method to identify the 5′-ends of the mitochondrial genes.

Supplementary Data

Supplementary data are available at Journal of Insect Science online.