Mitochondrial genomes organization in alloplasmic lines of sunflower (Helianthus annuus L.) with various types of cytoplasmic male sterility

Plant material

Fertile line HA89 and isonuclear CMS lines—PET1 and PET2 of sunflower were obtained from the genetic collection of the N. I. Vavilov Institute of Plant Genetic Resources (VIR, Saint-Petersburg, Russia). The lines had the same nuclear genome (HA89), but they differed in chloroplast and mitochondrial genomes, inherited from their wild ancestors. The CMS sources were initially obtained by the interspecific hybridization of domesticated sunflower (H. annuus L.) with H. petiolaris Nutt (Leclercq, 1969; Whelan, 1980).

Mitochondrial DNA extraction, genome library construction and NGS sequencing

We extracted the organelle fraction with a reduced amount of nuclear DNA from leaves of 14-day sunflower seedlings, following the protocol of Makarenko et al. (2016). For every line, we used the same quantity of leaf tissue from five plants. The DNA isolation was performed with PhytoSorb kit (Syntol, Moscow, Russia), according to the manufacturer’s protocol. The NGS libraries preparations were made using one ng of DNA and Nextera XT DNA Library Prep Kit (Illumina, Mountain View, CA, USA), following the sample preparation protocol by Illumina. For the qualitative control of libraries, Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA) was used. The libraries quantitation was performed with the Qubit fluorimeter (Invitrogen, Carlsbad, CA, USA) and qPCR. Libraries for NGS sequencing were diluted up to the concentration of eight pM. Libraries were sequenced on different sequencing platforms. Fertile line and PET1 NGS libraries were sequenced with NextSeq 500 sequencer using High Output v2 kit (Illumina, Carlsbad, CA, USA). A total number of 13,240,057 150-bp paired reads were generated for fertile line and 14,758,067 reads—for PET1 line. PET2 library was sequenced with HiSeq2000 and MiSeq platforms using TruSeq SBS Kit v3-HS and MiSeq Reagent Kit v2 500-cycles (Illumina, San Diego, CA, USA). A total number of 4,471,774 125-bp and 4,931,318 250-bp paired reads were generated for the PET2 line.

Analysis of sequencing data

Quality of reads was determined with Fast QC. Trimming of adapter-derived and low-quality (Q-score below 25) reads was performed with Trimmomatic software (Bolger, Lohse & Usadel, 2014). Using the Bowtie2 tool v 2.3.3 (Langmead & Salzberg, 2012), sequencing reads were aligned to the reference sequence from NCBI databank (NC_023337.1). The Bowtie 2 alignments were done only for concordant paired reads (–no-mixed, --no-discordant options). Variant calling was made with samtools/bcftools software (Li, 2011) and manually revised using the IGV tool (Thorvaldsdóttir, Robinson & Mesirov, 2013). De novo assembly was performed with SPAdes Genome Assembler v 3.10.1 (Nurk et al., 2013) with different K values equal to 75, 85, 95, 127, and read coverage cutoff value equal to 30.0 (--cov-cutoff option). The potential ORFs were identified using ORFfinder. The graphical genome map was generated using the OGDRAW tool (Lohse, Drechsel & Bock, 2007). Transmembrane domains were predicted using the TMHMM Server v.2.0 (available online: http://www.cbs.dtu.dk/services/TMHMM-2.0/).

Validation of genome assembly—PCR and Sanger sequencing

The contigs obtained in de novo assembly were aligned to the reference sunflower mitochondrial genome (NC_023337.1) using BLAST. Validation of discovered rearrangements was made by PCR analysis and Sanger sequencing. PCR reactions were performed with LongAmp Taq PCR Kit (New England Biolabs, Ipswich, MA, USA) for reactions with expected amplicons more than 1.5 kb, and with Tersus Plus PCR kit (Evrogen, Moscow, Russia) for other reactions, including Sanger sequencing. For 28–29 cycles of PCR, we used 0.4 uM of primers (Table 1) and one ng of extracted DNA. The direct sequencing of purified amplicons was performed using the BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific, Waltham, MA, USA) and ABI Prism 3130xl Genetic Analyser (Applied Biosystems, Foster City, CA, USA).

RNA extraction and qRT-PCR

