genes MDPI Article Transcriptome Changes during Major Developmental Transitions Accompanied with Little Alteration of DNA Methylome in Two Pleurotus Species Jiawei Wen 1.2.t,Zhibin Zhang 1.*Lei Gong 1t,Hongwei Xun2,Juzuo Li1,Bao Qi3 Qi Wang3,Xiaomeng Li 1,Yu Li 3,*and Bao Liu 1 1 Key Laborat catio Normal University, nu.edu.cn (L):lixm441ner ducn (X.L):b nu.edu.cn(B.L) Of Ag the Min un 13 dible andMed Fungi Jilin Agricultural Univer .Changchun 130118,China:qibao3712@163.com(B.Q.); qwang2006@126.com(QW 4@nenu.edu.cn(ZZ):yuli966@126.com(Y.L):Tel:+8-431-8509-9367(Z.Z) 486.431-8453-33090YL5 t These authors contributed equally to this work updates Abstract:Pleurotus tuoliensis(Pt)and P.eryngii var.eryngii(Pe)are important edible mushrooms.The epigene tic and gene expression signature haracterizing major developmental transitions in thes in largely unk e,we rep ort global analyses of DNA methylation an um t t in hot h Pt and Pe th largely stable nes and trans elements (tes)across the developmental transitions.The repressive impact of dna methylation on expression of a small subset of genes is likely due to TE-associated effects rather than their own developmental dynamics.Global expression of gene orthologs was also broadly conserved between Pt and Pe,but discernible interspecific differences exist especially at the fruit body formation stage, s in trans-acting factors. The me hylome and transcriptom ones w or the litate furth rpinning gene expres d evelopmer tus and rolated a Keywords:Pleurotus tuoliensis;P.eryngii;DNA methylation;epigenetics;gene expression; development;interspecific divergence 1.Introduction nd Pleurotus eryngii (Pe)are two famous species of the Pleurotus complex tha ar encomp s the e oys genu e in East 5.D dra al chan driv h ental factors such as ter eriod and culture suhstrates 161 An study identified candidate genes related to mushroom formation in Pt,which was associated with reproductive growth,activation of specific transcription factors,upregulation of genes involved in the Gms2019,10,465doi10.390 Vgenes1006046
genes G C A T T A C G G C A T Article Transcriptome Changes during Major Developmental Transitions Accompanied with Little Alteration of DNA Methylome in Two Pleurotus Species Jiawei Wen 1,2,† , Zhibin Zhang 1,* ,† , Lei Gong 1,† , Hongwei Xun 2 , Juzuo Li 1 , Bao Qi 3 , Qi Wang 3 , Xiaomeng Li 1 , Yu Li 3,* and Bao Liu 1 1 Key Laboratory of Molecular Epigenetics of the Ministry of Education (MOE), Northeast Normal University, Changchun 130024, China; wjw913@sina.com (J.W.); gongl100@nenu.edu.cn (L.G.); lijz346@nenu.edu.cn (J.L.); lixm441@nenu.edu.cn (X.L.); baoliu@nenu.edu.cn (B.L.) 2 Jilin Academy of Agricultural Sciences, Changchun 130033, China; xunhw334@nenu.edu.cn 3 Engineering Research Center of the Ministry of Education (MOE) for Edible and Medicinal Fungi, Jilin Agricultural University, Changchun 130118, China; qibao3712@163.com (B.Q.); q_wang2006@126.com (Q.W.) * Correspondence: zhangzb554@nenu.edu.cn (Z.Z.); yuli966@126.com (Y.L.); Tel.: +86-431-8509-9367 (Z.Z.); +86-431-8453-3309 (Y.L.) † These authors contributed equally to this work. Received: 17 April 2019; Accepted: 12 June 2019; Published: 17 June 2019 Abstract: Pleurotus tuoliensis (Pt) and P. eryngii var. eryngii (Pe) are important edible mushrooms. The epigenetic and gene expression signatures characterizing major developmental transitions in these two mushrooms remain largely unknown. Here, we report global analyses of DNA methylation and gene expression in both mushrooms across three major developmental transitions, from mycelium to primordium and to fruit body, by whole-genome bisulfite sequencing (WGBS) and RNA-seq-based transcriptome profiling. Our results revealed that in both Pt and Pe the landscapes of methylome are largely stable irrespective of genomic features, e.g., in both protein-coding genes and transposable elements (TEs), across the developmental transitions. The repressive impact of DNA methylation on expression of a small subset of genes is likely due to TE-associated effects rather than their own developmental dynamics. Global expression of gene orthologs was also broadly conserved between Pt and Pe, but discernible interspecific differences exist especially at the fruit body formation stage, and which are primarily due to differences in trans-acting factors. The methylome and transcriptome repertories we established for the two mushroom species may facilitate further studies of the epigenetic and transcriptional regulatory mechanisms underpinning gene expression during development in Pleurotus and related genera. Keywords: Pleurotus tuoliensis; P. eryngii; DNA methylation; epigenetics; gene expression; development; interspecific divergence 1. Introduction Pleurotus tuoliensis (Pt) and Pleurotus eryngii (Pe) are two famous species of the Pleurotus eryngii complex that encompass the largest number of species in the oyster mushroom genus [1–3]. Both species are commercially important and widely cultivated especially in East Asia [4,5]. During sexual reproduction, basidiomycetes species undergo dramatic morphological changes driven by environmental factors such as temperature, photoperiod and culture substrates [6]. A previous study identified candidate genes related to mushroom formation in Pt, which was associated with reproductive growth, activation of specific transcription factors, upregulation of genes involved in the Genes 2019, 10, 465; doi:10.3390/genes10060465 www.mdpi.com/journal/genes
Gms2019,10,465 2of14 carbohydrate metabolism pathway and cold and light responses 17l.However.whole trans analysis of pe during developmental stages has not been reported.moreover.with circa 18 million years divergence,significant morphological variations(especially in fruit body)have evolved between the two species[].Thus,an important question to ask is what are the genetic bases and molecular mechanisms underpinning the differences in growth habit and phenotypic transformation during important developmental mycelim to primordim and to fruit body in the vith other for modification (e.ghisto arkers),DNA methylation is 9.101.DNA methyla nst ani d fu and is involved in the regulation of diverse biological pr mic in development,transposable elements(TEs)silencing and overall control of g ession [-14 In mammals,DNA methylation is confined to the symmetric CG context,whereas in plants and fungi,it may occur at all cytosine bases classified into three sequence contexts,CG,CHG and CHH,where H =A,T or C[15-17].In animals and plants,promoter and coding region of protein-coding genes often show moderate els of DNA methylation(prir ly in CG conte wh h associates with gen on neither an e or positive manner [1 st,it wa dhatthegernera ly d nd othe that i plants and e DNA methylation mainly retains the ancient function of te silencing (the s ome defense theory)in the fungal kingdom(Eumycetes).Recently,TE methylation has been further studied in Eumycota.For example,Montanini et al.found that incomplete TE methylation plays important roles in promoting genome plasticity [17].Notably,TE-associated methylation contributing to silencing of adjacent genes was also demo Ample stt s in nu that DNAn n is dynamic acros albeit to ad to ole.fruit rip g and rm development are accompanied with DNA hy omethyaltion 14.2 et al found that DNA methvlation in both genic-regions and TEs also showed moderate dynamic changes during appressoria formation in the pathogenic fungus Magnaporthe Oryza 26].Likewise,similar extent of methylation dynamics was observed during the sexual development of Cordyceps militaris [27]. n contrast,1 was reported that gene expression change rather than methylation reprogramming triggers fruit b e tr unravel the geneti and molecula at ole it tha ssed div atly hetw bagasse-degraded fungi,Aspergillus niger and Trichoderma reesei28).Orthologsinvolved in pathogenesis, including oxalate bios vnthesis and endo-polygalacturonases.were differentially expressed betweer two plant pathogens,Sclerotinia sclerotiorum and S.trifoliorum,contributing to their different host ranges.aoug the tw ngand ndicared that ticot orthols verg in thi body stage DNA ,11 (WGBS) ortant Ple methylome at three antal stas ordium and fruit bo y and we cond comparative analvses.our that maior developmental transitions in the two mushroom species are associated with extensive changes of transcriptome but little alteration of DNA methylome
Genes 2019, 10, 465 2 of 14 carbohydrate metabolism pathway and cold and light responses [7]. However, whole transcriptome analysis of Pe during developmental stages has not been reported. Moreover, with circa 18 million years divergence, significant morphological variations (especially in fruit body) have evolved between the two species [8]. Thus, an important question to ask is what are the genetic bases and molecular mechanisms underpinning the differences in growth habit and phenotypic transformation during important developmental transitions, i.e., from mycelium to primordium and to fruit body in the two mushrooms? Compared with other forms of epigenetic modification (e.g., histone markers), DNA methylation is a relatively stable and heritable marker [9,10], DNA methylation exists in most animals, plants and fungi, and is involved in the regulation of diverse biological processes such as genomic imprinting, organ development, transposable elements (TEs) silencing and overall control of gene expression [11–14]. In mammals, DNA methylation is confined to the symmetric CG context, whereas in plants and fungi, it may occur at all cytosine bases classified into three sequence contexts, CG, CHG and CHH, where H = A, T or C [15–17]. In animals and plants, promoter and coding region of protein-coding genes often show moderate levels of DNA methylation (primarily in CG context), which associates with gene expression in either a negative or positive manner [11,15]. In contrast, it was found that the general landscape of CG methylation in representative fungi (belonging to ascomycete, basidiomycetes and zygomycetes) is largely depleted in genic regions and mainly enriched in TEs and other repetitive sequences [18]. This suggests that, in contrast to the situation in plants and especially animals, DNA methylation mainly retains the ancient function of TE silencing (the genome defense theory) in the fungal kingdom (Eumycetes). Recently, TE methylation has been further studied in Eumycota. For example, Montanini et al. found that incomplete TE methylation plays important roles in promoting genome plasticity [17]. Notably, TE-associated methylation contributing to silencing of adjacent genes was also demonstrated in several Pleurotus species [8,19,20]. Ample studies in human and animals have established that DNA methylation is dynamic across development. Genome-wide DNA methylation reprogramming occurred in mouse primordial germ cells and pre-implantation embryo [21,22]. In plants, methylome dynamics are also known to occur albeit to a lesser extent and/or confined to very specific cell types. For example, fruit ripening and endosperm development are accompanied with DNA hypomethyaltion [14,23–25]. Jeon et al. found that DNA methylation in both genic-regions and TEs also showed moderate dynamic changes during appressoria formation in the pathogenic fungus Magnaporthe Oryza [26]. Likewise, similar extent of methylation dynamics was observed during the sexual development of Cordyceps militaris [27]. In contrast, it was reported that gene expression change rather than methylation reprogramming triggers fruit body development in Pleurotus ostreatus [20]. Comparative transcriptome analysis is a powerful tool to unravel the genetic and molecular bases underpinning divergence of growth and development as well as differential environmental adaption between related organismal species. For example, it was found that enzymes involved in biomass degradation and related transcription factors were expressed divergently between two bagasse-degraded fungi, Aspergillus niger and Trichoderma reesei[28]. Orthologs involved in pathogenesis, including oxalate biosynthesis and endo-polygalacturonases, were differentially expressed between two plant pathogens, Sclerotinia sclerotiorum and S. trifoliorum, contributing to their different host ranges, although the two fungi are morphologically similar [29]. In basidiomycetes, comparisons among Coprinopsis cinerea, Laccaria bicolor and Schizophyllum commune indicated that most orthologs showed divergent expression at the fruit body stage [30]. In this study, we generated single-base resolution DNA methylomes by whole-genome bisulfite sequencing (WGBS) for the two commercially important Pleurotus species, P. tuoliensis and P. eryngii, at three major developmental stages, mycelium, primordium and fruit body, and we conducted comparative analyses. Our results suggest that major developmental transitions in the two mushroom species are associated with extensive changes of transcriptome but little alteration of DNA methylome
Gs2019.10.465 3of14 2.Materials and Methods 2.1.Strains and Culture Condition All samples used for whole- ome bisulfite sequencing (WGBS)and RNA-se g were collected from identical batches of cultivated Pleurotus eryngii var.eryngii (Pe,strain ID:JKXB130DA)and Pleurotus tuoliensis (Pt,strain ID:JKBL130LB)at the three major developmental stages:mycelium, primordium and fruit body.Both monokaryotic mycelia were cultivated in PDA (Potato De extrose Agar)liquid medium in dark at 23 days witl naking culture(120 rpm).Mycelium fi by sterile gauze was read extraction.In or um,first v.For Pt cold siml k d ay. formed at 0C for 48 h.Th d P vere under hle light condition (3001x-10001y)at 14C for imordium initiation aftor 10 davs of cultivation.primordia with 2-3 cm were sampled from the bottles.finally.the remaining cultures were transferred to a light(15C)/dark(7C)photoperiod of 12h condition for about 15 days to induce fruit body formation.For RNA-seq,three biological replicates per condition were separately sampled.All samples were stored in liquid nitrogen until DNA and RNA extractions.Morphology of all sequencing samples were shown in Figure Sl 2.2.Whole-Genome Bisulfite Sequencing and Data Processing Genomic DNA (NA)from primordium and fruit body of P sis and P was extracte B.CA eared by son to th ng first Trimn nomatic [32]wa uality ads The the clean data were aligned to the corresr ponding draft genomes of Pt and Pe with 1-bp mismatch by Bismark program [33].Only uniquely mapped reads were used for further analysis. 2.3.Differential Methylation Analysis Cytosine sites in CG-contexts with more than four uniquely mapped reads were used to further analysis.To decide methylated site,bisulfite non-conversion rate of%evaluated from unmethylated ADNA was used as background methylation leve .Then,bin mial test was performed for each cyto ide if the observed methylation level signifi ntly higher than such ackground methylatio . above adjusted e metho ned a ger tha roach with lkb size with step size of 200-bp was used to identify DMRs.Only windows including at least 10 mCG in at least one tissue were for further DMRs identification.For each window,Fisher's exact test was performed to measure the difference of methylation level (reads count as agent)betweer adjacent developmental stages.P-values of above test were adjusted by Benjamini-Ho chberg methoc ues Windows with 0.01 and the changes of methylation level 0.15 were defined 5 rlap ng window with identical DMRs characteristics (hyper-or hypo-methylation) were am ikb of g enes)and g e bodies including at least 10 mcq/kb ir t loast ono (MP)and ne body (MGB) respectively,and were used for further identification of different methylation level described as above 2.4.RNA-seq and Data Process mRNA extracted from mycelium,primordium and fruit body of P.tuoliensis and P ervngii were sequenced by the Illumina Hiseq 2500 platform.Raw data were filtered by removing low-quality reads
Genes 2019, 10, 465 3 of 14 2. Materials and Methods 2.1. Strains and Culture Conditions All samples used for whole-genome bisulfite sequencing (WGBS) and RNA-seq were collected from identical batches of cultivated Pleurotus eryngii var. eryngii (Pe, strain ID: JKXB130DA) and Pleurotus tuoliensis (Pt, strain ID: JKBL130LB) at the three major developmental stages: mycelium, primordium and fruit body. Both monokaryotic mycelia were cultivated in PDA (Potato Dextrose Agar) liquid medium in dark at 23 ◦C for 7 days with shaking culture (120 rpm). Mycelium filtered by sterile gauze was ready for nucleic acid extraction. In order to acquire primordium, first, liquid cultures were transformed to cultivation bottles and grown in dark at 25 ◦C for 60 days with 70% humidity. For Pt, cold simulation was performed at 0 ◦C for 48 h. Then, cultivation bottles of Pt and Pe were under blue light condition (300lx–1000lx) at 14 ◦C for primordium initiation. After 10 days of cultivation, primordia with 2–3 cm were sampled from the bottles. Finally, the remaining cultures were transferred to a light (15 ◦C)/dark (7 ◦C) photoperiod of 12 h condition for about 15 days to induce fruit body formation. For RNA-seq, three biological replicates per condition were separately sampled. All samples were stored in liquid nitrogen until DNA and RNA extractions. Morphology of all sequencing samples were shown in Figure S1. 2.2. Whole-Genome Bisulfite Sequencing and Data Processing Genomic DNA (gDNA) from mycelium, primordium and fruit body of P. tuoliensis and P. eryngii was extracted by CTAB method and sheared by sonication to the size range of 200 to 300 bp. Library construction and Bisulfite treatment were processed as described in Hu et al. [31]. Sequencing was performed on the Illumina Hiseq 2500 platform (Illumina; San Diego, USA) with standard protocols. For data processing, first, Trimmomatic [32] was used to remove low-quality sequencing reads. Then, the clean data were aligned to the corresponding draft genomes of Pt and Pe with 1-bp mismatch by Bismark program [33]. Only uniquely mapped reads were used for further analysis. 2.3. Differential Methylation Analysis Cytosine sites in CG-contexts with more than four uniquely mapped reads were used to further analysis. To decide methylated site, bisulfite non-conversion rate of 0.3% evaluated from unmethylated λDNA was used as background methylation level. Then, binomial test was performed for each cytosine site to decide if the observed methylation level significantly higher than such background methylation level [34,35]. P-values of the above tests were then adjusted by Benjamini-Hochberg method to q-values. Cytosine sites with q-value lower than 0.01 and corresponding methylation level larger than 5% were defined as methylated CG context (mCG). Differentially methylated regions (DMRs) between adjacent developmental stages were identified as described in [14] with minor modification. A sliding-window approach with 1kb size with step size of 200-bp was used to identify DMRs. Only windows including at least 10 mCG in at least one tissue were for further DMRs identification. For each window, Fisher’s exact test was performed to measure the difference of methylation level (reads count as agent) between adjacent developmental stages. P-values of above test were adjusted by Benjamini-Hochberg method to q-values. Windows with a q-value < 0.01 and the changes of methylation level ≥ 0.15 were defined as DMRs. Overlapping window with identical DMRs characteristics (hyper- or hypo-methylation) were merged into a large region. Similarly, promoters (upstream 1kb of genes) and gene bodies including at least 10 mCG/kb in at least one tissue were defined as methylated promoter (MP) and methylated gene body (MGB), respectively, and were used for further identification of different methylation level described as above. 2.4. RNA-seq and Data Process mRNA extracted from mycelium, primordium and fruit body of P. tuoliensis and P. eryngii were sequenced by the Illumina Hiseq 2500 platform. Raw data were filtered by removing low-quality reads
Gms2019,10,465 4of14 Filtered reads were aligned to each draft genome by hisat2 1361 with default parameters.read with uniquely mapped gene or transposable element (TE)(annotation as described in [8))were kept for further analysis.The raw counts matrix was extracted by Stringtie with the provided python script (http://ccb.jhu.edu/software/stringtie/dl/prepDEpy)and DESeq2 [37]was used to identify differentiall expressed genes(DEGs;q-value<0.01 and log-transformed fold change>1).Gene and TE with TPM on reads)large than I were taken as ex pressed.Gene expr e copy ort 0log pa en Pt ar wise by BLAT [391.o ith the d the s30nhm d for further analysis er than 70% All the clean data including both WGBS and RNA-seq data have been deposited at the SRA data base http://www.ncbi.nlm.nih.gov/sra/under accession number of PRJNA548464. 2.5.Gene Ontology (GO)Enrichment Analysis To trace the biological process of differentially regulated orthologs,GO enrichment analysis were performed by using Clusterprofiler [40].Go terms with q-value<0.05 and including at least five genes orthologs)were deemed as significantly enriched. 2.6.Statistics All Statistical tests in this paper were performed using basic packages in R language(Version 3.4.3,https://www.-project.org/) 3.Results 3.1.DNA Metlnylation Landscapes in the Two Mushroom Species We generated single-base resolution DNA methylomes by whole-genome bisulfite sequencing (WGBS)in P.tuoliensis and P eryngii,at three developmental stages,mycelium (MY),primordium PR)and fruit bo y (FB).We focused on ion profiles of CG conte pencingdep es 181 xts with ere comm en ntruct CG-co (mCCs)a those with a-values lower than 0.0l and methvlation levels higher than%(Materials and methods) Based on this criterion,22.3%(1.24 M)and 19.0%(0.99 M)CG contexts were mCGs in Pt and Pe, espectively,across the tages.We found that different developmental stages showed similar mCC levels in both species(16. o,15.8%and15.6oinP16.5o,16. and 16.9%in Pe).Similar to previou eports in rungi /mCG levels displayed a typical bimodal distributions in both mushroom Of a Fe): the remain d in al rg ter exon intron and downstream ikh of ger ns enriched in TEs,the primary peak of mCGs was located in high methylation regions whereas the primary peak of gene-related features was located to low methylation regions(Figure 1C) Pt and Pe showed the same trend in mCGs differences among the genomic features,ie.,both were intergenic region>gene downstream 1kb region>promoter region exon=intron(Games-Howell post-hoc test)
Genes 2019, 10, 465 4 of 14 Filtered reads were aligned to each draft genome by Hisat2 [36] with default parameters. Read with uniquely mapped gene or transposable element (TE) (annotation as described in [8]) were kept for further analysis. The raw counts matrix was extracted by Stringtie with the provided python script (http://ccb.jhu.edu/software/stringtie/dl/prepDE.py) and DESeq2 [37] was used to identify differentially expressed genes (DEGs; q-value < 0.01 and log-transformed fold change >1). Gene and TE with TPM (Transcripts Per Million reads) large than 1 were taken as expressed. Gene expression information of both Pt and Pe were given in Supplementary dataset. for Single copy ortholog pairs between Pt and Pe were identified by OrthoMCL [38] after performing BLASTP with e-value < 1 × 10−5 . After aligning pairwise by BLAT [39], orthologs with the proportion of consensus region larger than 70% and the length of consensus region > 300 bp were reserved for further analysis. All the clean data including both WGBS and RNA-seq data have been deposited at the SRA data base http://www.ncbi.nlm.nih.gov/sra/ under accession number of PRJNA548464. 2.5. Gene Ontology (GO) Enrichment Analysis To trace the biological process of differentially regulated orthologs, GO enrichment analysis were performed by using Clusterprofiler [40]. Go terms with q-value < 0.05 and including at least five genes (orthologs) were deemed as significantly enriched. 2.6. Statistics All Statistical tests in this paper were performed using basic packages in R language (Version 3.4.3, https://www.r-project.org/). 3. Results 3.1. DNA Methylation Landscapes in the Two Mushroom Species We generated single-base resolution DNA methylomes by whole-genome bisulfite sequencing (WGBS) in P. tuoliensis and P. eryngii, at three developmental stages, mycelium (MY), primordium (PR) and fruit body (FB). We focused on DNA methylation profiles of CG context only because of its predominance in the two Pleurotus species [8]. We found from 60% to 78% of CG contexts with sufficient sequencing depths (at least five reads at a given stage) were common between the two mushroom species across all stages. Methylation of these common CG contexts with high reliability were chosen to construct CG-context DNA methylation landscapes. We defined methylated CG contexts (mCGs) as those with q-values lower than 0.01 and methylation levels higher than 5% (Materials and methods). Based on this criterion, 22.3% (1.24 M) and 19.0% (0.99 M) CG contexts were mCGs in Pt and Pe, respectively, across the stages. We found that different developmental stages showed similar mCG levels in both species (16.0%, 15.8% and 15.6% in Pt; 16.5%, 16.2% and 16.9% in Pe). Similar to previous reports in fungi [17], mCG levels displayed a typical bimodal distributions in both mushrooms. Of all mCGs, 80.4% (Pt) and 77.4% (Pe) showed mCG levels lower than 10%, whereas the remaining showed levels > 70% (Figure 1A). More than 50% of mCGs were located in intergenic regions wherein TEs abound (Figure 1B). Typical bi-modal distribution of mCGs was observed in all types of genomic features (intergenic region, promoter, exon, intron and downstream 1kb of gene body) (Figure 1C). For intergenic regions enriched in TEs, the primary peak of mCGs was located in high methylation regions; whereas the primary peak of gene-related features was located to low methylation regions (Figure 1C). Pt and Pe showed the same trend in mCGs differences among the genomic features, i.e., both were intergenic region > gene downstream 1kb region > promoter region > exon = intron (Games-Howell post-hoc test)
Gs2019,10,465 5of14 9 Pe g02 Pt 8.a 0 dow 0.0 025 0.5 07 0.25 0.5 07516 Figure 1.Landscape of DNA methylation in the two mushroom sp ecies,Pleurotus tuoliensis (designated Pt)and P.eryngii var.rymgi(designated Pe).(A)Distribution of DNA methylation levels of all covered prim intergenic,promoter,exon and down-1kb,in Pt and Pe.(C)Density of mCG levels mone the three tissues belonging to different categories of genomic features in Pt and Pe.Dashed lines denote the average methylation levels of different genomic features in each species. We next investigated possible dynamics of mCG levels across the three developmental stages in both species.Differentially methylated regions(DMRs)were identified in both developmental transitions,transition 1(from mycelium to primordium)and transition 2( rom primordium to frui body).We found te DMks than Pt in mush for Pe 005.ln ind 421 than that of h DMRs in ition in hoth iue<005).For both rand hy ssessed the highest proportions although they were also of variable ratios(54-70%for hyper-DMR and 77-83%for hypo-DMR,respectively;Figure S2).These results indicate biased distribution of mCGs in both Pt and Pe(enriched in TE regions and depleted in genic regions).Moreover,changes of mCG levels in the DMRs were apparent with developmental progression (phase 2>phase 1),which again mainly occurred in TE-enriched intergenic regions 3.2.Transcriptional Regulation of Gene Expression in the Two Mushroom Species Toinvestigate profiles of gene expression regulation in the two mushroom species(Pt and Pe). 4360 expres Materials and Me ods)were ide to c wo o or phe from n P t the th ()which is similar with previous studies in plant species 41.However.we note that the correlation of expression levels between the stages were significantly lower than those between biological replicates in each stages of the same species(Table S2),suggesting gene expression changes with development.Moreover,the correlations between gene expression decreased with progression of development
Genes 2019, 10, 465 5 of 14 Genes 2019, 10, x FOR PEER REVIEW 5 of 14 Figure 1. Landscape of DNA methylation in the two mushroom species, Pleurotus tuoliensis (designated Pt) and P. eryngii var. eryngii (designated Pe). (A) Distribution of DNA methylation levels of all covered CG sites (each having at least five reads) in all three tissues, mycelium, primordium and fruit body, in Pt and Pe. (B) Distribution of methylated CG sites (mCGs), in different categories of genomic features, intergenic, promoter, exon and down-1kb, in Pt and Pe. (C) Density of mCG levels among the three tissues belonging to different categories of genomic features in Pt and Pe. Dashed lines denote the average methylation levels of different genomic features in each species. We next investigated possible dynamics of mCG levels across the three developmental stages in both species. Differentially methylated regions (DMRs) were identified in both developmental transitions, transition 1 (from mycelium to primordium) and transition 2 (from primordium to fruit body). We found Pe possessed more DMRs than Pt in both transitions, although both mushrooms showed more DMRs in transition 1 than in transition 2 (Table S1; 180 versus 278 for Pt and 421 versus1428 for Pe; binomial test, p-values < 0.05). Intriguingly, the number of hyper-DMRs was more than that of hypo-DMRs in transition 2 in both mushrooms (binomial test, p-value < 0.05). For both hyper- and hypo-DMRs, intergenic regions possessed the highest proportions although they were also of variable ratios (54%–70% for hyper-DMR and 77%–83% for hypo-DMR, respectively; Figure S2). These results indicate biased distribution of mCGs in both Pt and Pe (enriched in TE regions and depleted in genic regions). Moreover, changes of mCG levels in the DMRs were apparent with developmental progression (phase 2 > phase 1), which again mainly occurred in TE-enriched intergenic regions. 3.2. Transcriptional Regulation of Gene Expression in the Two Mushroom Species To investigate profiles of gene expression regulation in the two mushroom species (Pt and Pe), 4360 expressed single-copy orthologs (Materials and Methods) were identified and used to compare expression changes at the two developmental transitions or phases (phase 1, from mycelium to primordium; phase 2, from primordium to fruitbody) in both species. As shown in Figure 2A, orthologs between the two species exhibited well-correlated patterns of gene expression at the three stages (Pearson r = 0.75–0.84), which is similar with previous studies in plant species [41]. However, we note that the correlation of expression levels between the stages were significantly lower than those between biological replicates in each stages of the same species (Table S2), suggesting gene expression changes with development. Moreover, the correlations between gene expression decreased with progression of development. Figure 1. Landscape of DNA methylation in the two mushroom species, Pleurotus tuoliensis (designated Pt) and P. eryngii var. eryngii (designated Pe). (A) Distribution of DNA methylation levels of all covered CG sites (each having at least five reads) in all three tissues, mycelium, primordium and fruit body, in Pt and Pe. (B) Distribution of methylated CG sites (mCGs), in different categories of genomic features, intergenic, promoter, exon and down-1kb, in Pt and Pe. (C) Density of mCG levels among the three tissues belonging to different categories of genomic features in Pt and Pe. Dashed lines denote the average methylation levels of different genomic features in each species. We next investigated possible dynamics of mCG levels across the three developmental stages in both species. Differentially methylated regions (DMRs) were identified in both developmental transitions, transition 1 (from mycelium to primordium) and transition 2 (from primordium to fruit body). We found Pe possessed more DMRs than Pt in both transitions, although both mushrooms showed more DMRs in transition 1 than in transition 2 (Table S1; 180 versus 278 for Pt and 421 versus1428 for Pe; binomial test, p-values < 0.05). Intriguingly, the number of hyper-DMRs was more than that of hypo-DMRs in transition 2 in both mushrooms (binomial test, p-value < 0.05). For both hyper- and hypo-DMRs, intergenic regions possessed the highest proportions although they were also of variable ratios (54–70% for hyper-DMR and 77–83% for hypo-DMR, respectively; Figure S2). These results indicate biased distribution of mCGs in both Pt and Pe (enriched in TE regions and depleted in genic regions). Moreover, changes of mCG levels in the DMRs were apparent with developmental progression (phase 2 > phase 1), which again mainly occurred in TE-enriched intergenic regions. 3.2. Transcriptional Regulation of Gene Expression in the Two Mushroom Species To investigate profiles of gene expression regulation in the two mushroom species (Pt and Pe), 4360 expressed single-copy orthologs (Materials and Methods) were identified and used to compare expression changes at the two developmental transitions or phases (phase 1, from mycelium to primordium; phase 2, from primordium to fruitbody) in both species. As shown in Figure 2A, orthologs between the two species exhibited well-correlated patterns of gene expression at the three stages (Pearson r = 0.75−0.84), which is similar with previous studies in plant species [41]. However, we note that the correlation of expression levels between the stages were significantly lower than those between biological replicates in each stages of the same species (Table S2), suggesting gene expression changes with development. Moreover, the correlations between gene expression decreased with progression of development