Sexual stage-induced long noncoding RNAs in the filamentous fungus Fusarium graminearum

Long noncoding RNA (lncRNA) plays important roles in morphological differentiation and development in eukaryotes. In filamentous fungi, however, little is known about lncRNAs and their roles in sexual development. Here we describe sexual stage-induced lncRNAs during the formation of perithecia, the sexual fruiting bodies of Fusarium graminearum. We have identified 547 lncRNAs whose expression was developmental stage-specific, with about 40% of which peaked during the development of asci, the sac-like structures containing meiospores. A large fraction of the lncRNAs were found to be antisense to mRNAs, forming 300 sense–antisense pairs. Although small RNAs (sRNAs) were produced from these overlapping loci, most of the antisense lncRNAs appeared not to be involved in gene silencing pathways. Genome-wide analysis of sRNA clusters identified many silenced loci at the meiotic stage. However, we found transcriptionally-active sRNA clusters, many of which were associated with lncRNAs. Also, we observed that many antisense lncRNAs and their respective sense transcripts were induced in parallel as the perithecia matured. To identify regulatory components for lncRNA expression, we analyzed mutants defective in the nonsense-mediated decay (NMD) pathway. A subset of the lncRNAs appeared to be targeted by the NMD before the perithecia formation, suggesting a suppressive role of the NMD in lncRNA expression during vegetative stage. This research provides fundamental genomic resources that will spur further investigations on developmental lncRNAs that may play important roles in shaping the fungal fruiting bodies.