Total RNA from the leaves of five samples of each line was extracted with guanidinium thiocyanate-phenol-chloroform reagent kit—ExtractRNA (Evrogen, Moscow, Russia). RNA quality and concentration were measured using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and the Qubit fluorimeter (Invitrogen, Carlsbad, CA, USA). According to the manufacturer’s instruction, 0.5 μg of total RNA was treated with DNAse I (Thermo Fisher Scientific, Waltham, MA, USA). First-strand cDNA was synthesized using MMLV RT kit (Evrogen, Moscow, Russia) and specific primers. The qPCR was performed with designed primers (Table 1) and PCR kit with EvaGreen dye (Syntol, Russia) on Rotor-Gene 6000 (Corbett Research, Mortlake, NSW, Australia).

Organization of fertile line and PET1 mtDNA

De novo assembly of the fertile line and PET1 mitochondrial genomes revealed between nine and 12 large contigs (10–115 kb long), depending on K value. The optimum K value was 95. The obtained large contigs covered up to 95% of the reference genome. The remaining 5% of the mitochondrial genome correspond to repeats, therefore they were the breakpoints of contigs formation. Among numerous repeats in the mitochondrial genome, only one large (12,933 bp long) and six small (203–728 bp long) repeats played a crucial role in genome assembly and prevented single scaffold formation. Four regions (415–1,192 bp long) with 99% chloroplast DNA identity also affected the mitochondrial genome assembly. Eventually, manual assembly based on predominant contigs obtained by SPAdes supplemented by analysis of reads alignment by Bowtie2 and validation of controversial regions performed by PCR analysis and Sanger sequencing, allowed summary sequencing of data in completed mitochondrial genomes. Circular mtDNA of HA89 and PET1 lines are presented in Figs. 1 and 2.

The mtDNA comparative analysis of sunflower fertile lines HA412 (NCBI accession NC_023337.1) and HA89 revealed two single nucleotide thymidine insertions: in positions 35690–35691 and 129368–129369 of NC_023337. Two SNP in the noncoding part of 301 kbp genome is a negligible difference. We did not amend the HA89 mitochondrion sequence in the NCBI GenBank and, in this paper, used the same positions of mtDNA for HA412 and HA89 lines for simplicity.

The PET1 mitochondrial genome had structural rearrangements as well as polymorphic sites compared to the fertile lines’ (HA412/HA89) mtDNA. Previously, using the restriction analysis and Sanger sequencing, large structural variations of mtDNA associated with the sterility of plants were detected in the PET1 CMS type of sunflower—a 11 kb inversion and a five kb insertion (Kohler et al., 1991). The results of the current study not only confirmed the presence of these reorganizations in the HA89 PET1 mitochondrial genome but also allowed detection of more precise genome changes: 11,852 bp inversion, 4,732 bp insertion, 451 bp deletion. The revealed insertion was 98% identical to the PET1 insertion that can be obtained from GenBank NCBI (accession Z23137.1). Comparing with the fertile line genome, we have identified 18 variants in HA89 (PET1) mtDNA. Among the nucleotide variations, eight were localized in SSR loci, two deletions (single and dinucleotide) and seven SNP, including one transition and six transversions, were predominantly located in noncoding regions (Table 2). The exceptions were nonsynonymous mutations in orf777 (Asp251Glu), nad6 (Ser232Tyr), rpl16 (Lys32Gln). The HA89 (PET1) complete mitochondrial genome sequence has been deposited in the NCBI databank (accession MG735191).

Organization of PET2 mtDNA

We assembled the PET2 mitochondrial genome following the described procedure that was used for identification of HA89 fertile and PET1 complete mtDNA sequences. The complete mitochondrion of HA89 (PET2) is 316,582 bp (Fig. 3) and, in comparison with the HA89 fertile line, contained large-scale reorganizations of the mtDNA structure as well as minor changes represented by variable sites. Among significant rearrangements, two translocations, two deletions and two insertions were determined.

Even in a single plant cell, the mitochondrial genome is represented by several DNA molecules with various structure (Sloan et al., 2012). The so-called “master chromosome”—the single mtDNA molecule is a rare type of mitochondrial genome organization (Gualberto et al., 2014). More often, a mitochondrion includes a set of sub-genomic forms (Yang et al., 2015). Because sub-genomes could form differing master chromosomes, the statement of translocations in plant mitochondrial genomes is equivocal. To compare complete mitochondrial genomes of HA89 (PET2) and fertile HA89 lines, two translocations of approximately 27.5 kb and 106.5 kb (positions 37,112–64,614 bp and 194,439–300,945 bp in HA89 fertile line mtDNA) could be established. Using specific PCR primers (Table 1) we could demonstrate that the mtDNA in sunflower can also form a 154-kb sub-genomic circle molecule (positions 36393–190650). In the sunflower genome there is a repeat region of 722 bp (36393–37114 = 189929–190650 positions of the fertile line) with 100% similarity, which makes the cyclization of sub-genome circular molecule possible. In the case of HA89 (PET2) this sub-genome circle molecule has the wrong insertion point in the genome as compared with the fertile analog and thus results in translocations appearance.

