Revealing mitogenome-wide DNA methylation and RNA editing of three Ascomycotina fungi using SMRT sequencing

Journal Pre-proofs Revealing mitogenome-wide DNA methylation and RNA editing of three Ascomycotina fungi using SMRT sequencing Chaoxia Wang, Jianhua Feng, Yujiao Chen, Dongmei Li, Li Liu, Yuqian Wu, Shujun Zhang, Simiao Du, Yaozhou Zhang PII: DOI: Reference:

S1567-7249(18)30147-8 https://doi.org/10.1016/j.mito.2020.01.003 MITOCH 1439

To appear in:

Mitochondrion

Received Date: Revised Date: Accepted Date:

12 June 2018 25 December 2019 3 January 2020

Please cite this article as: Wang, C., Feng, J., Chen, Y., Li, D., Liu, L., Wu, Y., Zhang, S., Du, S., Zhang, Y., Revealing mitogenome-wide DNA methylation and RNA editing of three Ascomycotina fungi using SMRT sequencing, Mitochondrion (2020), doi: https://doi.org/10.1016/j.mito.2020.01.003

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier B.V.

Revealing mitogenome-wide DNA methylation and RNA editing of three Ascomycotina fungi using SMRT sequencing Chaoxia Wang1,2*, Jianhua Feng2, Yujiao Chen6, Dongmei Li2, Li Liu6, Yuqian Wu7, Shujun Zhang3, Simiao Du5, Yaozhou Zhang3,4* 1. Management Center of Tianjin Modern Agricultural Science and Technology Innovation Base, Tianjin Academy of Agricultural Sciences, Tianjin, 300192, China. 2. Tianjin Lakeside Powergene Science Development Co. Ltd., Tianjin, 300309, China. 3. Research Center of Human Genome, Tianjin University, Tianjin, 300309, China. 4. Tianjin International Joint Academy of Biomedical, Tianjin, 300457, China. 5. Zheng Yuan Tang (Tianjin) Biotechnology Co. Ltd., Tianjin, 300457, China. 6. Tuke (Tianjing) Pharmaceutical Technology Co. Ltd., Tianjin, 300457, China. 7. Gui’an Precision Medicine Academy Co. Ltd., GuiZhou, 550029, China.

Abstract Beauveria bassiana, Cordyceps militaris and Ophiocordyceps sinensis (Ascomycotina) are traditional Chinese medicines. Here, mitogenomes of these three Ascomycotina fungi were sequenced and de-novo assembled using single-molecule real-time sequencing. The results showed that their complete mitogenomes were 31,258, 31,854 and 157,584 bp, respectively, with sequencing depth approximately 278,760×, 326,283× and 69,385×. Types of repeat sequences were mainly (AA)n, (AAT)n, (TA)n and (TATT)n. DNA methylation motifs were revealed in DNA modifications of these three fungi. We discovered new models of RNA editing through analysis of transcriptomes from B. bassiana and C. militaris. These data lay a solid foundation for further genetic and biological studies about these three fungi, especially for elucidating the mitogenome evolution and exploring the regulatory mechanism of adapting environment. Key words: Beauveria bassiana, Cordyceps militaris, Ophiocordyceps sinensis, mitogenome, DNA methylation, RNA editing

1. Introduction Mitochondria is an important organelle of metabolism and energy exchange, whose characteristics are specific gene components, genetic codon, replication mode and relative quick variation of sequences (Brown et al., 1979; Bruns and Szaro, 1992). For studying systematically mitogenomes of fungi, Fungal Mitochondrial Genome Project were raised and performed in 1997 (Paquin et al., 1997). Mitogenomes of some medicinal fungi were reported progressively, such as Oyster mushroom, Lucid ganoderma and Hericium erinaceus (Wang et al., 2008; Li et al., 2013 and 2015; Zhang et al., 2017). B. bassiana, C. militaris and O. sinensis are entomogenous fungi and traditional Chinese medicines, with some advantages on medicinal functions and industrial development (Kim and Yun, 2005; Li et al., 2006; Sung et al., 2007; Paterson, 2008; Khan et al., 2010; Zhang et al., 2012). The mitogenomes of B. bassiana Bb13 (29,961bp) and Bb147 (32,263bp) were obtained by PCR amplification (Xu et al., 2009; Ghikas, et al., 2010). Mitogenomes of 11 C. militaris strains were assembled through shotgun sequences using the Roche 454 platform, with lengths between 26,535 and 33,967 bp (Zheng et al., 2011; Sung, 2015; Zhang et al., 2015). The mitogenomes of O. sinensis 1229 and CCTCC AF 2017003 were reported through single-molecule real-time (SMRT) sequencing, with a size of 157,510 and 157,539 bp, respectively (Li et al., 2015; Kang et al., 2017). At present, few researches about their DNA methylations and transcriptomes are released. Corresponding author: Yaozhou Zhang, E-mail address: [email protected] Chaoxia Wang, E-mail address: [email protected]