INTRODUCTION
Genomes of eukaryotes-from simple yeast to animals-are pervasively transcribed from noncoding intergenic regions and in antisense orientation from genic regions (Guttman et al. 2009;). Long noncoding RNAs (lncRNAs) are loosely defined as noncoding transcripts longer than 200 nucleotides, which are mostly transcribed by RNA polymerase II and share common features with mRNAs other than protein-coding capacity (Kapranov et al. 2007). lncRNAs are versatile molecules that not only regulate gene expression, but also affect enzymatic activities and chromosome conformation (reviewed in Rinn and Chang 2012;Quinn and Chang 2015). Since the discovery of the XIST lncRNA required for X chromosome inactivation (Brockdorff et al. 1992), the roles of lncRNAs in developmental processes such as embryogenesis and tissue differentiation have been extensively studied in animals along with the advent of RNA-seq technologies (Fatica and Bozzoni 2014;Flynn and Chang 2014). Yet the full scope of the developmental roles of lncRNAs is far from understood.
In the highly divergent yeasts Saccharomyces cerevisiae (budding yeast) and Schizosaccharomyces pombe (fission yeast), the onset of sexual sporulation and the following meiotic divisions are tightly regulated by elaborate mechanisms involving lncRNAs (reviewed in Hiriart and Verdel 2013). In budding yeast, a promoter-derived lncRNA suppresses the expression of IME1 (inducer of meiosis 1), the master regulator for the sexual sporulation, by inducing heterochromatin formation in the promoter region of IME1 during vegetative growth (van Werven et al. 2012). In addition, the transcription of another lncRNA antisense to IME4 gene inhibits the expression of IME4 by antagonizing the sense transcription (Hongay et al. 2006;Gelfand et al. 2011). Although there is no such conserved or Kim  analogous regulatory mechanism in fission yeast, lncRNAs also play diverse roles in the sexual sporulation, for example, sequestering RNA elimination factors that repress meiotic gene expression (Harigaya et al. 2006;Hiriart et al. 2012;Yamashita et al. 2012), and contributing homologous chromosome pairing . Despite the growing evidence of the regulatory roles in yeasts, information on lncRNA expression and function during fruiting body formation in filamentous fungi is scarce.
RNA quality control mechanisms are crucial for the regulation of lncRNA expression in budding yeast. The nuclear exosome is engaged in RNA processing and degradation of transcripts including lncRNAs that are specifically expressed during sexual sporulation; the deletion of RRP6 encoding the exosome-associated exonuclease resulted in the accumulation of noncoding transcripts that otherwise remained silenced during vegetative growth (Davis and Ares 2006;Camblong et al. 2007;Lardenois et al. 2011). In human cells, promoter-derived transcripts were also ectopically expressed upon deletion of the exosome components including the homologous RRP6, suggesting the conserved role of the exosome for lncRNA expression in diverse eukaryotes (Preker et al. 2008). The nonsense-mediated decay (NMD) pathway is another quality control checkpoint for aberrant transcripts in the cytoplasm, and recently emerged as a key player for fine-tuning of both coding and noncoding gene expression (Smith and Baker 2015). A genome-wide survey of human lncRNA sequences showed that most lncRNAs harbor short open reading frames (ORFs) that would lead to activation of the NMD pathway (Niazi and Valadkhan 2012). In fact, subsets of lncRNAs found in budding yeast, the model plant Arabidopsis and animals are subject to degradation through the NMD pathway (Kurihara et al. 2009;Tani et al. 2013; Ruiz-Orera et al. The cytoplasmic exonuclease Xrn1 is the final enzyme responsible for the degradation of decapped and de-adenylated transcripts that have been recognized and processed by NMD components. The deletion of XRN1 in budding yeast also leads to the accumulation of more than a thousand cryptic noncoding transcripts termed 'XUTs' (Xrn1-sensitive unstable transcripts), most of which are distinct from the noncoding transcripts that arise by exosome depletion (van Dijk et al. 2011). Many XUTs are antisense to annotated genes and seemed to have repressive roles in sense transcription by modulating chromatin status of the promoter regions (van Dijk et al. 2011).
It has been argued that organismal complexity is correlated with expression dynamics of noncoding transcripts (Mattick et al. 2010;Liu et al. 2013;Gaiti et al. 2015;Quinn and Chang 2015). In the Kingdom Fungi, multicellular fruiting bodies have independently evolved at least twice in the diverging lineages (Knoll 2011;Nguyen et al. 2017). Given the key regulatory roles of lncRNAs in sexual sporulation and morphological transition in yeasts (Bumgarner et al. 2009;Chacko et al. 2015), lncRNAs may have exerted their roles in evolution of multicellularity and sexual development in filamentous fungi.
Fusarium graminearum is a plant pathogenic fungus infecting our staple crops, such as wheat and corn, and thus has been a model for studying the developmental process of perithecia, the sexual fruiting bodies of the fungus, as well as other interesting aspects of biology including host-pathogen interaction and mycotoxin production (Trail 2009;Ma et al. 2013). The fungus is probably the best organism for investigating the lncRNA catalog in fruiting body-forming fungi, as perithecia develop at sufficient synchronicity in culture media, enabling time-series transcriptome analyses with this microscopic organism (Hallen et al. 2007;Sikhakolli et al. 2012;Trail et al. 2017). Also, the genome sequence assembly is complete, featuring a total 4 chromosomes (Cuomo et al. 2007), and the genome has been annotated and curated-although it still lacks lncRNA annotations (King et al. 2015(King et al. , 2017. In addition, plentiful genetic resources have accumulated through large-scale functional studies of perithecia development (Son et al. 2011;Wang et al. 2011a;Yun et al. 2015;Kim et al. 2015b;Liu et al. 2016;Son et al. 2017;Trail et al. 2017).
The goal of the present study was to characterize lncRNAs that are specifically expressed in the fungal fruiting body undergoing sexual development and to investigate their developmental stage-specific expression and their regulation. We identified lncRNAs with a pipeline that constructs de novo transcript annotations by combining RNA-seq data from vegetative and sexual developmental stages, and then removes those with detectable protein-coding potential or monotonous expression profile. Hundreds of lncRNAs that exhibit dynamic expression patterns were found, thereby expanding the universe of genomes known to have significant noncoding roles in development-specifically, here, the multicellular development of fungal fruiting bodies.

Transcriptional reprogramming during perithecia formation
To obtain time-course transcriptome data during the sexual development of F. graminearum, we sequenced samples from hyphae, strands of cells that make up vegetative stages of most of fungi (S0), and from five successive sexual stages (S1-S5; Fig. 1A) that capture key morphological transitions during the perithecia development, defined at the formation of: S1-formation of perithecial initials (hyphae curls that give rise to the perithecial tissues), S2-perithecial walls, S3-paraphyses (sterile cells Kim et al, 2018 Page 6 supporting perithecia), S4-asci (sac-like structures in which ascospores develop), and S5-ascospore maturation (Trail and Common 2000). A total of 480 million RNA-seq reads were generated from 18 samples (6 stages × 3 replicates), and there were average 25 million mapped reads per sample (Supplemental Fig. S1). We validated our sampling scheme by perithecial morphology for 3 biological replicates, using the BLIND program (Anavy et al. 2014) that determined the sequence of the developmental time-course data without prior information other than gene expression data (Supplemental Differentially-expressed (DE) genes between any two successive developmental stages (>4-fold at 5% false discovery rate [FDR]) were mostly unique when compared to other pairwise stage comparisons (Fig. 1B). Overrepresented GO terms for the stage-specific DE genes reflected key biological processes during the morphological transitions (Supplemental Table S1). For example, the GO term 'lipid metabolism' had the highest representation in the 'S0 vs. S1' comparison, although statistically non-significant (Fig. 1C). The accumulated lipids in hyphae and perithecial initials are vital for paraphyses and asci development (Guenther et al. 2009). Perithecia dramatically increased in size and became more rigid during S2 and S3, which is accompanied by GO terms related to carbohydrate metabolisms (Fig. 1C). Finally, when the asci develop from the fertile layer during S3 and S4, meiosisrelated genes were significantly enriched for GO terms including 'meiotic cell cycle process' and 'ascospore formation' at 5% FDR (Fig. 1C). It is noteworthy that most of the DE genes were upregulated in the later developmental stages (Supplemental Fig. S3), indicating that gene activation is the most common means of gene regulation during the sexual development. showing the number of differentially expressed (DE) genes between two successive developmental stages (>4-fold; 5% FDR). Note that most of the DE genes were unique in each comparison. (C) Functional enrichment analyses for DE genes between two successive developmental stages. Fifty-three GO termsthat can be broadly categorized into 7 biological processes-were assessed for degree of functional enrichment (Supplemental Table S1), and were projected to two-dimensional semantic spaces. Only GO terms with P < 0.05 were depicted in each panel.

Identification of lncRNA in perithecia
To discover lncRNAs expressed during the perithecia development, we adopted an established protocol for novel transcripts identification, with some modifications (Weirick et al. 2016;Supplemental Fig. S4). First, we constructed de novo transcript annotations (28,872 transcripts expressed from 20,459 genomic loci), and identified potentially novel transcripts that were absent in the reference annotations (Table 1). For noncoding transcripts identification, coding potential of the novel transcripts was computed by using the CPAT program . To maximize both sensitivity and specificity for noncoding transcript detection, the program was trained on the F. graminearum genome dataset and the threshold was set to a CPAT score of 0.540 (Supplemental Fig. S4) (cf. 0.364 for humans, 0.440 for mice; Chakraborty et al. 2014). The transcripts with a low coding potential (CPAT score ≤ 0.540) were further scanned against Pfam and Rfam databases to filter out transcripts encoding protein domain(s) and harboring any known structural RNA motifs, respectively (E < 10 −10 ; Table 1). Finally we only retained transcripts that were differentially expressed in at least one developmental stage (5% FDR), yielding a total 547 lncRNA candidates (Table 1; Supplemental Table S2).
The identified lncRNAs were distributed across the 4 chromosomes, and generally shorter with fewer exons, when compared to mRNAs (transcripts with CPAT score > 0.540) ( Fig. 2A-C). Based on the relative position to mRNAs, lncRNAs can be classified as antisense lncRNA (ancRNA) or long intergenic ncRNA (lincRNA). There were 280 ancRNAs that overlapped more than 100 bp of an mRNA on the opposite strand, and 237 lincRNAs that were situated between annotated genes (Supplemental Table S2). The mean AU content of ancRNA and lincRNA sequences falls between the coding sequences of the mRNAs and the intergenic regions ( sequences are commonly observed in other eukaryotes (Nam and Bartel 2012;Pauli et al. 2012;Gaiti et al. 2015;Nyberg and Machado 2016). However, we identified only 5 lncRNAs that showed similarity with lncRNAs in other eukaryotes (E < 10 −10 ; Supplemental Table S3), suggesting either the poor status of lncRNA annotations in filamentous fungi or a high degree of sequence divergence in fungal lncRNAs. b Among total 11 transcript class codes, transcripts tagged with the class codes 'J', 'P', 'U', and 'X' were considered as novel transcripts (see Supplemental Methods).
c Noncoding transcripts were identified by the coding potential assessment tool program , and further filtered by Pfam and Rfam database searches.

Developmental expression of lncRNA
The sexual stage transcriptome data showed predominance of lncRNAs at the meiotic stage (S4) where the expression of many lncRNAs peaked (234 out of the 547 lncRNAs; Fig. 2E). We compared the degree of differential expression of lncRNAs to that of mRNAs (9,457 transcripts differentially expressed in at least one developmental stage at 5% FDR). The ratio of the maximum expression among six developmental stages to the mean expression over the remaining five stages was calculated for lncRNAs and the differentially expressed mRNAs. By this metric, lncRNA was prone to be more differentially expressed than mRNA (P = 2.2 × 10 −16 , KS-test statistic D = 0.39), with the third quartile value of the ratio measuring 2.93 for lncRNA and 1.88 for mRNA ( Fig. 2F). Also, we identified seven co-expressed clusters of lncRNAs that showed developmental-specific expression patterns, suggesting distinct roles of lncRNAs in different stages of perithecia development (Fig. 3A).
In addition to the sexual stage dataset, we obtained transcriptome data during spore germination to investigate the degree of lncRNA expression in vegetative stages. The dataset was comprised of four spore germination stages: G0-fresh spore, G1-polar growth, G2-doubling of long axis, and G3- polyadenylation sites were determined (Supplemental Fig. S6). Also, to examine if there is an intraspecific conservation in the lncRNA content and expression, we utilized degradome-seq data of another F. graminearum wild-type strain sampled at meiotic stage (Son et al. 2017). Degradation of lncRNA transcripts expressed at S4 was evident, which in turn confirmed the consistent lncRNA production in the two strains ( Fig. 4; also see below, Fig. 5B).
Fungal genomes are known for having shorter intergenic spaces, compared to other eukaryotic genomes. Therefore, fungal lncRNAs could be a transcriptional noise arising from neighboring genes. To test this, we examined global patterns of the expression correlation between lncRNA and neighboring genes, with close examination of some selected examples whose expression was confirmed by 3' RACE-PCR. The lncRNAs antisense to ORC1, ORC2, and CENP-T each had an upstream gene in close proximity on the same strand ( Fig. 4; Supplemental Fig. S6). However, the sexual stage expression between the lncRNAs and their respective upstream genes was not correlated (| r | < 0.50; Supplemental Table S4). Positive correlation was observed in expression levels between lncRNAs and their divergently transcribed genes as in the HIR1 and NSE4 loci (r > 0.8; Fig. 4; Supplemental Table S4), indicating prevalence of bidirectional promoters for lncRNA transcription (Neil et al. 2009;Xu et al. 2009;Pelechano et al. 2013). On the other hand, the expression of lncRNAs and their convergently transcribed genes tend not to be correlated, as in the CENP-T and RMD1 loci (| r | < 0.50; Fig. 4; Supplemental Table   S4). These patterns were globally observed in lncRNA-associated loci (Supplemental

Identification of sRNA-enriched loci associated with lncRNAs
We found 300 sense mRNA-ancRNA pairs with different orientations: 5'-to-5' partial (70), 3'to-3' partial (61), and full overlaps (169). One of the most common mechanisms involving antisense transcription is the RNAi pathway incorporating sRNAs generated from the double-stranded RNA regions. To investigate the degree and effect of sRNA production in the ancRNA loci, we analyzed the previously published sRNA-seq and degradome-seq data at meiotic stage (Son et al. 2017). As expected, sRNA reads were mapped at a higher frequency to ancRNAs than to mRNAs without overlapping antisense transcripts, the sense mRNAs or lincRNAs ( Fig. 5A; P < 1.2 × 10 −8 , Kruskal-Wallis test statistic H = 155), suggesting that the ancRNA loci may serve as a major source for endogenous sRNA production. However, the correlation of the degradome-seq and our RNA-seq data at S4 showed that sRNA-mediated endonucleolytic cleavage of ancRNAs and sense mRNAs was comparable to that of mRNAs without antisense transcripts (Fig. 5B), implying that the ancRNA loci were not preferentially targeted by RNAi machinery, post-transcriptionally.
It remains paradoxical that gene silencing induced by heterochromatin formation requires sRNA production via co-transcriptional processes, sometimes from lncRNAs (Motamedi et al. 2004;Bühler et al. 2006;Zhang et al. 2008;Bayne et al. 2010;Dang et al. 2016). To search for any lncRNAs associated with sRNA-enriched loci that could be indicative of transcriptional gene silencing events, we examined top 80 sRNA clusters ranked by the number of mapped reads, which accounted for 62% of mapped sRNA reads (Supplemental Table S5). Production of sRNAs in the top 80 clusters was dependent on Dicers and Argonauts, indicating that the sRNAs were produced by RNAi machinery (Fig. 5C) sRNA reads were mapped (Fig. 5D). We observed that the coding genes closest to the centers of sRNA clusters exhibited overall low expression (RPKM < 0.5 in 22 out of the 80 clusters) (Fig. 5D). A significant portion of sRNAs were also derived from noncoding transcripts and genic regions overlapped with antisense transcripts, some of which were identified as lncRNAs ( Fig. 5D; Supplemental Fig. S7).

Co-expression of lncRNAs and their sense transcripts
We could not find strong evidence for sRNA-mediated transcriptional and posttranscriptional gene silencing in the ancRNA loci. Interestingly, we did observe gene expression correlation in many ancRNAs and sense mRNA pairs across the sexual stages ( Fig. 6; 85 out of the 300 pairs with Pearson's correlation | r | > 0.70; Fisher's exact test, P < 0.05), most of which were positively correlated (76 out of the 85 pairs). We asked whether the ancRNAs are antisense to genes involved in a specific biological process. Notably, the positively correlated sense mRNAs were enriched for the GO term 'DNA metabolism' at 5% FDR (Supplemental Table S6). This observation might be a consequence of the underlying structure of the dataset, in which the GO term was also the most overrepresented GO term for all the sense mRNAs overlapping the 280 lncRNAs (although this observation was not statistically significant, applying a 5% FDR; P = 0.006).
Kim Similarly, defects in ascospore formation were observed in heterokaryotic ∆dbp2 strain containing wildtype nuclei (Fig. 7A), suggesting that maintaining RNA homeostasis by NMD pathway is crucial for meiotic divisions and ascospore formation.
To identify lncRNA whose expression was controlled by NMD activity, we analyzed the transcriptome data from the wild-type (WT) and the ∆xrn1 strain at hyphae stage (S0) and meiotic stage (S4). The transcriptome data of the ∆xrn1 strain were distinct, but most similar, to that of WT at the corresponding stages, indicating drastic effects of XRN1 on gene expression levels as a major component for RNA turnover (Fig. 7B). After expression values were normalized with 124 ribosomal protein genes that were known to be relatively insensitive to NMD activity ( Table S7). After the normalization, the XUTs identified at S0 and S4 showed a median 5.0-and 5.8-fold increase in ∆xrn1, respectively ( Fig. 7D; Supplemental Table S7). Many XUTs were previously annotated transcripts or isoforms of them (70%; 781/1,122), and only 11% of XUTs (122/1,122) were predicted to be noncoding transcripts (CPAT score ≤ 0.540) (Supplemental Table S7).
Among the noncoding XUTs, we identified 25 lncRNAs whose expression was elevated upon the Xrn1 depletion at S0 and showed increasing patterns across the sexual stages ( Fig. 7E; Supplemental Fig. S9).

DISCUSSION
Here we profiled transcriptomes of vegetative and sexual stages that span the entire life cycle of F. graminearum. lncRNAs are usually an order of magnitude less abundant than mRNAs (Mercer et al. 2011;Cabili et al. 2015), so the unprecedented depth of the sequencing data we generated (total 938 million mapped reads) enabled us to capture scantly expressed noncoding transcripts. This study has revealed global properties of lncRNAs during the perithecia development, characterized by dynamic and developmental stage-specific expression. In the last step of our lncRNA identification pipeline, we only included differentially expressed noncoding transcripts to discern lncRNAs of biological significance.
This filtering step allowed us to identify low-abundance lncRNAs that alone could be argued to be transcriptional noise (ENCODE Project Consortium 2007;Bakel et al. 2010). For many antisense lncRNAs, we detected parallel induction with their respective sense transcripts across the sexual development and identified a subset of lncRNAs that are sensitive to nonsense-mediated decay activity before sexual induction.
Our lists of F. graminearum lncRNAs contain many confident lncRNA annotations, but are still far from complete. It is primarily due to difficulty in unequivocal determination of whether a transcript is coding or noncoding (Anderson et al. 2015). Although protein products from most of the bonafide lncRNAs with predicted ORFs have not been detected in cells (Bánfai et al. 2012;Gascoigne et al. 2012 These novel transcripts were initially filtered out by the CPAT program, however another annotation tool, CPC2 (Kang et al. 2017) classified them as lncRNAs, and the lncRNA annotations were also supported by the lack of cross-species conservation of the deduced polypeptides (E ≥ 10 −10 ; Supplemental Table S8).
Another source of false-negatives in our dataset may have arisen from poly-A-based library preparation that excluded some lncRNAs lacking poly-A tail. This may account for the inclusion of fewer noncoding transcripts in our XUT identification process. More precise identification of lncRNAs, especially for those undergoing deadenylation by the RNA quality control machineries, should be accompanied by cap-analysis gene expression (CAGE)-seq (Bogu et al. 2016;Wery et al. 2016;Liu et al. 2017). However, most lncRNAs identified in this study were also detected in the degradome-seq ( We were unable to test the role of the exosome components in lncRNA expression, as the major exosome component (RRP6) and the functionally-related helicase genes (DBP2 and MTR4) all seemed to be essential in F. graminearum. In fission yeast, the nuclear exosome mediates gene silencing of coding and noncoding loci during vegetative growth, and the sexual sporulation is triggered by disassembly of heterochromatin on the silenced loci (Zofall et al. 2012). Therefore, it will be interesting to see if the regulation of lncRNA expression is achieved by modulation of chromatin status before and after the sexual induction in F. graminearum.
Antisense transcription can influence synthesis, expression kinetics, and stability of sense transcripts through a variety of mechanisms (Böhmdorfer and Wierzbicki 2015;Quinn and Chang 2015).
Nevertheless, there has been growing evidence that lncRNAs activate repressed genes by modulating local chromatin structures, which can facilitate coordinated gene expression in budding yeast (Cloutier et al. 2013(Cloutier et al. , 2016 and in other eukaryotes (Ørom et al. 2010;Wang et al. 2011b;Li et al. 2013;Boque-Sastre et al. 2015). In favor of this regulatory phenomenon, we hypothesize that the F. graminearum lncRNAs, which showed parallel induction with sense transcripts across the sexual development, may play a role in regulation of DNA synthesis and degradation at later developmental stages, according to "guilt-by-association" (Guttman et al. 2009;Gaiti et al. 2015). However, lncRNAs often exhibited cell type-specific expression (Bitton et al. 2011;Bumgarner et al. 2012;Li et al. 2013), and even allelespecific expression within the same nucleus (Rosa et al. 2016). Therefore, we cannot rule out the crassa. The expression of an lncRNA antisense to the circadian clock gene, frequency, promotes the sense transcription by generating a more transcriptionally permissive chromatin that has been silenced by sRNA-mediated heterochromatin formation (Li et al. 2015).
This study presents genome-wide characterization of lncRNAs during the fruiting body development. The transcriptome landscape including lncRNAs during the life cycle of F. graminearum provides fundamental genomic resources to the fungal community. The detailed molecular study of newly identified lncRNAs, with its established tools for rapid genetic analyses and ample genetic resources, will contribute to the understanding of how fungi utilize noncoding genomes for laying out their multicellular body plan.

Data generation and processing
The F. graminearum genome assembly (Cuomo et al. 2007)  32 (King et al. 2015) of a wild-type strain PH-1 (accessions: FGSC9075/NRRL31084) were used throughout this study. For total RNA extraction, synchronized fungal tissues were collected from carrot agar cultures at the previously defined developmental stages during perithecia formation (Trail and Common 2000;Hallen et al. 2007;Sikhakolli et al. 2012;Trail et al. 2017). For transcriptome data of spore germination stages, asexual spores (macroconidia) were spread on Bird agar medium (Metzenberg 2004) overlaid with a cellophane membrane, and sampled at the indicated spore germination stages (Supplemental Fig. S5). Strand-specific cDNA libraries were constructed from poly-A captured RNAs,  .

Differential expression and functional enrichment analyses
Read counts for gene loci were calculated using the htseq-count program (v0.6.1; Anders et al. 2015 Table S1). Functional enrichment analyses for DE genes were performed using the GOseq R package (v1.24.0), including only those genes annotated by one or more GO terms (Young et al. 2010). To assess enrichment of GO terms, the Wallenius approximation (an extension of the hypergeometric distribution) and Benjamini-Hochberg method were used to calculate the FDR-corrected P value (Supplemental Table S1).

Coexpression network analysis
The weighted gene correlation network analysis (WGCNA) R package (v1.51; Langfelder and correlation based on not just the direct correlation value of pairs of genes, but also the weighted correlations of all of their shared neighbors in the network space (Zhang and Horvath 2005). The softthresholding power 26 was selected, which is the lowest power for which the scale-free topology model fit index reaches 0.80. A range of treecut values were tested for cluster detection and the value was set to 0.18 (corresponding to correlation of 0.82). All other WGCNA parameters remained at their default settings.

Identification of small RNA clusters
The sRNA-seq and degradome-seq data were obtained from NCBI GEO (GSE87835) and NCBI SRA (PRJNA348145), respectively (Son et al. 2017). In filamentous fungi, the size of a majority of sRNAs ranges from 17 to 27 nt with a strong 5'U preference (70-82%; Lee et al. 2010;Son et al. 2017).
Subsequently, the number of 5'U-sRNA reads aligned to different genomic features (e.g. coding regions) were counted for each sRNA cluster, using the htseq-count program (v0.8.0; Anders et al. 2015). The degradome-seq dataset was processed according to the previous study (Son et al. 2017).

Expression correlation analysis
The expression value matrix for the 18 RNA-seq data (6 stages × 3 replicates) were rearranged by the BLIND program (Supplemental Fig. S2). Pearson's correlation and the associated P value by using an R script (Gaiti et al. 2015). For sense mRNAs only, that showed positive correlation with antisense lncRNAs, we performed functional enrichment analyses, using the GOseq R package (v1.24.0; Young et al. 2010), as did for DE genes between two successive developmental stages (see above).

XUT identification
Generation and confirmation of gene-deletion mutants and RNA-seq methodology for ∆xrn1 were described in Supplemental methods. For ∆xrn1 RNA-seq data, transcripts were independently assembled and merged to the de novo annotations, using the 'merge' function in the StringTie program (v1.3.0; Pertea et al. 2015), to incorporate novel transcripts that were only expressed in ∆xrn1 due to loss of NMD activity. Transcript abundance is globally affected by XRN1 deletion in budding yeast, however the expression of ribosomal protein genes shows only slight increases in ∆xrn1 and remains constant after lithium treatment that inhibits 5'→3' exonuclease activities (Dichtl et al. 1997;van Dijk et al. 2011).
Thus, before XUT identification, expression values were normalized in such a way that ribosomal protein genes are expressed at the same levels in WT and ∆xrn1 by multiplying 0.66 and 1.96 to RPKM values of WT at S0 and S4, respectively (Supplemental Fig. S11). Differential expression analyses were performed, using the 'stattest' function with an option argument: 'libadjust=FALSE', in the Ballgown R package The RNA-seq data generated in the present work have been deposited in NCBI's Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo), and are accessible through GEO series accession number GSE109095 that is composed of the two datasets, the sexual stage and ∆xrn1 dataset (GSE109094) and the spore germination stage dataset (GSE109088).

ACKNOWLEDGMENTS
We thank Brad Cavinder, Kevin Childs, and Nicholas Panchy (Michigan State Univ., MI, USA) for helpful discussion on analyzing data and for providing Perl and Python scripts. This work was supported by the National Science Foundation (Grant IOS-1456482 to F.T.) and Michigan AgBioResearch (F.T.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. performed the experiments; W.K., J.W., J.P.T. and F.T. analyzed the data; W.K. wrote the paper; and J.P.T. and F.T. revised the paper.

DISCLOSURE DECLARATION
The authors declare no conflicts of interest.