Deletion of 711 bp (35682–36393 positions in fertile line mtDNA) resulted in the absence of orf777 in the HA89 (PET2) mitochondrial genome. The orf777 codes for a putative protein with unknown function, which shows no similarities with other mitochondrial proteins. So its elimination can hardly be the molecular reason for CMS development. The other deletion of 3,780 nucleotides (190659–194439 positions) affects only noncoding sequences. Notably, the deletions are associated with sub-genome integration regions (36393–37114 and 189929–190650 positions). So there is a possibility, that deletions in these regions could impair master chromosome assembly, which results in the translocations formations, as described above.

More significant mtDNA reorganizations are two revealed insertions—5,050 bp and 15,885 bp. Among them, the 5,050 bp insertion is most likely not involved in the origin of the CMS phenotype. This insertion was found in the intergenic region atp6-cox2 (275230–275231 positions of fertile line mtDNA), and it does not lead to the formation of new ORFs directly in the place of insertion into the mitochondrial genome. Moreover, there are no large (more than 300 nucleotides) ORFs within the sequence of 5,050 bp insertion. The exception is orf645 putatively translated into a polypeptide of 215 amino acids. Twenty-three amino acids at the N-terminus of this 215 amino acid protein were similar to N-terminus of ribosomal protein S3. We determined transcripts of orf645 in PET2 CMS line by qPCR, using specific primers (Table 1). However, mRNA of orf645 was absent in the fertile line and PET1 CMS line. Most often the molecular reason for CMS phenotype development is the emergence of chimeric ORFs with transmembrane domains, such as ATPase complex subunits, respiratory-related proteins, etc. (Gillman, Bentolila & Hanson, 2007; Yang, Huai & Zhang, 2009). Consequently, even if orf645 is translated in vivo, its role in male sterility development most likely is negligible, as there are no similarities with respiratory/ATP synthesis-related proteins.

The 15,885 bp insertion in the intergenic region nad4L-ccmFc is more complicated than the 5,050 bp insertion and includes different ORFs. First of all, it should be highlighted that most of the insertion (9,482 bp) has 100% similarity to another PET2 mtDNA region (126260–135741 positions). The repeating part of the insertion could be divided into two parts—6,097 bp (35686–41782 positions of PET2 mitochondrion) are common (99–100% identity) for Helianthus mitochondrion: 269147–275243, 273418–279514 and 126260–132356 positions of fertile line, PET1 and PET2 CMS lines, respectively. Such a repeat predominantly consists of noncoding sequences, except atp6 gene. The other 3,385 bp (41783–45167 positions of PET2 mitochondrion) have 100% similarity to the part of 5,050 bp insertion (132357–135741 positions of PET2 mitochondrion). However, this part of the insertion does not contain coding sequence. The next part of the 15,885 bp insertion counts 4,849 nucleotides which are unique for the PET2 mitochondrion (45168–50016 positions). The orf2565 is encoded in this part of insertion. The last 1,554 bp of the 15,886 bp insertion (50017–51570 positions of PET2 mitochondrion) have complex origin. A total of 390 nucleotides of this insertion (50017–50406 positions of PET2 mitochondrion) are similar to 114179–114568 positions of fertile line mtDNA, the next 271 bp (50406–50676 positions) are unique, then 500 bp (50677–51177 positions) are complement to 114587–115087 positions of the HA89 fertile mitochondrion and the rest 393 nucleotides (50677–51570) are also unique. From a functional point this region presents duplication of atp9 gene combined with 271 bp insertion and deletion of 12 bp, which resulted in two new ORFs formation—orf228 and orf285 (Reddemann & Horn, 2018). So the 15,885 bp insertion, especially its last 1,554 nucleotides, is the most important mitochondrial genome rearrangement, probably associated with PET2 CMS phenotype. It is also notable that 15,885 bp insertion has the same region in mtDNA of PET2 as the other reorganizations—711 and 3,780 bp deletions and 106.5 kb translocation.