DNA methylation is the most common epigenetic modification observed in the genomic DNA. Three DNA methylations, N6-methyl-adenine (6mA), N4-methyl-cytosine (4mC), and 5-methyl-cytosine (5mC), are present well-established as an epigenetic mark regulating gene expression, DNA replication, DNA repair, cell cycle progression and host-pathogen interaction (Jones, 2012; Roberts et al., 2015). M6A previously reported as a widespread DNA methylation in prokaryotes, but very limited in eukaryotes until the recent development of high-throughput sequencing (Collier et al., 2007; Fu et al., 2015; Koziol et al., 2016; Wu et al., 2016; Mondo et al., 2017). Among 16 diverse fungal genomes with SMRT sequencing technology, up to 2.8% of all adenines in early-diverging fungi were methylated into 6-methyladenine (6mA), far exceeding levels (0.05%-0.2%) observed in other eukaryotes and more derived fungi (Mondo et al., 2017). The methylation levels in fungi were greatly different. The methylation level in some fungi was very low or even zero. 6mA level in yeast genome was only 0.013% by LC-MS/MS, while 0.05% in SMRT-sequencing yeast genome (Kim et al., 2014; Liang et al., 2018). There was none methylated cytosine in -CCGG- sequence of yeast and Neurospora mitogenomes by restriction enzyme analysis (Groot and Kroon, 1979). DNA methylation level of Aspergillus flavus genome was negligible through bisulfite sequencing in conjunction with a biological replicate strategy (Liu et al., 2012). The levels of 6mA and 4-methylcytosine (4mC) in O. sinensis mitogenome were 0.017% and 0.13% by SMRT sequencing (Kang et al., 2017). Until now, the panoramic analysis of B. bassiana and C. militaris mitogenomic methylation has not been reported. Full-length transcripts are not released with assembled by SMRT sequencing, in favor of research on RNA modification, such as diversity of transcription isoforms, alternative splicing and poly-cistrons (Sharon et al., 2013; Gordon et al., 2015; Gao et al., 2018; Wang et al., 2018). PacBio full-length transcriptome profiling of insect mitogenome uncovered the isoform diversity of 12S rRNA and the polycistronic transcripts (ATP8/6 and Nad4/4L) (Gao et al., 2016). In this study, we sequenced mitogenomes of B. bassiana, C. militaris and O. sinensis and transcriptomes of the former two fungi using SMRT sequencing and performed comparative analysis with their mitogenomes reported. Moreover, DNA-methylation pattern and RNA-editing model were explored to broaden a better understanding of epigenetic and transcriptional modification in fungi mitogenome.

2. Materials and Methods 2.1 Fungal strains B. bassiana strain ZJ001 was provided by National Center for Silkworm Genetic Resources Preservation, Chinese Academy of Agricultural Sciences. The wild C. militaris collected from Changbai Mountains were cultured with silkworm pupae as previously described (Xiong et al. 2010). Through selection and domestication several years, C. militaris strain ‘Haining Yi Hao’ (HN001) was isolated by Zhejiang China gene Biomedical Co. Ltd, and identified by Institute of Microbiology, Chinese Academy of Sciences (Chen et al., 2005). Fruiting body of wild O. sinensis LJ001 was from Laji Mountain, Xining, Qinghai Province (Latitude 34.48◦N, Longitude 90.23◦E). The culture of B. bassiana and C. militaris was maintained at 23 ℃ on artificial medium. Their fruiting bodies were cultured for 90 days in our laboratory, then used to extract genomic DNA and total RNA. 2.2 DNA extraction and sequencing Genomic DNA was extracted by Universal Genomic DNA Kit purchased from CW

Biotechnology Co. Ltd (Beijing, China). The extracted DNA quality was checked by 0.7% agarose gel electrophoresis, assessed by Nanodrop1000 and quantified by using Qubit2.0 (Thermo Fisher Scientific). All extracted DNA was stored at -80 °C for further analysis. Large-insert PacBio library were prepared using the SMRTbell™ Template Prep Kit version 1.0 (Pacific Biosciences) according to manufacturer’s instructions. In brief, fungal DNA was sheared to 10~20 kb targeted size by using g-TUBEs (Covaris). The sheared genomic DNA was subject to DNA damage repair/end repair, blunt-end adaptor ligation followed by exonuclease digestion. The purified digestion products were loaded onto pre-cast 0.6% agarose for 7-20 kb size selection using the BluePippin Size Selection System (Sage Science), and the recovered size-selected library products were purified using 0.5× pre-washed PB AMPure beads (Beckman Coulter). The library concentration was determined with Qubit 2.0 Fluorometer (Life Technologies). Libraries sequenced using P6C4 polymerase and chemistry on PacBio RS II platform with 240 min movie times, at Tianjin Lakeside Powergene Science Development Co. Ltd (Tianjin, China). A total of 13 SMRT cells were used for each fungus. 2.3 RNA extraction and sequencing Total RNA was isolated using UNIQ-10 column Trizol total RNA extraction kit (Sangon Biotech), following the manufacturer's instructions, and followed by treatment of DNase I. mRNA was purified by poly T column separation and stored at -80 °C for further analysis. The Iso-Seq library was prepared according to the PacBio Isoform Sequencing protocol (Iso-Seq™). The RNA was reverse transcribed using the SMARTer® PCR cDNA Synthesis Kit, and PCR amplified using the KAPA HiFi PCR Kits. These cDNA products were purified using SMRTbell DNA Template Prep Kit 3.0 for library construction. Libraries were sequenced using P6C4 polymerase and chemistry on PacBio RS II platform with 240 min movie times at Tianjin Lakeside Powergene Science Development Co. Ltd. A total of 7 SMRT cells were used for each fungus. 2.4 Genome assembly and annotation De novo assembly of whole genome was carried out using the RS_HGAP_Assembly.3 protocol implemented in SMRT Analysis Portal v2.3.0 (Chin et al., 2013). All parameters were set default and exception that minimum subread length = 9,000; minimum seed read length = 11,000; genome size 35,000,000; target coverage = 30. The filtered reads were mapped back to the contigs using Blasr software (Chaisson and Tesler, 2012), and the contigs were polished by Quiver5 to generate a high-quality genome. Blast analysis against the assembled genome data, using the mitogenome of B. bassiana Bb147 (EU100742.1), C. militaris EFCC-C2 (KF432176.1) and O. sinensis (GenBank: KY622006.1) as the query, revealed one contig harboring mitochondrial genes, and this contig was complete mitochondrial DNA genome. The mitogenomes of these three fungi were first annotated automatically using the MFannot tool (http://megasun.bch.umontreal.ca/cgi-bin/mfannot/mfannotInterface.pl). The nucleotide sequences of PCGs were translated to protein sequences using the Mold, Protozoan and Coelenterate Mitochondrial code (transl_table D4). The MFannot program can predict protein-encoding genes, rRNAs, tRNAs, and intron types. And then, the tRNA genes were checked using tRNAscan-SE software (Lowe and Chan, 2016). The predictions were individually checked and then revised according to the annotation of reference strains. The mitochondrial genome map was generated with Circos software (Krzywinski et al., 2009). Analysis of comparative mitogenome were performed by Mauve 2.3.1 to illustrate positions of introns and major alignment gaps. (Darling et al., 2010).

