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)