Summarizing the functional changes data obtained, the 15,885 bp insertion consists of four coding sequences: duplicated atp6 gene and three new ORFs—orf2565, orf228 and orf285. Atp6 chimeric ORFs or new ORFs co-transcribed with atp6 are the quite common causes of CMS development in different plant species (Kim, Kang & Kim, 2007; Jing et al., 2012; Tan et al., 2018). So we proposed that atp6 gene and its colocalized area are of particular interest as the candidate sequence for CMS phenotype development. The annotation of H. annuus L. mitochondrion (NC_023337.1) presents the atp6 gene coding mRNA which consists of 906 nucleotides. But the fact is that in the 5′ adjoining sequence (150 bp) to atp6 start codon there is the other one initiating codon (ATG). Transcription from this codon results in new mRNA counting 1,056 nucleotides. The elongated transcript could be translated in putative ATP6 protein with additional 50 amino acids, wherein 37 of 50 amino acids are identical with N-terminus of coxI. The similar extension of protein was discovered in another CMS type of sunflower—ANT1 (Spassova et al., 1994), but the described protein had additional 87 aa at the C-terminus of the ATP6 protein. So we assumed that extended transcript, which could be produced in atp6 duplicated region, may cause CMS phenotype. To verify the hypothesis, the expression levels of atp6 and 5′ elongated atp6 transcripts were analyzed using qPCR with the same reverse primer but different forward primers (Table 1). The elongated atp6 transcript has been expressed at the same level as the atp6 gene (Gene ID: 18250997) both in PET1, PET2 CMS and fertile lines. So the atp6 transcript counting 1,056 nucleotides is the normal one, and there is a mistake in its annotation in NC_023337.1. Moreover, the atp9 gene also (Gene ID: 18250970) has wrong transcript annotation—it counts only 261 bp instead of proper one with 300 bp. Using qPCR and specific primers (Table 1) we established the same expression level for 261 bp and 300 bp transcripts in all three studied lines. Notable that the relative expression level (ΔCt) of atp6 gene to atp1 gene in PET2 CMS line had no significant difference as compared with fertile and PET1 CMS analogs, despite atp6 duplication in PET2 mitochondrion. Thereby atp6 duplication in PET2 mitochondrion does not involve in CMS appearance.

All three revealed ORFS (orf2565, orf228 and orf285) had demonstrated transcription activity in PET2 CMS line but not in the fertile or PET1 CMS analogs. The orf2565 translates to the putative polypeptide of 855 amino acids. Homology search in the NCBI database using BLAST pointed to similarity to DNA polymerase (type B). The impact of orf2565 on CMS phenotype development is quite doubtful, taking into account that the polypeptide has no transmembrane domains. The orf228 encodes a polypeptide of 76 aa, of which 75 are identical to C-terminus of ATPase subunit 9, as well as a new start codon (AUG) is formed due to the 271 bp insertion. Naturally, ATPase subunit 9 has two transmembrane domains—near N- and C-terminus, in orf228 there are also two predicted transmembrane domains (Reddemann & Horn, 2018). However, it should be noted that the N-terminus transmembrane domain of orf228 encoded polypeptide lacks four amino acids as compared with atp9 ones. This difference, in turn, could affect polypeptide (orf228) interaction with mitochondrial membrane and so results in mitochondrial membrane potential changes. The orf285 encodes a polypeptide of 95 aa, of which 18 are complement to N-terminus of ATPase subunit 9. In total orf228 and orf285 polypeptides share 93 of 99 amino acids of the atp9 protein. The polypeptide encoded by orf285 also has the transmembrane domain which formed by predominantly 20–42 aa. Thus both ORFs (orf228, orf285) could be the main reason the development of PET2 CMS type.

The comparative analysis of PET2 mitochondrion sequence with complete mtDNA of fertile line revealed 83 polymorphic sites—14 SSR, 13 small indels (1–29 bp) and 56 SNP (Table 2). Among nucleotide variations only two (synonymous) SNP were in the coding sequence of genes rps3 and atp6. Additional analysis of distribution of variants can offer insight into functional properties and evolution of sunflower mtDNA (Triska et al., 2017). Interestingly, PET1 and PET2 share 10 polymorphic sites as compared with the fertile line. The obtained HA89 (PET2) mitochondrial genome sequence has been deposited to the NCBI GenBank (accession MG770607.2). It is important to note that the sets of primers used for identification of PET2 insertions may be used for designing molecular markers for this type of CMS in sunflower.