Complete mitochondrial genome sequence of the “copper moss” Mielichhoferia elongata reveals independent nad7 gene functionality loss

Sample collection and DNA isolation

The M. elongata samples were collected from July 12 to 17, 2011 in the area near Mus-Khaya Peak (62°31′–36′N, 140°56′–141°07′E), Republic of Sakha (Yakutia) and deposited in MHA, the Herbarium of the Main Botanical Garden Russian Academy of Science, Moscow (Ignatova et al., 2011). This moss was originally identified in a cited paper as M. mielichhoferiana. However, a subsequent analysis of the nuclear rDNA 5.8S-ITS 2 region attributed this plant to that of morphologically hardly distinguishable M. elongata (Fig. S1). Rocks in the area are especially rich in MnS, with other heavy metals (Pb, Sn, As, Zn, Ag, etc., usually as sulfides) present in high concentrations. Consequently, many brooks have very acidic water and sulfur deposits along them. Siderite (iron carbonate) forms red outcrops (“Red rocks”) rich in iron and is always enriched with other heavy metals. When the outcrops are dry, M. elongata is the only moss that grows on this substrate or at least it is the only particularly abundant one.

DNA was extracted from specimens in the herbarium collection that had been gathered with a minimal soil amount and dried using ordinary herbarium techniques (in a paper envelope, under a tent, in the shade for several days until dry), and then stored in the herbarium at room temperature. A Nucleospin Plant DNA Kit (Macherey Nagel, Düren, Germany) was used for total DNA extraction from whole shoots of plants according to the manufacturers’ protocol. A yield of about 2 μg DNA was obtained according to measurements determined with a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA).

Library preparation and sequencing

A 500 ng sample of genomic DNA was fragmented using a Covaris S220 sonicator (Covaris, Woburn, MA, USA) and a library was prepared using TruSeq DNA sample preparation kit (Illumina, Mountain View, CA, USA). The concentration of the prepared library was measured with the Qubit fluorometer (Invitrogen, Carlsbad, CA, USA) and qPCR and fragment length distribution was determined with Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). The library was diluted to 10 pM and used for cluster generation on a cBot instrument with TruSeq PE Cluster Kit v3 reagents (Illumina, Mountain View, CA, USA). Sequencing was performed on a HiSeq2000 sequencer with read length of 101 from both ends of the fragments. About six million read pairs were obtained.

Mitogenome assembly and annotation