2.5 Analysis of repetitive sequences Local BLASTn search of mitogenome against itself was performed using a cut-off e-value of −7 10 (Altschul et al., 1990 and 1997). Repetitive sequences were analyzed by Repeatmasker 4.0.6 (http://www.repeatmasker.org/) and REPuter (http://bibiserv.techfak.uni-bielefeld.de/reputer/, Kurtz et al., 2001). REPuter was applied to identify the forward, reverse, complement and palindromic sequences. 2.6 Methylation modification and motif analysis DNA modification detection and motif analysis were performed using the PacBio SMRT analysis platform (protocol version = 2.3.0.140936, method = RS Modification and Motif Analysis.1, http://www.pacb.com/devnet/code.html). Briefly, raw reads were filtered using SFilter software, to remove short reads and reads derived from sequencing adapters. Filtered reads were aligned to the reference genome using BLASR(v1) (Chaisson et al. 2012). Modified sites were then identified through kinetic analysis of the aligned DNA sequence data (Flusberg et al. 2010). Modified sites were then grouped into motifs using Motif Finder v1.2.1. These motifs represent the recognition sequences of MTase genes active in the genome (Clark et al. 2011). 2.7 Transcription and splicing analysis Blast analysis against B. bassiana and C. militaris mitogenomes, using their transcriptome data as the query, revealed some transcripts harboring mitochondrial genes. 2.8 Data availability The accession numbers for B. bassiana, C. militaris and O. sinensis mitogenomes assembled are MH211586, MH211587 and MH400233, respectively.

3. Results 3.1 General bioinformatics analysis for three mitochondrial genomes The raw data of B. bassiana, C. militaris and O. sinensis mitogenomes were 8.7, 10.4 and 10.9 G, with 680,010, 978,246 and 1,194,940 reads, respectively. Through de novo assembly with RS_HGAP, the final mitogenomes of three fungi were assembled into one contig with a size of 31,257, 31,854 and 157,499 bp, and an N50 of 19.6, 19.4 and 13.7 kb, respectively (Table 1). The average sequencing depth were approximately 278,760×, 326,283× and 69,385×, respectively. The bioinformatics analysis provided the general genome information (Table 1), including GC% content (27.0%, 26.9% and 30.2%, respectively), the proportion of protein-coding genes (57.2%, 47.7% and 68.9%, respectively), the proportion of introns (17.7%, 24.2% and 61.0%, respectively), and the proportion of intergenic regions (15.6%, 13.8% and 15.1%, respectively). Table 1 General information statistics of mitochondrial genome B. bassiana

C. militaris

O. sinensis

Number of bases (bp)

8,713,489,035

10,393,415,752

10,933,990,048

Number of reads

680,010

978,246

1,194,840

The length of N50 (bp)

19,603

19,317

13,671

The length of mean reads (bp)

12,813

10,624

9,151

The scores of mean reads

0.85

0.85

0.84

Mapped reads

26,223

63,568

77,478

Mapped read length of insert (bp)

5,208

3,966

3,759

Average length (bp)

31,258

31,854

157,584

Average coverage

278,760.29

326,282.91

69,385.15

G+C content (%)

27

27

30.2

Protein-coding genes (% of total

57.2

47.7

68.9

rRNAs + tRNAs (%)

27.2

34.5

16

Inter-genic regions (%)

15.6

13.8

15.1

Introns (%)

17.7

24.2

61

mtDNA)

3.2 Mitogenome annotation and structure characteristics The mitogenomes of these three fungi were annotated by MFannot tool, all of which include 14 PCGs and two ribosomal RNAs (large subunit ribosomal RNA rnl and small subunit ribosomal RNA rns) arranged in the same order (cox1-nad1-nad4-atp8-atp6-rns-cox3-nad6-rnl-nad2-nad3atp9-cox2-nad4L-nad5-cob), and varying amounts of tRNAs and open reading frames (ORFs) (Fig. 1). Furthermore, different numbers of introns existed in some PCGs (Supplementary Table S1).

Fig. 1. Circular map of the B. bassiana, C. militaris, O. sinensis mitogenomes. From outermost to innermost: ring 1 includes protein-coding genes (blue), rRNAs (grey), tRNAs (red) and ORFs (yellow). Rings 2 and 3 show DNA methylation sites in the forward and backward strand, respectively. Red lines indicate 4mC, blue lines show 6mA. Ring 4 and 5 show other DNA base modifications sites in the forward and backward strand, respectively. Ring 6 with grey background are GC contents in the sliding window is over 200 bp. The ribbons inside the circle connect the multi-collinearity of genomes with significant (e-value < 10-7) similarity.

Analysis of repetitive sequences in these three mitogenomes were performed by REPuter online. The longest repetitive sequences were all forward, the length of which were 45, 117 and

128 bp in B. bassiana, C. militaris and O. sinensis, respectively. The lengths of repetitive sequences in O. sinensis mitogenome had large differences and were mainly 18-35 bp (Fig. 2). The main length of repetitive sequences in B. bassiana and C. militaris mitogenomes was 15-17 bp. This result indicated a positive correlation between the length of repetitive sequences and the size of mitogenome. Types of repetitive sequences were mainly (AA)n, (AAT)n, (TA)n and (TATT)n by Repeatmasker (Table 2; Supplementary Table S2). (AA)n and (AAT)n in B. bassiana mitogenome preferred to be distributed in exon regions. In O. sinensis mitogenome, these four types were mainly found in intergenic regions and intron regions, less in exon regions. No obvious preference about the distributions of repetitive sequences was in C. militaris mitogenome.

Fig. 2. Box plot representation of length distribution of repetitive sequences in B. bassiana, C. militaris, and O. sinensis mitogenomes. The crosses represent the 95% confidence interval around the median. Table 2 Distributions of repetitive sequences Types

(AA)n

B. bassiana

C. militaris

O. sinensis

Num.

Location

Num.

Location

Num.

Location

4

nad2, nad4, atp6,

1

cox2-E1,

11

nad3, nad2-E1, intergenic region (2), orf481,

intergenic region

orf308, orf413, atp6-I2, cob-E6,

intergenic region

cox1-I6-orf322, cob-I3-orf415 (AAT)n

6

intergenic region and

5

atp8, intergenic region, nad1, cox1-E1

(2)a,

cob-I1-orf305,

9

nad5-I3-orf113, rnl-I5, rnl-I5-orf539, nad6-E2,

cox2-E1, rnl-I1, nad4,

intergenic region, cox1-I13, cox1-I5-orf354,

nad1

cox1-I4 (2),

cox2-E1, (TA)n

3

rnl-I1-orf462/rps3,

2

intergenic region (2) (TATT)n

0

nad4,

4

intergenic (3), rnl-E7,

10

rnl-E3, intergenic region (4), atp6-I1, nad1-I4

intergenic region 2

cox1-I1-orf314, nad1

(2), nad1-I1, cox1-I3, a

The number in the bracket represents the amounts of types.

3.3 Methylation analysis In B. bassiana mitogenome, a total of 584 modification sites were determined with an average coverage of ~3782 × . There were 264 modification sites located on the forward strand versus 320 located on the backward strand (Fig. 1; Fig. 3). Eleven 6mA (0.097% of all adenine)

and 22 4mC (0.471% of all cytosine) modification sites were identified in B. bassiana mitogenome. In C. militaris mitogenome, we identified 8 putative 6mA (0.068%), 14 4mC (0.377%) and further 476 rather unspecific “modified bases” where the type of modification motifs was not recognized by the software. There were 261 modification sites on the forward strand versus 237 located on the backward strand in C. militaris mitogenome. In O. sinensis mitogenome, a total of 4600 modification sites were determined at an average coverage of ~8061 × . There were 2294 modification sites located on the forward strand versus 2306 located on the backward strand, of which 178 were 4mC (0.68%) and 85 6mA (0.16%).

Fig. 3. Numbers of methylation sites in three mitogenomes

3.4 Sequence feature and distribution statistics of methylation sites The five upstream and downstream bases of methylation sites were proceeded in three mitogenomes (Fig. 4). The sequence feature of 6mA was TmAGGTA in B. bassiana and C. militaris mitogenomes, while D(G/A/T)mAG in O. sinensis mitogenome. The sequence feature of 4mC was all ANNNNmCNGTA in three mitogenomes (Fig. 3; Supplementary Table S3).

Fig. 4. The accumulated nucleotide composition of all potential 6mA and 4mC identified by SMRT and computation. Sequence logo is generated by WebLogo (http://weblogo.berkeley.edu/). The height represents the significance of enrichment of each nucleotide at that position. (A), (C) and (E) showed sequence feature of 6mA, while (B), (D) and (F) showed sequence feature of 4mC in B. bassiana, C. militaris and O. sinensis mitogenomes, respectively.

The distributions of these methylation sites were analyzed separately within coding regions of genes (CDS, exon regions) and noncoding regions (NRs, including intergenic regions, intron regions and tRNA) (Fig. 5). Two-thirds 6mA and 4mC were distributed in intergenic regions (between tRNA and nad1/nad6/cox3, or between tRNA and tRNA) and intron regions of different genes (e.g., nad1-2, nad5, cox1-2, cob and rnl) in B. bassiana and O. sinensis mitogenomes; Less one-third 6mA and 4mC were in the exon regions of PCGs, such as cob-E1 (Figs. 5A, 5B, 5E and 5F; Supplementary Table S6). In C. militaris mitogenome, the methylation sites in exon regions were slightly more than those located within NRs (Figs. 5C and 5D; Supplementary Table S4 and S5).

Fig. 5. Distributions of 6mA and 4mC modified sites in three mitogenomes. (A), (C) and (E) showed distributions of 6mA, while (B), (D) and (F) showed distributions of 4mC in B. bassiana, C. militaris and O. sinensis mitogenomes, respectively.

3.5 Transcription analysis The data of B. bassiana transcriptome was 16 Kb with 11 full-length transcripts (Supplementary Table S7). Aligned with B. bassiana mitogenome by blastn, four transcripts were aligned to rns ， of which the size were from 814 to 1280 bp and the depth 3 or 4. These four

transcripts could be classified into two transcription isoforms. The other seven transcripts aligned to rnl were grouped into three transcription isoforms, with the length of 259 to 1904 bp and depths of 1-5 (Fig. 6). RNA editing was found in two transcripts (c47637/f1p3/1546 and c5751/f1p1/1366), the splicing sites of which were 26844 and 28742 in B. bassiana mitogenome, and the splicing pattern was GT-AG (Fig. 6).

Fig. 6. Transcription isoforms and splice junctions of rnl in B. bassiana. There were three transcription isoforms of rnl. The splicing model was GT-AG.

The data of C. militaris transcriptome was 22 Kb with 18 full-length transcripts (Supplementary Table S7). Aligned with C. militaris mitogenome by blastn software, 10 transcripts were aligned to rns，of which the size was from 618 to 1811 bp, and depths between 2 and 30. These ten transcripts could be divided into four transcription isoforms. The other eight transcripts were aligned to rnl, with the length between 693 and 1537 bp, depths of 2-7, and five transcription isoforms (Fig. 7). Three splicing sites were 13850 and 15043, 15221 and 16461, and 20641 and 13195, with patterns of AA-TG, AA-CG and AA-AA, respectively.

Fig. 7. Transcription models and splice junctions of rnl in C. militaris. There were five transcription isoforms of rnl. The splicing models were AA-TG, AA-CG and AA-AA.

In B. bassiana mitogenome, the number of methylated sites in rns and rnl were separately 1 and 4, with depths of transcripts 13 and 23, which showed more methylation caused higher

transcription level (Table 3). On the contrary, the number of methylated sites in rns and rnl in C. militaris mitogenome were severally 2 and 7, with depths of transcripts 89 and 33, showing that methylation could suppress transcription level. Table 3 Numbers of methylation sites and depths of transcripts in rns and rnl rns Strains

rnl

Methylation

Transcript

Methylation

Transcript

num.

depth

num.

depth

B. bassiana

1

13

4

23

C. militaris

2

89

7

33

3.6 Structure variation Compared with other reported B. bassiana mitogenomes, some structural variations existed in B. bassiana ZJ001 mitogenome (Fig. 8). The size of ZJ001 mitogenome was 2442 bp more than that of K4 (KT201148.1) mitogenome, due to cox1 and cox2 in ZJ001 mitogenome with one intron of 1271 bp and 1120 bp, respectively. The size of ZJ001 mitogenome was 1314 bp more than that of E17 (KT201149) mitogenome, because of no intron in cox1 and the intron in cox2 harboring ORF281 in E17 mitogenome and ORF295 with 42 more bases in ZJ001 mitogenome. However, the size of ZJ001 mitogenome was 1005 bp less than that of Bb147 (EU100742.1) mitogenome, mainly on account of one intron in nad1 of Bb147 with 1075 bp.

Fig. 8. Alignments of four B. bassiana mitogenomes. Numbered scale bars indicate distance in base pairs, and vertical bars indicate sequence similarity. Size differences among four mitogenomes were mainly due to presence/absence of introns in certain genes and length variations at certain regions. Major alignment gaps (G1– G3) are marked. G1, the cox2 intron present in E17, ZJ001 and Bb147 but absent in K4; G2, the cox1 intron present in ZJ001 and Bb147 but absent in K4 and E17; G3, the nad1 intron only present in Bb147 but absent in K4, E17and ZJ001.

Some structural variations were found between C. militaris HN001 mitogenome and other reported C. militaris mitogenomes (Fig. 9). Compared with EFCC-C2 mitogenome, no intron with 1,214 bp in cox3 and one less intron in rnl were existed in HN001 mitogenome. The size of HN001 mitogenome was 5318 bp more than that of CM06 mitogenome, because of cox1、cox2 and cob all with one intron (1052, 1125 and 1258 bp, respectively), one more intron with 1193bp in rnl and one more orf108 found in HN001 mitogenome.

Fig. 9. Alignment of six C. militaris mitogenomes. Numbered scale bars indicate distance in base pairs, and vertical bars indicate sequence similarity. Size differences among six mitogenomes were mainly due to presence/absence of introns in certain genes and length variations at certain regions. Major alignment gaps (G1– G9) are marked. G1, the cox1 intron present in most strains except CM06; G2, orf108 present in CM09-9-24, CMB, HN001 and F02 but absent in CM06 and EFCC-C2; G3, the orf108-trnR intergenic region which is 366 bp present in CM09-9-24, CMB, HN001 and F02 but absent in CM06 and EFCC-C2; G4, the rnl-i4 present in most strains except CMB; G5, G6 and G7, separately rnl-i1, cox2 intron and cob intron present in CMB, HN001, F02 and EFCC-C2 but absent in CM06 and CM09-9-24; G8 and G9, separately cox3 intron and rnl-i3 both only present in EFCC-C2.

The size of O. sinensis LJ001 mitogenome was only 45 bp more than that of KY622006 mitogenome, mainly due to insertions and deletions of single base and one insertion with 21 bp in nad2-i1 from LJ001 mitogenome (Fig. 10; Supplementary Table S6). Insertions and deletions of single base preferred to be in intergenic regions and intron regions, leading to different ORFs in intron regions. The same location of ORF280 in cob-i4 from LJ001 mitogenome was ORF236 in KY622006 mitogenome, on account of one more cytosine at 147,693 in LJ001 mitogenome. The mutation of TAA to GAA in LJ001 mitogenome leaded to 1061 bp nad5 longer than that in KY622006 mitogenome. The structure of nad5 protein with more 397 amino acid residues in LJ001 mitogenome formed some pleated sheets and helixes, except 12 transmembrane helixes the same as that of nad5 protein in KY622006 mitogenome (Fig. 11).

Fig. 10. Alignments of two O. sinensis mitogenomes. Numbered scale bars indicate distance in base pairs, and vertical bars indicate sequence similarity.

Fig. 11. Three-dimensional space structures of nad5 protein. (A) Structure of nad5 protein in O. sinensis KY622006. (B) Structure of nad5 protein in O. sinensis LJ001. (C) Overlapping structures of nad5 protein.

4 Discussion DNA methylation of fungi was rarely reported in previous researches, especially 5-methylcytosine (Antequera et al., 1984; Magill and Magill, 1989). DNA methylations of 17 eukaryotic genomes were reported including 5 fungi, 5 plants and 7 animals (Zemach et al., 2010). Early unicellular animals and fungi reproduced asexually were more preferred to lose DNA methylations. Loss of this pathway had probably occurred early in fungal evolution and later in several plant and animal lineages. The levels of DNA methylation were negligible in some fungi, such as Saccharomyces cerevisiae (Proffitt et al., 1984), Schizosaccharomyces pombe (Capuano et al., 2014) and Aspergillus flavus (Liu et al., 2012). Such low levels of DNA methylation were due to diverse limited sensitivity and specificity of MSREs, bisulfite sequencing and HPLC-MS/MS (Dahl and Guldberg, 2003). Furthermore, levels of DNA methylation correlated with selected samples of the day on account of dynamical DNA methylation procedure (Wang et al., 2015). According to the kinetic analysis of DNA synthesis, SMRT has higher sensitivity (Kim et al., 2014; Liang et al., 2018). In this study, we revealed a mitogenome-wide epigenetic DNA modification pattern in B. bassiana, C. militaris and O. sinensis mitogenomes by SMRT sequencing, further exploring epigenetic modifications in fungal mitogenome. Mondo et al. (2017) firstly reported 6mA in 16 fungi and revealed 6mA level (0.048-0.21%) in Ascomycota and Basidiomycota by SMRT sequencing. Our study quantified 6mA modifications in B. basisana (0.097%), C. militaris (0.068%) and O. sinensis (0.16%) mitogenomes, which enriched DNA methylation research of Ascomycota. Moreover, 4mC levels were also presented in the three mitogenomes (0.47%, 0.38% and 0.68%, respectively) (Fig. 1 and Fig. 3). 6mA and 4mC levels in this study were much higher than those in O. sinensis mitogenomes reported by Kang et al. (2017). Most of these methylated sites were enriched at intergenic and intronic regions, and others were in exons, rnl, rns and tRNAs. The hypermethylation in promoters and introns, and hypomethylation in exons, might be closely related to the expression of PCGs or the binding affinity of transcription factors to mitogenome (Volpe, 2005). DNA methylation is a crucial component of epigenetic regulation that controls many important aspects of biological processes such as host-pathogen interactions, DNA damage and DNA repair (Tollefsbol, 2011). Mitochondria is lack of protection of histone proteins, without relative effective system of DNA repair (Byun et al., 2013). Therefore, mitochondria were easily suffered from reactive oxygen species (ROS), which was a decisive factor that changed DNA methylation (Valinluck et al., 2004). DNA methylation levels in B. bassiana and C. militaris mitogenomes were lower than that in O. sinensis mitogenome in our study. The reason was inferred that the former two fungi were collected from artificial culture and grew fast, while the wild O. sinensis sample grew slowly in high-altitude, low-temperature and low-oxygen

environments (Li et al., 2011; Yu et al., 2011). Mitochondria is crucial for responses to hypobaria, hypothermia, and hypoxia due to their central role in energy production and consumption (Chandel and Schumacker, 2000; Chitra and Boopathy, 2013). These hyper-methylations in O. sinensis may affect the enzymatic activity of oxidative phosphorylation and then the mitochondria respiratory rate; or they may influence mitochondria biogenesis. We supposed that hyper-methylations in O. sinensis mitogenome might be a genetic feature to adapt to the cold and low PO2 environment at high altitude and be closely related to the expression of electron transport subunits, and thereby modulate the enzymatic activity of OXPHOS and the mitochondria respiratory rate, and then change the ability to capture O2 and produce energy. In summary, this study uncovered potential advantages of PacBio SMRT in detecting DNA methylation and RNA editing. The main diversity between mitogenomes of B. bassiana, C. militaris and O. sinensis and those mitogenomes reported by others, was differences of introns probably caused by differences of strains and growing environment (Figs. 8 and 9). We believed that their methylome characterization would help us to completely understand the mitochondrial genomic function and need further researches.

References 1.

Antequera F, Tamame M, Villanueva JR, et al. DNA methylation in the fungi[J]. Journal of Biological Chemistry, 1984, 259(13): 8033-8036.

2.

Altschul,S.F. et al. Basic local alignment search tool. J. Mol. Biol., 1990, 215, 403-410.

3.

Altschul,S.F. et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res., 1997, 25, 3389-402.

4.

Brown WM, George M, Wilson AC. Rapid evolution of animal mitochondrial DNA[J]. Proceedings of the National Academy of Sciences, 1979, 76(4): 1967-1971.

5.

Bruns TD, Szaro TM. Rate and mode differences between nuclear and mitochondrial small-subunit rRNA genes in mushrooms[J]. Molecular Biology and Evolution, 1992, 9(5): 836-855.

6.

Byun HM, Panni T, Motta V, et al. Effects of airborne pollutants on mitochondrial DNA methylation[J].

7.

Capuano F, Mülleder M, Kok R, et al. Cytosine DNA methylation is found in Drosophila melanogaster but

Particle and fibre toxicology, 2013, 10(1): 18. absent in Saccharomyces cerevisiae, Schizosaccharomyces pombe, and other yeast species[J]. Analytical chemistry, 2014, 86(8): 3697-3702. 8.

Chaisson MJ, Tesler G. Mapping single molecule sequencing reads using basic local alignment with successive refinement (BLASR): application and theory[J]. BMC bioinformatics, 2012, 13(1): 238.

9.

Chandel NS, Schumacker PT. Cellular oxygen sensing by mitochondria: old questions, new insight[J]. Journal of applied physiology, 2000, 88(5): 1880-1889.

10.

Chen G, XU C, GONG C, et al. Pharmacology of cultivated Haining strain of silkworm Cordeceps militaris[J]. Chinese Journal of Applied and Environmental Biology, 2005, 11(4): 453.

11.

Chin CS, Alexander DH, Marks P, et al. Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data[J]. Nature methods, 2013, 10(6): 563.

12.

Chitra L, Boopathy R. Adaptability to hypobaric hypoxia is facilitated through mitochondrial bioenergetics: an in vivo study[J]. British journal of pharmacology, 2013, 169(5): 1035-1047.

13.

Clark TA, Murray IA, Morgan RD, et al. Characterization of DNA methyltransferase specificities using single-molecule, real-time DNA sequencing[J]. Nucleic acids research, 2011, 40(4): e29-e29.

14.

Dahl C, Guldberg P. DNA methylation analysis techniques[J]. Biogerontology, 2003, 4(4): 233-250.

15.

Darling A E, Mau B, Perna N T. progressiveMauve: multiple genome alignment with gene gain, loss and rearrangement[J]. PloS one, 2010, 5(6): e11147.

16.

Flusberg BA, Webster DR, Lee JH, et al. Direct detection of DNA methylation during single-molecule, real-time sequencing[J]. Nature methods, 2010, 7(6): 461-465.

17.

Fu Y, Luo GZ, Chen K, et al. N6-methyldeoxyadenosine marks active transcription start sites in Chlamydomonas. Cell, 2015, 161: 879-892.

18.

Gao S, Ren Y, Sun Y, et al. PacBio full-length transcriptome profiling of insect mitochondrial gene

19.

Gao S, Tian X, Chang H, et al. Two novel lncRNAs discovered in human mitochondrial DNA using PacBio

expression[J]. RNA biology, 2016, 13(9): 820-825. full-length transcriptome data[J]. Mitochondrion, 2018, 38: 41-47. 20.

Ghikas DV, Kouvelis VN, Typas MA. Phylogenetic and biogeographic implications inferred by mitochondrial intergenic region analyses and ITS1-5.8 S-ITS2 of the entomopathogenic fungi Beauveria bassiana and B. brongniartii[J]. BMC microbiology, 2010, 10(1): 174.

21.

Gordon SP, Tseng E, Salamov A, et al. Widespread polycistronic transcripts in fungi revealed by single-molecule mRNA sequencing[J]. PloS one, 2015, 10(7): e0132628.

22.

Groot GSP, Kroon AM. Mitochondrial DNA from various organisms does not contain internally methylated cytosine in-CCGG-sequences[J]. Biochimica et Biophysica Acta (BBA)-Nucleic Acids and Protein Synthesis, 1979, 564(2): 355-357.

23.

Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond[J]. Nature Reviews Genetics, 2012, 13: 484-492.

24.

Kang X, Hu L, Shen P, et al. SMRT Sequencing Revealed Mitogenome Characteristics and Mitogenome-Wide DNA Modification Pattern in Ophiocordyceps sinensis[J]. Frontiers in microbiology, 2017, 8: 1422.

25.

Khan MA, Tania M, Zhang D, et al. Cordyceps mushroom: a potent anticancer nutraceutical[J]. The Open Nutraceutical Journal, 2010, 3: 179-183.

26.

Kim HO, Yun JW. A comparative study on the production of exopolysaccharides between two entomopathogenic fungi Cordyceps militaris and Cordyceps sinensis in submerged mycelial cultures[J]. Journal of Applied Microbiology, 2005, 99(4): 728-738.

27.

Kim KE, Peluso P, Babayan P, et al. Long-read, whole-genome shotgun sequence data for five model organisms[J]. Scientific data, 2014, 1: 140045.

28.

Koziol MJ, Bradshaw CR, Allen GE, et al. Identification of methylated deoxyadenosines in vertebrates reveals diversity in DNA modifications. Nature Structural & Molecular Biology, 2016, 23: 24-30.

29.

Krzywinski M, Schein J, Birol I, et al. Circos: an information aesthetic for comparative genomics[J]. Genome research, 2009, 19(9): 1639-1645.

30.

Kurtz S, Choudhuri J V, Ohlebusch E, et al. REPuter: the manifold applications of repeat analysis on a genomic scale[J]. Nucleic acids research, 2001, 29(22): 4633-4642.

31.

Li J, Zhang J, Chen H, et al. Complete mitochondrial genome of the medicinal mushroom Ganoderma

32.

Li SP, Yang FQ, Tsim KWK. Quality control of Cordyceps sinensis, a valued traditional Chinese

lucidum[J]. PLoS One, 2013, 8(8): e72038. medicine[J]. Journal of Pharmaceutical and Biomedical Analysis, 2006, 41(5): 1571-1584. 33.

Li Y, Hu XD, Yang RH, et al. Complete mitochondrial genome of the medicinal fungus Ophiocordyceps sinensis[J]. Scientific reports, 2015, 5: 13892.

34.

Li Y, Wang XL, Jiao L, et al. A survey of the geographic distribution of Ophiocordyceps sinensis[J]. The Journal of Microbiology, 2011, 49(6): 913-919.

35.

Liang Z, Yu G, Liu J, et al. The N6-adenine methylation in yeast genome profiled by single-molecule technology[J]. Journal of Genetics and Genomics, 2018, doi: 10.1016/j.jgg.2018.03.003.

36.

Liu SY, Lin JQ, Wu HL, et al. Bisulfite sequencing reveals that Aspergillus flavus holds a hollow in DNA methylation[J]. PloS one, 2012, 7(1): e30349.

37.

Lowe TM, Chan PP. tRNAscan-SE on-line: search and contextual analysis of transfer RNA genes. Nucleic Acids Research, 2016, 44: W54-57.

38.

Magill JM, Magill CW. DNA methylation in fungi[J]. Developmental genetics, 1989, 10(2): 63-69.

39.

Mondo SJ, Dannebaum RO, Kuo RC, et al. Widespread adenine N6-methylation of active genes in fungi[J]. Nature genetics, 2017, 49(6): 964-968.

40.

Paquin B, Laforest M J, Forget L, et al. The fungal mitochondrial genome project: evolution of fungal mitochondrial genomes and their gene expression[J]. Current genetics, 1997, 31(5): 380-395.

41.

Paterson RRM. Cordyceps: A traditional Chinese medicine and another fungal therapeutic biofactory?[J]. Phytochemistry, 2008, 69(7): 1469-1495.

42.

Proffitt JH, Davie JR, Swinton D, et al. 5-Methylcytosine is not detectable in Saccharomyces cerevisiae DNA[J]. Molecular and cellular biology, 1984, 4(5): 985-988.

43.

Roberts RJ, Vincze T, Posfai J, et al. REBASE-a database for DNA restriction and modification: enzymes, genes and genomes[J]. Nucleic Acids Research, 2015, 43: 298-299.

44.

Sharon D, Tilgner H, Grubert F, et al. A single-molecule long-read survey of the human transcriptome[J]. Nature biotechnology, 2013, 31(11): 1009.

45.

Sung GH, Hywel-Jones NL, Sung JM, et al. Phylogenetic classification of Cordyceps and the clavicipitaceous fungi[J]. Studies in Mycology, 2007, 57: 5-59.

46.

Sung GH. Complete mitochondrial DNA genome of the medicinal mushroom Cordyceps militaris (Ascomycota, Cordycipitaceae)[J]. Mitochondrial DNA, 2015, 26(5): 789-790.

47.

Tollefsbol TO. Epigenetics: The new science of genetics[M]. Handbook of Epigenetics. 2011: 1-6.

48.

Valinluck V, Tsai HH, Rogstad DK, et al. Oxidative damage to methyl-CpG sequences inhibits the binding of the methyl-CpG binding domain (MBD) of methyl-CpG binding protein 2 (MeCP2)[J]. Nucleic acids research, 2004, 32(14): 4100-4108.

49.

Volpe P. The language of methylation in genomics of eukaryotes[J]. Biochemistry (moscow), 2005, 70(5): 584-595.

50.

Wang M, Wang P, Liang F, et al. A global survey of alternative splicing in allopolyploid cotton: landscape, complexity and regulation[J]. New Phytologist, 2018, 217(1): 163-178.

51.

Wang Y, Zeng F, Hon C C, et al. The mitochondrial genome of the Basidiomycete fungus Pleurotus ostreatus (oyster mushroom)[J]. FEMS microbiology letters, 2008, 280(1): 34-41.

52.

Wang Y, Wang Z, Liu C, et al. Genome-wide analysis of DNA methylation in the sexual stage of the insect pathogenic fungus Cordyceps militaris[J]. Fungal biology, 2015, 119(12): 1246-1254.

53.

Wu TP, Wang T, Seetin MG, et al. DNA methylation on N6-adenine in mammalian embryonic stem cells. Nature, 2016, 532: 329-333.

54.

Xiong C, Xia Y, Zheng P, et al. Developmental stage-specific gene expression profiling for a medicinal fungus Cordyceps militaris[J]. Mycology, 2010, 1(1): 25-66.

55.

Xu J, Huang B, Qin C, et al. Sequence and phylogenetic analysis of Beauveria bassiana with mitochondrial genome[J]. Mycosystema, 2009, 28(5): 718-23.

56.

Yu L, Wang X, Ting N, et al. Mitogenomic analysis of Chinese snub-nosed monkeys: Evidence of positive selection in NADH dehydrogenase genes in high-altitude adaptation[J]. Mitochondrion, 2011, 11(3): 497-503.

57.

Zemach A, McDaniel IE, Silva P, et al. Genome-wide evolutionary analysis of eukaryotic DNA methylation[J]. Science, 2010, 328(5980): 916-919.

58.

Zhang C, Li C, Ye W, et al. The complete mitochondrial genome of Hericium coralloides (Hericiaceae, Basidiomycota)[J]. Mitochondrial DNA Part B, 2017, 2(2): 385-386.

59.

Zhang Y, Li E, Wang C, et al. Ophiocordyceps sinensis, the flagship fungus of China: terminology, life strategy and ecology[J]. Mycology, 2012, 3(1): 2-10.

60.

Zhang Y, Zhang S, Zhang G, et al. Comparison of mitochondrial genomes provides insights into intron dynamics and evolution in the caterpillar fungus Cordyceps militaris[J]. Fungal Genetics and Biology, 2015, 77: 95-107.

61.

Zheng P, Xia Y, Xiao G, et al. Genome sequence of the insect pathogenic fungus Cordyceps militaris, a valued traditional Chinese medicine[J]. Genome biology, 2011, 12(11): 1.

Beauveria bassiana, Cordyceps militaris and Ophiocordyceps sinensis (Ascomycotina) are traditional Chinese medicines. Here, mitogenomes of these three Ascomycotina fungi were sequenced and de-novo assembled using single-molecule real-time sequencing. The results showed that their complete mitogenomes were 31,258, 31,854 and 157,584 bp, respectively, with sequencing depth approximately 278,760×, 326,283× and 69,385×. Types of repeat sequences were mainly (AA)n, (AAT)n, (TA)n and (TATT)n. DNA methylation motifs were revealed in DNA modifications of these three fungi. We discovered new models of RNA editing through analysis of transcriptomes from B. bassiana and C. militaris. These data lay a solid foundation for further genetic and biological studies about these three fungi, especially for elucidating the mitogenome evolution and exploring the regulatory mechanism of adapting environment. 62.

The highlights of this manuscript include (but not limit): 1. Mitogenomes of B. bassiana, C. militaris and O. sinensis were sequenced and de-novo assembled using single-molecule real-time sequencing. 2. DNA methylation motifs were revealed in DNA modifications of the three fungi mitogenomes. 3. New models of RNA editing were discovered through analysis of transcriptomes from B. bassiana and C. militaris. 63.