Raw sequencing reads were preprocessed with Trimmomatic software (Bolger, Lohse & Usadel, 2014) to remove adapters and low-quality data from further analysis. The whole genome assembly was then accomplished using the Spades assembler (Bankevich et al., 2012). A Blast database was generated from the assembled contigs, and a Blast search was performed against the Physcomitrella patens MG sequence (Terasawa et al., 2007) using the standalone NCBI BLAST-2.2.29+ (Altschul et al., 1990). The longest hit was the M. elongata complete MG. Iterative mapping was carried out using Geneious R10 software (https://www.geneious.com; Kearse et al., 2012) to verify the assembled genome. The resulting sequence had almost 100× coverage depth. The correctness of the genome boundaries was verified by PCR amplification followed by Sanger sequencing. Initial reads mapping to the genome sequence with Bowtie 2 (Lingmead et al., 2009) was applied as an additional genome structure verification step.

Genome annotation based on sequence similarity was performed using Geneious software. The MG sequence of Bartramia pomiformis which gave a maximum score in a BLAST search against a M. elongata MG query was applied as a reference. The annotated genome sequence was submitted to GenBank (accession number: MF417767). A circular genome map was drawn using the CGView Server (Grant & Stothard, 2008; http://stothard.afns.ualberta.ca/cgview_server).

SSR analysis

Simple sequence repeats (SSRs) were detected and located in the MG of M. elongata using GMATo v1.2 software (Wang, Lu & Luo, 2013).

Phylogenomic analysis

Phylogenetic reconstruction was conducted by selecting only functional protein-coding sequences (CDS) present in MGs of all bryophytes under investigation. A total of 33 of these CDS are known, including atp1, atp4, atp6, atp8, atp9, cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad9, rpL2, rpL5, rpL6, rpL16, rps1, rps2, rps4, rps7, rps11, rps12, rps13, rps14, rps19, sdh3, sdh4, and tatC. These were extracted from the MG sequences of 39 mosses, liverwort Treubia lacunose available in GenBank (http://www.ncbi.nlm.nih.gov), and the M. elongata sequenced in this work. The GenBank files were imported into Geneious R10 and merged to export a fasta dataset file. All sequences from this dataset were aligned using the default option implemented in MAFFT (Katoh & Standley, 2013). The final alignment was adjusted manually in BioEdit 7.2.5 (Hall, 1999).

Phylogenetic reconstruction was performed using the Bayesian method with the program MrBayes v3.2.6 (Ronquist et al., 2012). For Bayesian analyses, we used a parallel MPI version of MrBayes (Altekar et al., 2004). Two simultaneous runs of Metropolis Coupled Markov Chain Monte Carlo (MC3), both with one cold and seven heated chains were performed for 10 million generations. Two starting trees were chosen randomly. The general time reversible evolutionary model (GTR+I+G) with four rate categories was used. Posterior probabilities (PP) for trees and parameters were saved every 1,000 generations and parameters for each data partition were sampled independently from each other; the first 25% of the trees was discarded in each run. Bayesian PPs were used as branch support values.

Structure of the M. elongata mitogenome

The MG of M. elongata is 100,342 bp in length and has a typical circular structure (Fig. 1). The nucleotide composition of this genome has a GC content of 39.8%. The MG of M. elongata contains 66 genes including genes for three rRNAs (rrn18, rrn26, and rrn5), 24 tRNAs, and 39 conserved mitochondrial proteins (15 ribosomal proteins, four ccm proteins, eight nicotinamide adenine dinucleotide dehydrogenase subunits, five ATPase subunits, two succinate dehydrogenase subunits, one apocytochrome b, three cytochrome oxidase subunits, and one twin-arginine translocation complex subunit). Besides the functional genes, a single pseudogene, nad7, resides in the genome (Table 1).

Structure of nad7 gene in bryophytes

The lack of a functional gene copy of the nad7 gene has been reported previously in the MG of hornworts and the majority of liverworts (Groth-Malonek et al., 2007; Li et al., 2009; Xue et al., 2010). Evolution and losses of the functionality of the gene copies within mosses also deserve special attention and scrutiny. Pseudogenization of the nad7 gene is currently described for Tetraphis pellucida and Buxbaumia aphylla (Bell et al., 2014; Liu, Medina & Goffinet, 2014), whereas all other sequenced bryophyte MGs have a functional gene, that consists of three exons separated by two introns. The only known exception is the nad7 locus structure in MG of Hypnum imponens (NC 024516), its functional gene consists of only two exons and one intron sequences. The intron 2 of the gene was lost and exons 2 and 3 were merged together in one exon sequence. The low conservation of the pseudogene sequences has created difficulties in constructing a reliable nucleotide alignment and unambiguously judging whether exons 2 and 3 are completely deleted in either these chondriomes or whether some exon remnants are still preserved. We performed a Tblastn search of these exons amino acid sequence of the nad7 gene in T. pellucida and B. aphylla and confirmed the absence of exon 2 in B. aphylla and exon 3 in both species. The same finding is evident from Fig. S2 with the alignment of nad7 from B. aphylla, T. pellucida, M. elongata, and six other moss species. It agrees with the earlier data provided by Bell et al. (2014) on the structure of the T. pellucida MG. In addition, B. aphylla and T. pellucida pseudogenes have deletions in the sequences of the first gene exon, although at different locations. The main difference in the nad7 pseudogene primary structure in these bryophytes is two deletions in the sequence of exon 2 in B. aphylla, whereas T. pellucida has an intact exon 2 sequence. By contrast, the nad7 pseudogene of M. elongata completely lacks the second exon and has intact exon 1 sequence and exon 3 with frame shift mutation as a result of 2 bp insertion located at 190 bp from 5′ end of the exon (Fig. 2).

SSR analysis of the M. elongata mitochondrial genome

Following more stringent criteria (Zhao, Zhu & Liu, 2016) of perfect SSR locus identification (minimal number of repeating units ≥10 for mononucleotides, ≥5 for dinucleotides, ≥4 for trinucleotides, and ≥3 for tetra-, penta- and hexanucleotides) 73 SSR loci were identified in the MG of M. elongata (Table 2; Fig. 3). Most microsatellites refer to mono- and dinucleotides classes (35 and 28 loci, respectively). Trinucleotides are the least frequent SSRs group in the genome (one locus). No hexanucleotide microsatellite repeats occur in the genome. Among all the SSRs, 87.67% are composed only of A/T bases. The total length of the SSR loci is 852 bp, which comprises approximately 0.85% of the genome length.

Phylogenetic analysis

The alignment of 33 mitochondrial protein CDS of 40 moss taxa and hepatic Treubia lacunosa (Haplomitriopsida, Treubiidae, Treubiales, Treubiaceae) consists of 24,827 positions. The Bayesian phylogenetic tree inferred from this data with the hepatic T. lacunosa as an outgroup is shown in Fig. 4. Most nodes of the tree have very high PP supports. Two exceptions are two nodes among the Orthotrichaceae.