Long noncoding RNA repertoire and targeting by nuclear exosome, cytoplasmic exonuclease, and RNAi in fission yeast

Long noncoding RNAs (lncRNAs), which are longer than 200 nucleotides but often unstable, contribute a substantial and diverse portion to pervasive noncoding transcriptomes. Most lncRNAs are poorly annotated and understood, although several play important roles in gene regulation and diseases. Here we systematically uncover and analyze lncRNAs in Schizosaccharomyces pombe. Based on RNA-seq data from twelve RNA-processing mutants and nine physiological conditions, we identify 5775 novel lncRNAs, nearly 4× the previously annotated lncRNAs. The expression of most lncRNAs becomes strongly induced under the genetic and physiological perturbations, most notably during late meiosis. Most lncRNAs are cryptic and suppressed by three RNA-processing pathways: the nuclear exosome, cytoplasmic exonuclease, and RNAi. Double-mutant analyses reveal substantial coordination and redundancy among these pathways. We classify lncRNAs by their dominant pathway into cryptic unstable transcripts (CUTs), Xrn1-sensitive unstable transcripts (XUTs), and Dicer-sensitive unstable transcripts (DUTs). XUTs and DUTs are enriched for antisense lncRNAs, while CUTs are often bidirectional and actively translated. The cytoplasmic exonuclease, along with RNAi, dampens the expression of thousands of lncRNAs and mRNAs that become induced during meiosis. Antisense lncRNA expression mostly negatively correlates with sense mRNA expression in the physiological, but not the genetic conditions. Intergenic and bidirectional lncRNAs emerge from nucleosome-depleted regions, upstream of positioned nucleosomes. Our results highlight both similarities and differences to lncRNA regulation in budding yeast. This broad survey of the lncRNA repertoire and characteristics in S. pombe, and the interwoven regulatory pathways that target lncRNAs, provides a rich framework for their further functional analyses.


INTRODUCTION
Genomes are more pervasively transcribed than expected from their protein-coding sequences. For example, over 80% of the human genome is transcribed but only ∼2% encodes proteins (Hangauer et al. 2013). So-called long noncoding RNAs (lncRNAs), which exceed 200 nucleotides (nt) in length but lack long open reading frames, make up a sub- Derrien et al. 2012;Pauli et al. 2012;Hon et al. 2017). In general, lncRNAs are transcribed by RNA polymerase II and seem to be capped and polyadenylated (Guttman et al. 2009;Cabili et al. 2011;Derrien et al. 2012), although the patterns of transcription and RNA processing can radically differ between mRNAs and lncRNAs (Tuck and Tollervey 2013;Quinn and Chang 2016;Mukherjee et al. 2017;Schlackow et al. 2017). Some lncRNAs engage with ribosomes, which can trigger nonsense-mediated decay (NMD) to dampen their expression but may also produce functional peptides in some cases (Malabat et al. 2015;Quinn and Chang 2016;Wery et al. 2016;de Andres-Pablo et al. 2017).
Given the profusion, diversity and low expression of lncRNAs, their full description and annotation is still ongoing and evolving (Atkinson et al. 2012;Mattick and Rinn 2015;St. Laurent et al. 2015). A conceptually simple way to classify lncRNA genes is by their position relative to neighboring coding genes. For example, long intervening noncoding RNAs (lincRNAs), transcribed from intergenic regions that do not overlap any mRNAs, have been the subject of much research in mammalian cells (Rinn and Chang 2012;Ulitsky and Bartel 2013;Schlackow et al. 2017). Antisense lncRNAs are transcribed in the opposite direction to mRNAs with which they completely or partially overlap; they can affect sense transcript levels via diverse mechanisms (Pelechano and Steinmetz 2013;Mellor et al. 2016). Bidirectional lncRNAs, on the other hand, emerge close to the transcriptional start site of coding genes but run in the opposite direction; most eukaryotic promoters initiate divergent transcription leading to widespread bidirectional lncRNAs, although transcriptional elongation is often only productive in the sense direction (Grzechnik et al. 2014;Quinn and Chang 2016). Bidirectional transcription has been proposed to drive the origination of new genes (Wu and Sharp 2013) and modulate gene-expression noise (Wang et al. 2011).
Pervasive transcription is potentially harmful as it can affect the expression of coding genes (Mellor et al. 2016), and nascent RNAs can compromise genome stability (Li and Manley 2006). Cells therefore apply RNA surveillance systems to keep the expression of lncRNAs in check . Many lncRNAs are actively degraded shortly after transcription, suggesting that they could reflect transcriptional noise, byproducts of coding transcription, or that the act of transcription, rather than the lncRNA itself, is functionally important . Different RNA-decay pathways preferentially target distinct sets of lncRNAs, which has been used for lncRNA classification in budding yeast. The "cryptic unstable transcripts" (CUTs) accumulate in cells lacking Rrp6 (Neil et al. 2009;Xu et al. 2009), a nuclear-specific catalytic subunit of the RNA exosome (Houseley and Tollervey 2009;Kilchert et al. 2016). Conversely, the "stable unannotated transcripts" (SUTs) are less affected by Rrp6 (Neil et al. 2009;Xu et al. 2009). The human "promoter upstream transcripts" (PROMPTs) are analogous to CUTs (Preker et al. 2008). The "Xrn1-sensitive unstable transcripts" (XUTs) accumulate in cells lacking the cytoplasmic 5 ′ exonuclease Xrn1 and are often antisense to mRNAs (Houseley and Tollervey 2009;van Dijk et al. 2011;Wery et al. 2018). "Nrd1-unterminated transcripts" (NUTs) and "Reb1p-dependent unstable transcripts" (RUTs) are controlled via different mechanisms of transcriptional termination linked to RNA degradation (Schulz et al. 2013;Colin et al. 2014). These classes, while useful, are somewhat arbitrary and overlapping as most lncRNAs can be targeted by different pathways, especially if one pathway is absent in mutant cells .
The fission yeast Schizosaccharomyces pombe provides a potent complementary model system to study gene regulation. In some respects, RNA metabolism of fission yeast is more similar to metazoan cells than budding yeast. For example, RNA interference (RNAi) (Castel and Martienssen 2013), RNA uridylation (Schmidt et al. 2011), alternative-polyadenylation features (Liu et al. 2017), and Pab2/PABPN1-dependent RNA degradation (Lemay et al. 2010;Lemieux et al. 2011;Beaulieu et al. 2012) are conserved from fission yeast to humans, but not in budding yeast. Several genome-wide studies have uncovered widespread lncRNAs (Dutrow et al. 2008;Wilhelm et al. 2008;Rhind et al. 2011;Eser et al. 2016), and over 1500 lncRNAs are currently annotated in the PomBase model organism database (McDowall et al. 2014). Studies with natural isolates of S. pombe revealed that only the most highly expressed lncRNAs show purifying selection (Jeffares et al. 2015), but the regulation of many lncRNAs is affected by expression quantitative trait loci (Clément-Ziza et al. 2014). Over 85% of the annotated S. pombe lncRNAs are expressed below one copy per cell during proliferation, and over 97% appear to be polyadenylated . Ribosome profiling showed that as many as 24% of lncRNAs are actively translated . As in other organisms, a large proportion of the S. pombe lncRNAs are antisense to mRNAs and have been implicated in controlling the meiotic gene expression program (Ni et al. 2010;Bitton et al. 2011;Chen et al. 2012;Wery et al. 2018). Diverse chromatin factors function in suppressing antisense and other lncRNAs in S. pombe, including the HIRA histone chaperone (Anderson et al. 2009), the histone variant H2A.Z (Zofall et al. 2009;Clément-Ziza et al. 2014), the Clr4/Suv39 methyltransferase together with RNAi (Zhang et al. 2011), the Spt6 histone chaperone (DeGennaro et al. 2013;Wery et al. 2018), and the CHD1 chromatin remodeler (Hennig et al. 2012;Pointner et al. 2012;Shim et al. 2012). Analogous to CUTs in budding yeast, different types of RNAs accumulating in rrp6 mutants have also been described (Zhou et al. 2015). Several lncRNAs have been functionally characterized in S. pombe: meiRNA controls meiotic differentiation and chromosome pairing (Ding et al. 2012;Yamashita et al. 2016), stress-induced lncRNAs activate expression of the downstream fbp1 gene during glucose starvation (Hirota et al. 2008;Oda et al. 2015), adh1AS is an antisense inhibitor of the adh1 gene during zinc limitation (Ehrensberger et al. 2013), prt recruits the exosome to control phosphate-tuned pho1 expression (Ard et al. 2014;Shah et al. 2014), nc-tgp1 inhibits the phosphate-responsive permease tgp1 gene by transcriptional interference (Ard et al. 2014), nam1 regulates entry into meiotic differentiation (Touat-Todeschini et al. 2017), and SPNCRNA.1164 regulates expression of the atf1 transcription-factor gene in trans during oxidative stress (Leong et al. 2014). Naturally, these studies only scratch the surface, with the noncoding transcriptome and any of its functions remaining relatively poorly defined in fission yeast and other organisms.
Transcriptome analyses under selective conditions, such as in RNA-processing mutants, have proven useful to define lncRNAs in budding yeast. Here we analyze transcriptome sequencing under multiple genetic and physiological perturbations in fission yeast to maximize the detection and initial characterization of lncRNAs. Some of these RNA-seq samples have previously been analyzed with respect to mRNA processing and expression (Lemieux et al. 2011;Marguerat et al. 2012;Schlackow et al. 2013;Lemay et al. 2014;Bitton et al. 2015a). They interrogate pathways, such as RNAi and Pab2/PABPN1, which are conserved in humans but not in budding yeast. We identify 5775 novel, unannotated lncRNAs, in addition to the previously annotated lncRNAs. The expression of lncRNAs is more extensively regulated in stationary phase, quiescence and, most notably, meiotic differentiation than the expression of mRNAs. Many lncRNAs comprise unstable transcripts that are repressed by three partially overlapping RNA-processing pathways. Analogous to budding yeast, we classify the unstable lncRNAs targeted by Rrp6, Dcr1, and Exo2 into CUTs, DUTs, and XUTs, respectively. We further analyze the positions and expression of all novel and annotated lncRNAs with respect to neighboring mRNAs, and other biological characteristics such as translation and nucleosome patterns. Both similarities and notable differences to lncRNAs in budding yeast are discussed. This extensive study provides a framework for functional characterization of lncRNAs in fission yeast and beyond.

Detection of novel lncRNAs
To broadly identify lncRNAs in fission yeast, we examined strand-specific RNA-seq data acquired under multiple genetic and physiological conditions. We analyzed the transcriptomes of twelve RNA-processing mutants to facilitate detection of RNAs that may be rapidly degraded (Supplemental Table  S1). This mutant panel affects proteins for key pathways of RNA processing and degradation: Rrp6, a 3 ′ -5 ′ exonuclease of the nuclear RNA exosome (Harigaya et al. 2006;Lemay et al. 2014); Dis3, a 3 ′ -5 ′ exo/endonuclease of the core RNA exosome (Wang et al. 2008); Ago1 (Argonaute), Dcr1 (RNase III-like Dicer), and Rdp1 (RNA-dependent RNA polymerase) of the RNAi pathway (Volpe et al. 2002); Exo2, a cytoplasmic 5 ′ exonuclease (ortholog of XRN1) (Houseley and Tollervey 2009); Ski7, a cytoplasmic cofactor which links the Ski complex to the exosome (Lemay et al. 2010); Cid14, a poly(A) polymerase of the TRAMP complex which is a cofactor of the nuclear RNA exosome (Wang et al. 2008); Pab2, a poly(A)-binding protein targeting RNAs to the nuclear exosome (PABPN1 ortholog) (Lemieux et al. 2011); Pan2, a deadenylase of the Pan2-Pan3 complex (Wolf and Passmore 2014); and Upf1, an ATP-dependent RNA helicase of the NMD pathway (Rodríguez-Gabriel et al. 2006). We also analyzed transcriptomes under nine physiological conditions to sample key cellular states (Supplemental Table S2): two timepoints of stationary phase after glucose depletion (100% and 50% cell viability); two timepoints of quiescence/cellular ageing (days 1 and 7 after nitrogen removal); and five timepoints of meiotic differentiation (0 to 8 h after triggering meiosis).
We designed a simple segmentation heuristic to detect novel lncRNAs longer than 200 nt, favoring sensitivity at the cost of specificity (Materials and Methods). Applying this approach to the RNA-seq data covering the genetic and physiological perturbations, we identified 5775 novel, unannotated lncRNAs in addition to the ∼1550 previously annotated lncRNAs. We assigned systematic names, SPNCRNA.2000 to SPNCRNA.7774, to these novel lncRNAs (Supplemental Table S3). Of these unannotated lncRNAs, 159 fully or partially overlapped on the same strand with 167 of the 487 novel lncRNAs recently reported (Supplemental Table S4; Eser et al. 2016). On the other hand, 214 of these 487 lncRNAs could not be validated using our data (Supplemental Table  S4). In our samples, these 214 lncRNAs typically showed no sequencing signals in the sense strand, but strong signals on the opposite strand for these regions (Materials and Methods).
Naturally, annotation of weakly expressed lncRNAs based on RNA-seq data, and especially the locations of their TSS and TTS, is somewhat arbitrary and should be understood as an approximation. Nevertheless, our data reveal many more S. pombe lncRNAs than have been previously recognized, and the broad annotations provide a resource for further studies. We defined three confidence classes for all annotated and novel lncRNAs, based on RPKM values from the sequenced samples (Materials and Methods). This approach resulted in 2090 high-confidence lncRNAs (1064 annotated, 1026 novel), 5004 medium-confidence lncRNAs (452 annotated, 4552 novel), and 254 low-confidence lncRNAs (57 annotated, 197 novel). So most novel lncRNAs are medium-confidence, while most annotated lncRNAs are high-confidence using our criteria. Supplemental Table S5 provides the confidence class for all annotated and novel lncRNAs. Examples of novel lncRNAs are provided in a browser view in Supplemental Figure S1. This figure suggests that some lncRNAs might consist of transcribed regions encompassing several smaller RNAs or variable, condition-specific isoforms. To assess these complexities for specific lncRNAs, we also provide a web tool to view RNA-seq data for all lncRNAs and mRNAs in the different conditions (http:// bahlerlab.info/ncViewer). Our data will also be included in a specific track of the JBrowse genome brower in PomBase.
The novel lncRNAs were generally shorter (mean 797 nt) than the annotated lncRNAs (mean 1233 nt) and mRNAs (mean 2148 nt) (Supplemental Fig. S2a). While some sequence library protocols can generate spurious antisense RNAs (Perocchi et al. 2007), our protocol is resilient to this artifact as it relies on ligating two RNA oligonucleotides to fragmented mRNAs. We mostly analyzed poly(A)-enriched samples, because almost all S. pombe lncRNAs are polyadenylated . As a control, we compared the data to samples depleted for ribosomal RNA (rRNA) in rrp6Δ and exo2Δ mutants, which confirmed that the great majority of lncRNAs are polyadenylated (Supplemental Fig. S2b). We cannot exclude the possibility, however, that some nonpolyadenylated lncRNAs are missing from our annotation.
None of the novel lncRNAs overlapped with mRNAs on the same strand using current PomBase annotations, and only 21 overlapped using CAGE (Li et al. 2015) and poly(A) data (Mata 2013), respectively, for transcription start and termination sites (Supplemental Table S5). Thus, the novel lncRNAs do not represent alternative transcription start or termination sites of known mRNAs. On the other hand, 3650 of 5138 protein-coding regions (71%) overlapped by at least 10 nt with antisense lncRNAs, either annotated or novel. The 1461 (28.4%) of coding regions not associated with antisense lncRNAs were enriched for the 20% shortest mRNAs (P ∼ 1.6 × 10 −16 ), including those encoding ribosomal proteins (P ∼ 1.1 × 10 −5 ), as well as for several features associated with high gene expression, including high mRNA levels (P ∼ 0.004) (Pancaldi et al. 2010), stable mRNAs (P ∼ 2.1 × 10 −9 ) (Hasan et al. 2014), and mRNAs that show high RNA polymerase II occupancy (P ∼ 2.4 × 10 −5 ) and high ribosome density (P ∼ 1.2 × 10 −41 ) (Lackner et al. 2007). These enrichments suggest that highly expressed genes are either protected from antisense transcripts or interfere with antisense transcription. Figure 1 shows the relative expression changes of all mRNAs, annotated lncRNAs and novel lncRNAs in the genetic and physiological conditions. For all conditions, differential expression was determined relative to three reference samples (exponentially proliferating wild-type and auxotrophic control cells in minimal medium with supplements; Supplemental Table S2). This panel of reference samples included different genetic markers and supplements used for the other conditions. Only 403 RNAs, including 184 novel lncRNAs, were differentially expressed between wild-type and auxotrophic control cells grown in rich yeast extract (YE) medium compared to minimal media (>2-fold change, P < 0.05). Thus, the growth medium and auxtrophic markers only minimally affected lncRNA expression. All differential expression data are available in Supplemental Table S3.
About 50% of the mRNAs were differentially expressed, both induced and repressed, in the different physiological A C B FIGURE 1. Hierarchical clustering of gene expression in different RNA-processing mutants and physiological conditions. Expression profiles are shown for (A) all 5177 mRNAs, (B) 1573 annotated lncRNAs, and (C) 5775 novel, unannotated lncRNAs. Changes in RNA levels in response to the different genetic and physiological conditions (indicated at bottom) relative to control cells grown in minimal medium are colorcoded as shown in the color legend at bottom right (log 2 expression ratios). The rrp6 and exo2 samples indicated by asterisks have been depleted for rRNA instead of poly(A) purification used for all other samples. Data obtained from these two types of RNA-seq libraries are compared in Supplemental Figure S2b. Details on strains and conditions are provided in Supplemental Tables S1 and S2, and all expression data are provided in Supplemental Table S3. conditions, but much less so in the different genetic conditions (Figs. 1A, 2). The mRNAs up-regulated in the exosome, pab2, and RNAi mutants were enriched for Gene Ontology (GO) terms related to meiosis, consistent with reported roles of the corresponding proteins in meiotic gene silencing (St-André et al. 2010;Yamanaka et al. 2010Yamanaka et al. , 2013. Notably, mRNAs up-regulated in the exo2 mutant were strongly and specifically enriched for middle meiotic genes (Mata et al. 2002). Figure 1B and C show the expression of the previously annotated and novel lncRNAs, respectively. Compared to mRNAs, much higher proportions and numbers of lncRNAs featured strong differential expression, mostly induced, under the different genetic and physiological conditions (Figs. 1B,C, 2). The following mutants led to the most pronounced effects on lncRNA expression: nuclear exosome (rrp6-ts, rrp6Δ, dis3-54), RNAi (ago1Δ, dcr1Δ, rdp1Δ), and the cytoplasmic exonuclease (exo2Δ). In nuclear exosome mutants, most lncRNAs were strongly derepressed, in stark contrast to the mRNAs which showed a higher proportion of repressed transcripts in this condition (Figs. 2,3A). Many of the novel lncRNAs in particular were also strongly derepressed in RNAi and exo2 mutants and in meiotic cells (Fig. 3A), most notably during late meiotic stages (Mei 6-8 h; Figs. 1C, 2). On the other hand, the novel lncRNAs showed generally lower expression levels in both genetic and physiological conditions than the annotated lncRNAs, and much lower than the mRNAs (Fig. 3B). Remarkably, only 54 novel lncRNAs were neither up-nor down-regulated in any of the conditions (>2-fold change, P < 0.05). Moreover, based on recent NET-seq data (Shetty et al. 2017), we estimated that ∼87% of the novel lncRNA genes are actively transcribed in proliferating or Spt5-depleted cells (Materials and Methods). Such active transcription included genes for which we could not detect the corresponding lncRNAs by RNA-seq under the standard condition. These results argue against sequencing FIGURE 2. Histograms showing numbers and proportions of induced (orange), repressed (blue), and all other (gray) transcripts for the different RNA classes as indicated. Differentially expressed genes were defined as those being ≥2-fold induced (mean of two biological repeats) or repressed and showing significant changes (P < 0.05) compared to reference as determined by DESeq2. A B FIGURE 3. Gene expression in major groups of genetic and environmental conditions. (A) (Left graph) Box plot of expression ratios (condition relative to control) of all mRNAs, annotated and novel lncRNAs in nuclear exosome (rrp6Δ, rrp6-ts), RNAi (ago1Δ, dcr1Δ, rdp1Δ), and cytoplasmic exonuclease (exo2Δ) mutants. (Right graph) As left but for quiescence, stationary phase, and meiosis conditions. The horizontal dashed lines indicate twofold induction and repression. (B) As in panel A, but for expression levels (RPKM scores). All expression data are provided in Supplemental Table S3. artifacts and support that the turnover of these transcripts is modulated in response to different biological conditions. Together, their low expression levels and induction under specialized conditions can explain why the novel lncRNAs have not been identified in previous studies (Wilhelm et al. 2008;Rhind et al. 2011;Eser et al. 2016).
In conclusion, these findings indicate that the novel lncRNAs comprise many cryptic transcripts that are degraded by three main RNA-processing pathways during mitotic proliferation: the nuclear exosome, the RNAi machinery, and the cytoplasmic exonuclease Exo2. The other RNA processing factors analyzed here seem to play only minor or redundant roles in lncRNA regulation. Many lncRNAs become induced under specific physiological conditions when they might play specialized roles.

Classification of lncRNAs into CUTs, DUTs, and XUTs
In budding yeast, different groups of lncRNAs have been named according to the RNA-processing pathways controlling their expression. For example, CUTs are targeted for degradation by the nuclear exosome (Neil et al. 2009;Xu et al. 2009), and XUTs are targeted by the cytoplasmic exonuclease Xrn1 (ortholog of S. pombe Exo2) (van Dijk et al. 2011). We introduce an analogous classification of lncRNAs to provide a framework for analysis. Based on the panel of RNA-processing mutants tested here, the nuclear exosome, the cytoplasmic exonuclease, and the RNAi machinery are the three main pathways targeting lncRNAs (Figs. 1, 3). These three pathways thus provide a natural way to classify the lncRNAs as CUTs and XUTs (corresponding to the budding yeast classes of the same names) and DUTs (Dicer-sensitive unstable transcripts). DUTs define a novel class not applicable to budding yeast which lacks the RNAi machinery.
Using a fuzzy clustering approach (Materials and Methods), we classified both the novel and previously annotated lncRNAs that were significantly derepressed in at least one of the following mutants: rrp6Δ (CUTs), dcr1Δ (DUTs), or exo2Δ (XUTs). Of 7308 lncRNAs, 2068 remained unclassified because they were not significantly derepressed in any of the three mutants (1896 lncRNAs) or could not be assigned to a single class (172 lncRNAs). The remaining lncRNAs were classified into 2732 CUTs (493 annotated, 2239 novel), 1116 XUTs (181 annotated, 935 novel), and 1392 DUTs (209 annotated, 1183 novel). The resulting three classes were quite distinct (Fig.  4A). These class associations are provided in Supplemental  Tables S3 and S5. As expected, the lncRNAs of a given class were most highly expressed on average in the mutant used to define the class, although they were also more highly expressed in the other two mutants than in wild-type cells (Fig. 4B). This pattern was also evident when clustering the three lncRNA classes separately based on the expression changes in the different RNA-processing mutants: lncRNAs of a given class were most highly derepressed in the mutant used to define this class, and in mutants affecting the same pathway, but they also tended to be derepressed in mutants of other pathways (Fig. 4C). These results show that the lncRNAs of a given class are not B A C FIGURE 4. Classification of lncRNAs into CUTs, DUTs, and XUTs. (A) The lncRNAs significantly induced (DESeq2; adjusted P ≤ 0.05) in rrp6Δ, dcr1Δ, or exo2Δ mutants were clustered into CUTs, DUTs, and XUTs, respectively, using the Mfuzz R package (default parameters, three clusters specified). The clustering shows the 5586 uniquely classified lncRNAs after filtering those with a membership score <0.7. The red/blue colors indicate the mean RPKM values in the three mutants as indicated, scaled by subtracting the mean of the row and division by the standard deviation of the row (z-score). The assigned clusters are indicated at left. (B) Box plots of RPKM values (log 2 ) of all CUTs (left), DUTs (middle), and XUTs (right) in control (MM) and mutant cells as indicated. The data for the biological repeats 1 and 2 are plotted separately. (C) Hierarchical clustering of genetic conditions for all CUTs, DUTs, and XUTs as indicated. Changes in RNA levels in response to the different genetic conditions (indicated at bottom) relative to control cells grown in minimal medium are color-coded as shown in the legend at bottom right (log 2 expression ratios). All expression and classification data are provided in Supplemental Table S3. exclusively derepressed in the mutant used to define this class. Thus, lncRNAs can be degraded by different pathways, although one pathway is typically dominant for a given lncRNA which is used here for classification.

Relationships between RNA-processing pathways targeting lncRNAs
To further dissect the functional relationships among the nuclear exosome, RNAi, and cytoplasmic exonuclease pathways, we attempted to construct double and triple mutants of rrp6Δ, dcr1Δ, and exo2Δ. We could not obtain an rrp6Δ exo2Δ mutant from 24 tetrads dissected from the corresponding cross. Among the tetrads analyzed with three viable spores, the nonsurviving spore was always of the rrp6Δ exo2Δ genotype, with a significant difference between observed and expected frequencies of wild-type, single, and double mutant spores (χ 2 test, P ∼ 10 −4 ; based on 61% surviving spores). We conclude that the rrp6Δ exo2Δ double mutant is not viable. This synthetic lethality indicates that the nuclear exosome and cytoplasmic exonuclease together exert an essential role.
Conversely, the exo2Δ dcr1Δ and rrp6Δ dcr1Δ double mutants were viable, although the latter showed stronger growth defects than either single mutant (Supplemental Fig. S3). While the exo2Δ dcr1Δ cells were elongated like the exo2Δ cells (Szankasi and Smith 1996), the rrp6Δ and rrp6Δ dcr1Δ cells were of normal length, with the double mutant looking more irregular and sick (Supplemental Fig. S3). These findings suggest that the nuclear exosome and RNAi machineries can back each other up to some extent. We also attempted to construct an rrp6Δ exo2Δ dcr1Δ triple mutant by mating of the exo2Δ dcr1Δ and rrp6Δ dcr1Δ double mutants. This mating only produced ∼17% viable spores, suggesting that the RNAi machinery is required for spore survival. As expected, we did not obtain any triple mutant among 24 tetrads dissected from this cross (χ 2 test, P ∼ 10 −4 ; based on 17% surviving spores). Thus, both the rrp6Δ exo2Δ double mutant and the rrp6Δ exo2Δ dcr1Δ triple mutant are not viable. Consistent with lncRNAs being targeted by multiple pathways (Fig. 4B,C), these results point to some redundancy in function between the different RNA degradation pathways (Discussion).
We analyzed the transcriptomes of the exo2Δ dcr1Δ and rrp6Δ dcr1Δ double mutants by RNA-seq. The rrp6Δ dcr1Δ double mutant showed a greater number of derepressed lncRNAs than either single mutant, especially among the novel lncRNAs (Figs. 2, 4C). The 5647 lncRNAs that were significantly derepressed in the rrp6Δ dcr1Δ double mutant (expression ratio >2, P < 0.05) included 2653 CUTs and 1317 DUTs, but also 884 XUTs and 793 unclassified lncRNAs. This result again shows that the nuclear exosome and RNAi have partially redundant roles and can back each other up with respect to many RNA targets. Moreover, these two nuclear pathways can also degrade most XUTs that are further targeted by the cytoplasmic exonuclease. These findings highlight a prominent role of the joint activity of the nuclear exosome and RNAi to suppress a large number of lncRNAs.
In contrast, the exo2Δ dcr1Δ double mutant showed fewer derepressed XUTs and DUTs than either single mutant (Figs. 2,4C). The 2272 lncRNAs that were significantly derepressed in the double mutant included only 675 XUTs and 658 DUTs, but 786 CUTs and 153 unclassified lncRNAs. This genetic suppression is consistent with the slightly diminished growth defect of the exo2Δ dcr1Δ double mutant compared to the exo2Δ single mutant (Supplemental Fig. S3). Taken together, the complex findings from single and double mutants indicate that the three RNA-processing pathways can all target a wide range of lncRNAs, with differential preference for some targets, and channeling of lncRNAs to alternative pathways in the mutants. The higher numbers of derepressed lncRNAs in rrp6Δ and dcr1Δ mutants (Figs. 2, 4C) indicate that most lncRNAs are degraded in the nucleus, most notably by the nuclear exosome.

Classification of lncRNAs by neighboring mRNA positions
We also classified all known and novel lncRNAs into the main types based on their positions relative to the neighboring mRNAs: bidirectional, antisense, and intergenic (Fig.  5A). The criteria for these assignments, and overlaps between different classes, are specified in Materials and Methods. In total, we defined 1577 Bidirectional (539 annotated, 1038 novel), 4474 Antisense (575 annotated, 3899 novel), and 1189 Intergenic (356 annotated, 833 novel) lncRNAs, besides 108 that overlapped mRNAs in sense direction (103 annotated, 5 novel; reflecting recent changes in TSS and TTS annotations in PomBase). Data for these lncRNA classes are provided in Supplemental Table S5. The bidirectional lncRNAs were enriched for CUTs (Fig. 5A). The Antisense lncRNAs were enriched for XUTs and, most notably, for DUTs, while the Intergenic lncRNAs were enriched for other, not classified lncRNAs (Fig. 5A). Similar trends were also evident when analyzing the lncRNA types the other way round (Fig. 5B): While most CUTs, DUTs, and XUTs were antisense, only the DUTs and XUTs were enriched for antisense lncRNAs, while the CUTs were enriched for bidirectional lncRNAs.

Expression patterns of lncRNAs
The expression of lncRNAs could be affected by the expression of neighboring mRNAs, or it could actively control the expression of neighboring mRNAs. To analyze the relationship between Bidirectional and Antisense lncRNAs and their associated mRNAs, we calculated the correlation coefficients for expression levels of each lncRNA-mRNA pair across all genetic and physiological conditions. Figure 5C shows that the expression of the Bidirectional lncRNAs tended to positively correlate with the expression of their associated mRNAs, for both the genetic and physiological conditions. The Antisense lncRNAs, on the other hand, revealed substantial differences for genetic versus physiological conditions. They mainly correlated negatively with sense mRNA expression in the physiological, but not in genetic conditions (Fig. 5C). Thus, accumulation of Antisense lncRNAs in the different RNA-processing mutants is generally not sufficient to repress mRNA levels. These striking contrasts between the two lncRNA classes and different types of conditions provide clues about the lncRNA-mRNA expression relationships for Bidirectional and Antisense lncRNAs (Discussion).
The up-regulation of lncRNAs under specific physiological conditions (Figs. 1, 3A) could reflect specialized roles under these conditions. We checked whether the different lncRNAs showed distinct regulatory patterns by clustering the classes separately based on their expression changes in the different physiological conditions. When clustering CUTs, DUTs, and XUTs, the physiological conditions always clustered into three main groups that showed distinct patterns of lncRNA expression (Supplemental Fig.  S4): late meiosis, stationary phase (triggered by glucose limitation), and quiescence/early meiosis (both triggered by nitrogen limitation). These three groups reflect the major physiological states interrogated by the different conditions. The expression changes of CUTs, DUTs, and XUTs, however, did not substantially differ across the physiological conditions. Most lncRNAs in all three classes, and among novel lncRNAs in particular, were strongly induced in late meiosis (Mei 8 h; Fig. 2; Supplemental  Fig. S4). Many lncRNAs were also induced during stationary phase and quiescence/early meiosis, with CUTs being relatively most conspicuous in these conditions ( Fig. 2; Supplemental Fig. S4).
We also analyzed the lncRNAs classified by their positions relative to neighboring mRNAs. As expected from Figure 5A,B, the nuclear exosome mutants showed the highest proportion of derepressed Bidirectional lncRNAs but were also prominent in derepressing Antisense and Intergenic lncRNAs (Fig. 6A). The rrp6Δ dcr1Δ double mutant led to derepression of a large number of lncRNAs, particularly Intergenic and Antisense lncRNAs (Fig. 6A). Among the physiological conditions, late meiosis (Mei 8h) showed the highest proportions of induced lncRNAs for all three classes; this response was most notable for Antisense lncRNAs, many of which were highly induced, while much fewer Intergenic lncRNAs were induced (Fig.  6A,B). Overall, the Bidirectional, Antisense and Intergenic lncRNAs showed stronger class-specific expression signatures in the physiological conditions than did the CUTs, XUTs, and DUTs ( Fig. 6B versus Supplemental Fig. S4). But expression patterns across the different physiological conditions were not sufficiently distinct to predict class membership based on these patterns.

Nucleosome profiles of lncRNA regions
Transcribed regions are often accompanied by distinct patterns of nucleosome distribution. As a proxy for nucleosome distributions around lncRNA regions, we sequenced A B C FIGURE 5. Analyses of lncRNAs by positions relative to mRNAs. (A) We grouped the annotated and novel lncRNAs into three main positional types as represented schematically: 1577 bidirectional, 4474 antisense, and 1189 intergenic RNAs, leaving only 108 lncRNAs (105 annotated, three novel) that overlapped mRNAs in sense direction (Materials and Methods). Pie charts of the corresponding proportions of CUTs, DUTs, XUTs, and other lncRNAs are provided beneath each positional type, and also for all annotated and known lncRNAs. Significantly enriched slices are indicated with asterisks (R prop.test function, P < 10 −6 ). (B) Pie charts of the proportions of bidirectional, antisense, intergenic, and sense lncRNAs for all (annotated and known) lncRNAs, and among the CUTs, DUTs, and XUTs. Significantly enriched slices are indicated with asterisks (R prop.test function, P < 10 −4 ). (C) Pearson correlation coefficients for RPKM expression data of each bidirectional lncRNA-mRNA pair (left) and each antisense lncRNA-mRNA pair (right). The correlation data are shown separately for all genetic and physiological conditions as indicated on x-axes. For antisense lncRNAs, the difference between the distributions in genetic versus physiological conditions is highly significant (P Wilcoxon = 4.6 × 10 −170 ). All classification data are provided in Supplemental Table S5. mononucleosomal DNA to compare the chromatin organization between protein-coding and noncoding transcribed regions in proliferating wild-type cells. We only analyzed lncRNA regions that did not overlap with any mRNA regions to minimize confounding data from coding transcription (although some mRNA transcription start sites may remain among the Bidirectional lncRNA data). Figure 7 shows that both Bidirectional and Intergenic lncRNAs initiated in nucleosome-depleted regions, at the 5 ′ -end of a positioned nucleosome. This feature was shared with mRNA regions, although these regions showed much higher nucleosome densities and a higher order of subsequent nucleosomes. In general, intergenic regions showed lower nucleosome densities than did mRNA regions, as has been observed before (Lantermann et al. 2010). We conclude that there are both similarities (around the transcription start site) and  in proliferating wild-type cells. Bidirectional lncRNA loci that overlap mRNAs in antisense direction were not included. We omitted 70 intergenic lncRNA loci that showed unusually high histone occupancies (Supplemental Table S8); these loci are mostly located in centromeric or subtelomeric regions. The top graphs show average nucleosome profiles for the different types of transcribed regions, aligned to the transcription start sites (TSS). The lower graphs show heatmaps for the first two kilobases of all transcribed regions analyzed, ordered by transcript length from top (longest RNAs) to bottom (shortest RNAs). Sequencing scores are color-coded as shown in legend at bottom right. differences in chromatin organization between coding and noncoding regions.

Translation of lncRNAs
Ribosome profiling before and during meiotic differentiation has revealed that as much as 24% of the annotated lncRNAs are actively translated . Such translation typically increases during meiosis and involves short open reading frames (ORFs), often more than one per lncRNA. In addition, translation has been detected in numerous unannotated regions of the genome . We analyzed the ribosome-profiling data from Duncan and Mata (2014), covering the whole genome of proliferating and meiotic cells, to assess translation of the different classes of annotated and novel lncRNAs defined above. Details of all 771 translated noncoding regions are provided in Supplemental Table S6. Table 1 shows the numbers and percentages of actively translated lncRNAs, both for all translated regions of at least one codon and for a conservative set of translated regions with at least ten codons. Overall, the novel lncRNAs showed a much lower proportion of translated transcripts than the annotated lncRNAs. This result likely reflects their lower expression levels (Fig. 3B). Such low expression makes it harder to obtain sufficient ribosome-profiling reads to determine translation for most lncRNAs, especially because no ribosome profiling data were available for most conditions in which the novel lncRNAs became derepressed. Nevertheless, 66 or 148 novel lncRNAs were found to be actively translated using the more or less conservative cutoff, respectively. Among the main classes, the Bidirectional lncRNAs showed the largest numbers and proportions of translated RNAs, which were highly enriched (Table 1; P = 1.7 × 10 −26 ). Moreover, up to 30% of the 108 previously annotated lncRNAs overlapping mRNAs in sense direction were actively translated. This finding suggests that some of these RNAs are alternative "untranslated regions" of the mRNAs with translation of short upstream ORFs (Duncan and Mata 2014).
It is not clear to what extent any stable, functional peptides are generated by all this translational activity. The engagement of lncRNAs with ribosomes could trigger NMD and degradation via the cytoplasmic exonuclease (Malabat et al. 2015;Quinn and Chang 2016;de Andres-Pablo et al. 2017). We therefore checked whether the translated lncRNAs were enriched among those being derepressed in the upf1Δ mutant which is defective for NMD (Rodríguez-Gabriel et al. 2006). There were no significant overlaps between the translated RNAs and the novel, annotated or all lncRNAs derepressed in upf1Δ cells. Moreover, despite including several translated lncRNAs, the Bidirectional lncRNAs showed no significant overlap with RNAs derepressed in upf1Δ cells. This finding is consistent with Bidirectional lncRNAs being mainly targeted by the nuclear exosome rather than the cytoplasmic exonuclease (Fig. 5A,B). Conversely, the translated XUTs and Antisense lncRNAs were both significantly enriched for RNAs derepressed in upf1Δ cells (P = 6.7 × 10 −24 and 1.5 × 10 −25 , respectively). This result is consistent with the cytoplasmic exonuclease being the major pathway for NMD-mediated RNA degradation. Together, these findings suggest that engaging with ribosomes triggers the degradation of some XUTs and Antisense lncRNAs, but not of Bidirectional lncRNAs which are thus more likely to produce peptides.

DISCUSSION
This study has uncovered 5775 novel lncRNAs, 1026 of which are high confidence based on sequence-read coverage. Our annotation is optimized for sensitivity at the cost of specificity. From our web viewer (http://bahlerlab.info/ncViewer), it is clear that some lncRNAs might consist of transcribed regions encompassing several smaller RNAs or variable, condition-specific isoforms. Annotations of lncRNAs are currently approximate in essentially all systems analyzed, and RNA-seq data have clear limitations for mapping of TSS and TTS. Nevertheless, the broad, approximate annotations of lncRNAs (or transcribed regions) under diverse conditions are a useful framework for future research on lncRNA function.
Compared to mRNAs and previously annotated lncRNAs, the novel lncRNAs are subject to stronger and more widespread differential expression, mostly induction, in response to multiple genetic and physiological conditions. Analysis of lncRNA expression across a broad panel of RNA-processing mutants indicates that the nuclear exosome, the RNAi pathway and the cytoplasmic exonuclease are key pathways targeting lncRNAs. Analogous to budding yeast, we have classified lncRNAs into CUTs, DUTs, and XUTs, defined by the pathway which preferentially degrades them. Notably, mRNAs are much less affected by the absence of these RNA-processing pathways than are lncRNAs (Figs. 1-3), suggesting that lncRNAs are important targets of these pathways. Unstable lncRNAs have been most extensively described at a genome- wide level in budding yeast which guided our classification. The lncRNA classes defined here show both similarities and differences to the classes defined in budding yeast, as highlighted below. The ∼2000 budding yeast CUTs are derepressed upon deletion of the nuclear-specific exosome subunit Rrp6 and transcribed divergently from mRNAs with which they positively correlate in expression (Neil et al. 2009;Xu et al. 2009). These findings are similar to our results (Figs. 4, 5). A difference, however, is that budding yeast CUTs are greatly stabilized by loss of Trf4, a key component of the exosometargeting TRAMP complex (Wlotzka et al. 2011;Frenk et al. 2014), while fission yeast CUTs are only marginally affected by loss of the TRAMP subunit Cid14 (Figs. 1, 2), as shown before (Zhou et al. 2015). This result indicates a TRAMP-independent mechanism of exosomal degradation of CUTs in fission yeast. Indeed, the poly(A)-binding protein Pab2, functioning in a complex called MTREC, targets meiotic or unspliced mRNAs and lncRNAs in fission yeast (McPheeters et al. 2009;St-André et al. 2010;Yamanaka et al. 2010;Lee et al. 2013;Egan et al. 2014;Zhou et al. 2015). Our data support a Pab2-dependent mechanism for targeting CUTs, showing a stronger derepression of CUTs in pab2 mutants than in cid14 mutants, although both mutants show much smaller effects than the exosome mutants rrp6 and dis3 (Figs. 1-3). Also in human cells, many lncRNAs are targeted for exosomal degradation by PABPN1, an ortholog of Pab2 (Beaulieu et al. 2012;Meola et al. 2016). Reminiscent of yeast CUTs, exosome depletion in mammals has revealed lncRNAs divergently transcribed from promoter regions of protein-coding genes Grzechnik et al. 2014;Quinn and Chang 2016). Whether the evolutionary conservation of this principle reflects functional importance of CUTs or their transcription, or simply that they are nonfunctional by-products of the basic mechanics of transcription, remains an open question.
The ∼850 SUTs in budding yeast are detectable in proliferating wild-type cells, and are processed differently from CUTs (Neil et al. 2009;Tuck and Tollervey 2013). SUTs could be considered analogous to the previously annotated lncRNAs in fission yeast, which can be readily detected in proliferating wild-type cells (Wilhelm et al. 2008;Rhind et al. 2011) and show less variable expression than novel lncRNAs across the different conditions (Fig. 3A). CUTs and SUTs almost exclusively originate from nucleosomedepleted regions at the ends of coding genes ). Our nucleosome-profiling data also suggest a strong tendency of Bidirectional and Intergenic lncRNAs to initiate in nucleosome-depleted regions upstream of positioned nucleosomes (Fig. 7). These nucleosome data are only approximate, as we did not determine the profiles under the different genetic or physiological conditions when lncRNAs are more highly expressed.
The ∼1700 budding yeast XUTs are derepressed upon deletion of the cytoplasmic exonuclease Xrn1 (ortholog of S. pombe Exo2), and are mostly antisense to mRNAs and anti-correlate with sense expression (van Dijk et al. 2011). These findings resemble our results (Figs. 4,5), and related recent findings in S. pombe (Wery et al. 2018). The targeting of XUTs by a cytoplasmic exonuclease implies their efficient export to the cytoplasm. However, the proposed inhibitory functions of XUTs on coding transcription are likely mediated cotranscriptionally, and so the relevance of their cytoplasmic export is unclear (Hansen et al. 2013;Tuck and Tollervey 2013). In budding yeast, XUTs are targeted by the NMD pathway before being degraded by Xrn1, and this pathway can be regarded as the last filter to dampen lncRNA expression (Malabat et al. 2015;Wery et al. 2016). Consistent with the cytoplasmic exonuclease acting as a backup system, we find that many CUTs and DUTs are also targeted by Exo2 and the NMD factor Upf1 (Figs. 2, 4C). Our XUTs and Antisense lncRNAs that engage with ribosomes, but not other classes of lncRNAs, are significantly enriched for RNAs derepressed in upf1Δ cells, supporting a role of the NMD in Exo2 degradation. Although there are overlaps of lncRNA expression in the absence of Upf1 and Exo2, much fewer lncRNAs are derepressed in upf1 than in exo2 mutants (Figs. 2, 4C). The absence of two other regulators of cytoplasmic RNA degradation, Ski7 and Pan2 (Lemay et al. 2010;Wolf and Passmore 2014), has only subtle effects on lncRNA derepression (Figs. 2, 4C). These results suggest that Exo2 plays the major role and can degrade lncRNAs independently of these other factors. However, a limitation of the current study is that RNA-seq only measures steadystate RNA levels, which integrates transcription and degradation. Findings from budding yeast indicate that mRNA levels can be adjusted by buffering mechanisms, allowing compensation of increased degradation by increased transcription or vice versa (Haimovich et al. 2013;Sun et al. 2017). Xrn1, the Exo2 ortholog, is required for this buffering. So it is possible that the weak derepression phenotypes of cytoplasmic RNA degradation mutants, other than exo2, reflect that lncRNA levels are efficiently buffered in these mutants. More work is required to investigate whether the buffering system is conserved in fission yeast and whether lncRNAs are subject to it.
Budding yeast has no analogous lncRNA class to the DUTs defined here, because the RNAi pathway is missing (Harrison et al. 2009). Our results show that RNAi is important to control lncRNA expression in fission yeast, most notably Antisense lncRNAs that are derepressed in late meiosis (Figs. 5,6). In fission yeast, RNAi can dampen RNA expression via either transcriptional or post-transcriptional mechanisms (Castel and Martienssen 2013;Smialowska et al. 2014). Our data do not allow to distinguish between these two possibilities. Although RNAi is not required for antisensemediated transcriptional repression at three meiotic mRNAs ), our and other results (Bitton et al. 2011) indicate a prominent global role of RNAi to suppress many Antisense lncRNAs. About 75% of all DUTs are Antisense lncRNAs. RNAi plays an even more important role than Exo2 in repressing Antisense lncRNAs, but also targets Bidirectional and Intergenic lncRNAs (Fig. 5A,B). It is not clear whether the RNAi machinery is involved to a similar extent in controlling lncRNAs in multicellular organisms.
NUTs are another class of unstable lncRNAs in budding yeast, which substantially overlap with CUTs and XUTs (Schulz et al. 2013). NUTs are detected upon depletion of Nrd1, a member of the Nrd1-Nab3-Sen1 (NNS) complex that promotes transcriptional termination of lncRNAs (Schulz et al. 2013). NNS also plays an important role in recruiting TRAMP for degradation of its targets (Tudek et al. 2014). No NNS complex was identified in fission yeast, and depletion of Seb1 impairs poly(A)-site selection but not RNA abundance (Lemay et al. 2016;Wittmann et al. 2017). Thus, a class corresponding to NUTs does not appear to exist in fission yeast. This conclusion is also consistent with Pab2, rather than TRAMP, being more important for exosomemediated degradation of lncRNAs.
The different lncRNA classes based on RNA-processing pathways, while useful, are fairly arbitrary and overlapping. The lncRNAs are targeted by multiple redundant or coordinating pathways in an intricate backup system, although one pathway is often dominant for a given lncRNA. This property has also been documented in budding yeast (Marquardt et al. 2011). Moreover, RNA-processing mutants can lead to cellular re-routing of RNA degradation. Accordingly, there are substantial overlaps between different lncRNA classes in both budding and fission yeast. We find that cells require either the nuclear exosome or cytoplasmic exonuclease to survive, with the absence of both pathways being lethal. While these pathways function in other aspects of RNA metabolism, this synthetic lethality points to the importance of dampening the extensive lncRNA expression, which are much more affected in the corresponding mutants than are mRNAs (Fig. 2). It is possible that the cytoplasmic exonuclease can serve as a backup to degrade transcripts that escaped degradation by the nuclear exosome. On the other hand, cells survive without the cytoplasmic exonuclease and RNAi or without the nuclear exosome and RNAi. Surprisingly, cells lacking both the cytoplasmic exonuclease and RNAi show fewer derepressed XUTs and DUTs than cells lacking only one of these pathways (Fig. 2). This suppression might reflect that lncRNAs that cannot be degraded by RNAi are effectively targeted by the nuclear exosome. Consistent with this possibility, absence of both the nuclear exosome and RNAi leads to poor growth and large numbers of derepressed lncRNAs ( Fig.  2; Supplemental Fig. S3). These findings indicate partially redundant roles for the nuclear exosome and RNAi pathways, which can back each other up with respect to many RNA targets. These two nuclear pathways can also degrade most XUTs that are further targeted by the cytoplasmic exonuclease. The RNAi and exosome pathways in fission yeast have overlapping functions to repress aberrant transcripts (Bühler et al. 2008;Zofall et al. 2009;Zhang et al. 2011) as well as meiotic mRNAs and other genomic regions (Yamanaka et al. 2010Sugiyama and Sugioka-Sugiyama 2011). Our study highlights that the nuclear exosome and RNAi pathways also cooperate to suppress thousands of lncRNAs.
The expression levels of most lncRNAs are highly induced in nondividing states (stationary phase and quiescence) and during meiotic differentiation, most notably in late meiosis when over 3000 Antisense lncRNAs are induced (Figs. 1, 6). These results raise the possibility that lncRNAs function during these conditions. It is known that unstable lncRNAs, normally targeted for rapid degradation, can become stabilized and functional under specialized conditions (Camblong et al. 2007;Houseley et al. 2008). Environmentally regulated changes to RNA quality-control activities can alter transcriptomes and mediate stress responses (Joh et al. 2016). RNAprocessing pathways might become down-regulated under certain physiological conditions, allowing lncRNAs to accumulate. The mRNA levels of relevant RNA-processing genes do not strongly change in response to our physiological conditions, although mRNAs encoding nuclear-exosome components decrease ∼2.7-fold during meiosis (Bitton et al. 2015a). Many meiotic mRNAs are repressed in mitotic cells by the RNAi and exosome pathways and derepressed during meiosis (Yamanaka et al. 2010Sugiyama and Sugioka-Sugiyama 2011). Derepression of lncRNAs during meiosis and other specialized conditions could involve similar regulation. Indeed, our findings indicate that Exo2 also plays an important role in repressing many lncRNAs, but also many middle meiotic genes (Mata et al. 2002) that are derepressed in exo2 mutants and meiosis. These results put Exo2 on the map as an important new regulator of meiotic gene expression.
In addition to derepression, the induction of lncRNAs could involve increased transcription (Castelnuovo et al. 2014). RNA-processing factors likely regulate RNA levels via coordinated interplays between transcription and degradation (Haimovich et al. 2013;Sun et al. 2017), and changes in this coordination could lead to the accumulation of different lncRNAs in different physiological conditions. Our data cannot distinguish between lncRNA regulation at the level of transcription or RNA decay.
Antisense transcripts are the most widespread class of lncRNAs. It has recently been reported that differences in poly(A) sites between convergent genes generate distinct antisense landscapes in budding and fission yeast (Liu et al. 2017). In our data, over 70% of coding sequences produce at least one Antisense lncRNA from the other strand. This finding complements and extends previous analyses on antisense transcription in fission yeast (Dutrow et al. 2008;Wilhelm et al. 2008;Zofall et al. 2009;Ni et al. 2010;Bitton et al. 2011;Rhind et al. 2011;Zhang et al. 2011;Chen et al. 2012;Marguerat et al. 2012;DeGennaro et al. 2013;Clément-Ziza et al. 2014;Eser et al. 2016;Wery et al. 2018). Antisense lncRNAs include CUTs, DUTs, XUTs and other lncRNAs, with XUTs and especially DUTs being strongly enriched (Fig. 5A,B). Thus, several RNA processing pathways can be involved in controlling Antisense lncRNAs. Previous studies have reported repressive effects of Antisense lncRNAs on their sense mRNAs (Ni et al. 2010;Bitton et al. 2011;Chen et al. 2012;Marguerat et al. 2012;Leong et al. 2014;Wery et al. 2018). Accordingly, we find a strong global tendency toward anticorrelation between Antisense lncRNA-mRNA expression levels under physiological conditions (Fig.  5C). In contrast, Antisense lncRNA expression shows a slight tendency toward positive correlation with mRNA expression under the genetic conditions (Fig. 5C). Thus, stabilization of antisense lncRNAs in the absence of different RNAprocessing factors appears generally not to be sufficient for mRNA repression. This finding suggests that Antisense lncRNAs often control mRNA expression at the level of transcription (e.g., by transcriptional interference or altered chromatin patterns) rather than functioning as transcripts. Alternatively, many Antisense lncRNAs might simply reflect opportunistic transcription, enabled by down-regulation of the dominant, overlapping mRNAs, with the anticorrelated expression reflecting passive, indirect effects. The ∼29% of protein-coding regions not associated with Antisense lncRNAs are enriched for highly expressed genes, suggesting that these genes are either protected from, or interfere with, antisense transcription. Despite the global anticorrelation (Fig. 5C), large numbers of Antisense lncRNAs go against this trend, indicating that the expression relationships between lncRNAs and mRNAs involve multiple processes and cannot be explained by a few regulatory or indirect mechanisms. This conclusion is consistent with the diverse findings on antisense lncRNA processes in other organisms (Pelechano and Steinmetz 2013;Mellor et al. 2016).

Conclusions
This study increases the number of lncRNAs annotated in fission yeast by almost fivefold, revealing both similarities and striking differences to lncRNA characteristics and regulation budding yeast. The novel lncRNAs are typically very lowly expressed but become derepressed in response to different genetic and physiological perturbations. In stark contrast, the mRNAs and annotated lncRNAs show less widespread changes in expression, especially in the genetic perturbations. The nuclear exosome, RNAi machinery, and cytoplasmic exonuclease are the dominant RNA-processing pathways degrading lncRNAs, used to define the CUTs, DUTs, and XUTs, respectively. Bidirectional lncRNAs are enriched for CUTs and translating ribosomes, and positively correlate with divergent mRNA expression. Antisense lncRNAs are enriched for DUTs and XUTs, are mostly derepressed in late meiosis, and negatively correlate with sense mRNA expression in physiological, but not in genetic conditions. Intergenic lncRNAs are enriched for lncRNAs that are not classified as CUTs, DUTs, or XUTs. The transcripton of Intergenic and Bidirectional lncRNAs initiates from regions that in wild-type cells are nucleosome-depleted, just upstream of a positioned nucleosome. Given their low expression and other features, it seems likely that any regulatory functions mediated by most lncRNAs are in cis and cotranscriptional.
Our findings highlight a substantial role of RNAi, in coordination with the nuclear exosome, in controlling a large number of lncRNAs typified by the new class of DUTs. Moreover, the findings reveal a prominent new function of the Exo2 cytoplasmic exonuclease, together with RNAi, in dampening the expression of both lncRNAs and mRNAs that become derepressed during meiosis. The nuclear exosome and cytoplasmic exonuclease together play an essential role for cell viability. The three RNA-processing pathways show overlapping roles and can target most lncRNAs with different affinities, forming an intricate, intertwined RNAsurveillance network. Besides these biological insights, this study provides broad data on diverse lncRNA characteristics and a rich resource for future studies on lncRNA functions in fission yeast and other organisms.

S. pombe strains
All strains, physiological and growth conditions (Edinburgh minimal media, EMM2, or Yeast Extract media, YE), and biological repeats used for RNA-seq in the current study are detailed in Supplemental Tables S1 and S2. Biological repeats are based on independent cell cultures of the different strains. The number of different samples was decided to broadly interrogate key genetic and physiological conditions, balancing biological insight with costs. The PCR-based approach (Bähler et al. 1998) was used for gene deletions of exo2, pan2, and strains used to generate double mutants. Double mutants among dcr1, rrp6, and exo2 were created by crossing the corresponding single mutants (Supplemental Table S1). Strain h − dcr1::nat ura4 − was generated using the PCR-based approach (Bähler et al. 1998). Random spore analysis was used to create the other single mutant strains with the correct mating-types. Strains were crossed and incubated on malt extract agar (MEA) for 2-3 d at 25°C. Tetrads were treated with zymolyase (0.5 mg/mL, MP Biomedicals Europe) and incubated at 37°C for at least 4 h to release spores. Spores were germinated on YE agar plates before being replica plated to selective EMM2 plates as appropriate. All deletion junctions were PCR verified (Bähler et al. 1998). Crosses and selection by random spore analysis were as follows: h + ade6-M216 leu1-32 ura4-D18 his3-D1 rrp6::ura4 was crossed with h − ura4-D18 with selection on EMM2 plates to create h + rrp6::ura4 ura4-D18; h − exo2::kanMX6 ade6-216 was crossed with h + ura4-D18 with selection on YE + kanamycin plates, and on EMM2 plates with or without uracil, to select for h + exo2:: kanMX6 ura-D18 and h − exo2:: kanMX6 ura-D18. Tetrad analysis was used to analyze the meiotic products resulting from the crosses in all combinations of the rrp6Δ, dcr1Δ, and exo2Δ single mutant strains. Strains were crossed and incubated on MEA plates for 2-3 d at 25°C. The resulting tetrads were dissected using a micromanipulator (Singer Instruments), and spores germinated on YE plates after 5 d of growth. Haploid colonies arising from germinated spores were then streaked to selective plates to test for KAN, NAT, and URA markers. All deletion junctions in double-resistant colonies were PCR-verified (Bähler et al. 1998).

Growth conditions
All mutant cell cultures were harvested at mid-log phase (optical density, OD 595 = 0.5). For stationary-phase experiments, wild-type cells were grown in EMM2 at 32°C. A sample representing 100% survival was harvested when cultures reached a stable maximal density. Colony forming units (CFUs) were measured every 24 h after this initial time-point, and another sample harvested when cultures reached 50% survival. For quiescence experiments, cells were grown in EMM2 at 32°C an OD 600 of 0.2, before being centrifuged, washed twice in EMM2 without nitrogen (NH 4 Cl), and cultured in EMM2 without nitrogen at 32°C. Cells under nitrogen starvation reached an OD 600 of 0.8 within 24 h, and were harvested at 24 h and 7 d after nitrogen removal. For meiotic timecourses, pat1-114 diploid cells were grown to mid-log phase before being shifted to EMM2 without nitrogen. Cells were incubated at 25°C overnight to synchronize them in G1 phase. Meiosis was induced by addition of NH 4 Cl to a final concentration of 0.5 g/L and incubation at 34°C (0 h time point). Cells were harvested by centrifugation of 50 mL cultures at 2300 rpm for 3 min, and pellets were snap-frozen and stored at −80°C prior to RNA extraction.

RNA-seq experiments and initial analyses
RNA was extracted from harvested cells using a hot-phenol method (Bähler and Wise 2017). The quality of total extracted RNA was assessed on a Bioanalyser instrument (Agilent). Strand-specific RNA-seq libraries were prepared using an early version of the Illumina TruSeq Small RNA Sample Prep Kit. For poly(A)-enriched samples, library preparation and sequencing protocols were as previously described (Lemieux et al. 2011), and for samples depleted for rRNAs (rrp6Δ, exo2Δ), as previously described (Bitton et al. 2014). RNA-seq libraries were sequenced on an Illumina HiSeq 2000 instrument, using single-end runs with 50 bp reads (The Berlin Institute for Medical Systems Biology, Germany). Reads were aligned to the fission yeast genome with the exonerate software (Slater and Birney 2005). The few reads which mapped equally well to multiple genomic locations were assigned at random to one of these locations. Reads containing up to five mismatches (not clustered at read ends) were kept for further analysis. For a previous study , we performed mapping quality tests to show that five mismatches work well for the S. pombe genome which has few repetitive regions. Between 20 and 50 million mappable reads were obtained for each library (∼80%-85% of total reads were mappable). Expression scores were calculated for annotated features using the genome annotation available in PomBase on 9th May 2011 (Wood et al. 2012), and "in house" Perl scripts as previously described . Reads per kilobase of transcript per million reads mapped (RPKMs) for annotated features correlated strongly between biological replicates (r Pearson > 0.98). Mapping and expression score pipelines were performed essentially as previously described (Lemieux et al. 2011). The reads were mapped only to the genome, and the reads were aligned in Exonerate using the default (ungapped) mode.

Segmentation of sequence data to define novel lncRNAs
Custom scripts for segmentation of RNA-seq data were written in R and Perl. A simple heuristic was designed to detect novel lncRNAs from RNA-seq data. This segmentation heuristic was optimized for its ability to detect the 1557 annotated lncRNAs, and validated by visual inspection of RNA-seq data. The following RNA-seq data from initial sequencing runs were pooled (two biological repeats each): rrp6-ts, exo2Δ, dis3-54, pab2Δ, nmt1-mtr4 (Lemay et al. 2014;Bitton et al. 2015a), ago1Δ, rdp1Δ, dcr1Δ, pan2Δ, upf1Δ, Stat 100%, Stat 50%, Quies 24 h, Quies 7d, Meiotic pool (Schlackow et al. 2013), YE1, and Reference (control) (Supplemental Tables S1 and S2). Segments were delimited from the pooled data using a 10 hits/bp cutoff. This approach improves the signal-to-noise ratio to increase sensitivity at the cost of specificity. However, nearly 100% of the novel lncRNAs were retained also when using only single samples, without pooling, applying the same threshold of >10 reads/bp for at least one sample. Thus, the final count of lncRNAs was not inflated by pooling of multiple conditions in which lncRNAs showed ≤10 reads. The advantage of having multiple diverse conditions was to uncover lncRNAs expressed only under specific conditions, and were therefore missed in previous analyses. Segments <100 bp apart and differing in pooled read density (average hits/bp of segment) by <10-fold were joined together. The ability to detect the annotated 1557 lncRNAs was judged on the percentage coverage of each segment overlapping a lncRNA, and the percentage coverage of each lncRNA overlapping a segment. To optimize these two values, we varied the following parameters: (i) hits/bp cutoff, (ii) fusing distance, (iii) imposing rule on whether to fuse based on fold difference in expression of consecutive segments, and (iv) varying this fold difference in expression at which consecutive segments were fused. Imposition of the criterion that consecutive segments should only be fused if their fold difference in expression meets a certain threshold aimed to detect lncRNAs near, but distinct from, mRNAs, while discarding data which likely represent misannotated untranslated regions. With the optimized segmentation procedure, annotated lncRNAs were covered at 92% by the detected segments.
Using the PomBase genome annotation (May 2011), segments overlapping annotations on the same strand, including untranslated regions, were removed. We discarded segments of <200 bp as lncRNAs are defined by an arbitrary minimal length cutoff of 200 nt (Mattick and Rinn 2015), reflecting RNA-seq library protocols that exclude small RNAs. The remaining consecutive segments >200 bp defined 5775 novel lncRNAs. Overall, 214 of the 487 novel lncRNAs reported by Eser et al. (2016) could not be validated using our segmentation (Supplemental Table S4). This analysis revealed that TSS uniquely called by Eser et al. (2016) show no signal in our RNA-seq data in any of the conditions, yet they often correspond to very strong signal from the opposite strand. These patterns raise the possibility that many of the 214 lncRNAs exclusively called by Eser et al. (2016) might be artifacts based on "leak-through" of the opposite-strand signal, which can result from some sequence library protocols (Perocchi et al. 2007).
To further verify the robustness of our segmentation, we applied the RNA-seq segmentation algorithm published by Eser et al. (2016) to annotate the TSS for all RNAs (coding and noncoding) in each of three of our key data sets: dcr1, exo2, rrp6 mutants. We used a detection threshold calculated by a bimodal distribution model, min-length of 200 and max-gap equal 80. The detection thresholds (reads per bp) with our dcr1, exo2, and rrp6 mutant data were 11.1, 7.6, and 9.7, respectively. We checked which of these TSS called by Eser's method overlap with previously annotated or novel ncRNAs. This analysis revealed that 4252 of 9100 called TSS in dcr1 mutants overlap with 3332 RNAs (2587 annotated, 745 novel), 4991 of 9810 called TSS in rrp6 mutants overlap with 3939 RNAs (3004 annotated, 935 novel), and 4426 of 9377 called TSS in exo2 mutants overlap with 3519 RNAs (2671 annotated, 848 novel). The data for these overlaps are provided in the last three columns of Supplemental Table S3. We then asked whether Eser's method detected much fewer novel ncRNAs than mRNAs, which would signify differences in the annotation quality between the two categories. However, Eser's method finds similar proportions of mRNAs and novel ncRNAs in the three mutant data sets (dcr1 mutant: 29.2% mRNAs, 22.6% ncRNAs; exo2 mutant: 30.7% mRNAs, 24.0% ncRNAs; rrp6 mutant: 36.2% mRNAs, 25.4% ncRNAs). This analysis shows that the novel ncRNAs are nearly as likely to be called by Eser's method as are the mRNAs, although the latter are much more highly expressed. We therefore conclude that novel ncRNA annotations can be detected almost as effectively by Eser's method as mRNAs.
We also compared the XUTs reported by Wery et al. (2018) to our annotations. Wery et al. define 1628 XUTs mapping to Chromosomes I-III, of which 1228 overlap with 1150 lncRNAs (575 annotated, 575 novel), and of which we classified 299, 315 and 84 as XUTs, CUTs, and DUTs, respectively. These differences reflect the different samples used to define the XUTs. While the XUTs annotated by Wery et al. (2018) are induced in the exo2 mutant, many are even more induced in rrp6 or dcr1 mutants, and we therefore defined them as CUTs or DUTs, respectively. The overlaps of the XUTs defined by Wery et al. (2018) with annotated or novel lncRNAs are indicated in Supplemental Table S3.
We defined three confidence classes based on RPKM values from all sequenced samples as follows: High confidence lncRNAs had ≥10 RPKM in at least one sample; medium confidence lncRNAs had <10 RPKM but ≥1 RPKM in at least one sample; and low confidence lncRNAs had <1 RPKM in all samples. We have set up a web tool to view the RNA-seq data for all lncRNAs and mRNAs in the different conditions (http://bahlerlab.info/ncViewer).

Analyses of RNA expression
Novel lncRNAs defined by the segmentation process described above, together with all annotated transcripts, were analyzed using the Bioconductor DESeq2 package (Love et al. 2014). Differentially expressed genes were defined as those being >2-fold induced (average of two biological repeats) or repressed and showing significant changes (adjusted P < 0.05) compared to three reference samples as determined by DESeq2. For hierarchical clustering of expression data, log 2 ratios were clustered in R with the pheatmap package, using the Euclidian distance measure and the ward.D or ward.D2 clustering options. For expression correlation analyses, we evaluated the similarity of expression levels of Bidirectional and Antisense lncRNAs to the expression levels of their neighboring mRNAs using normalized expression values across the entire data set. Vectors of mRNA-lncRNA pairs were generated and Pearson's correlation coefficients were computed. For lncRNAs associated with multiple mRNAs, only the nearest mRNA was considered. Functional enrichments of gene lists were performed using the AnGeLi tool which applies a two-tailed Fisher's exact test (Bitton et al. 2015b).
To test for nascent transcription of lncRNA genes, we analyzed recently published data from native elongating transcript sequencing (NET-seq) (Shetty et al. 2017). Normalized NET-seq data of proliferating cells (with or without depleted of Spt5) in "WIG" format were downloaded from GEO (accession: GSM2258030). The NET-seq targets the 3 ′ -end of nascent transcripts, and we systematically computed the normalized NET-seq signal across the entire lengths of the novel lncRNAs. The lncRNAs were considered present when signal was >0 in at least one replicate in any condition. The low threshold was required owing to the very low read numbers in the NET-seq data. This analysis suggested that 87.4% of all lncRNAs are transcribed in proliferating cells. Given the limited sensitivity of the NET-seq data, we expect this number to be a lower estimate.

Classification into CUTs, DUTs, and XUTs
Differential expression data from dcr1Δ, exo2Δ or rrp6Δ mutants were filtered to retain only transcripts that were significantly induced in ≥1 mutant compared to wild-type controls (expression ratio >2 and adjusted P-value <0.05). RPKM values from independent biological repeats for the differentially expressed RNAs were then standardized to have a mean value of 0 and a standard deviation of 1, followed by clustering using the Mfuzz clustering function in R (Kumar and Futschik 2007). The number of clusters "c" was set to 3 and the fuzzification parameter "m" to 1.25. To further reduce ambiguity when associating RNAs to clusters, the minimum membership value of a lncRNA belonging to a specific cluster was set to 0.7 (Kumar and Futschik 2007). For this classification, more recent PomBase annotations were used which contained only 1533 annotated lncRNAs (7308 annotated and novel lncRNAs in total).

Classification into bidirectional, antisense, and intergenic lncRNAs
To assess whether a given annotated or novel lncRNA overlaps with any mRNAs in either orientation, we systematically aligned the lncRNA coordinates relative to the annotation in Ensembl S. pombe, Assembly ASM294v2, release 33 (Flicek et al. 2014), enhanced by a modified annotation set that better delineates transcript boundaries. To this end, we exploited Transcription Start Sites (TSS) determined using Cap Analysis of Gene Expression (CAGE) (Li et al. 2015) and Transcription Termination Sites (TTS) defined using genome-wide polyadenylation site mapping (Mata 2013). For genes without these higher quality boundaries, we used the annotated TSS and TTS. All TSS and TTS used are provided in Supplemental Table S7. We called overlaps in either orientation when ≥1 nt was shared between transcripts. Using the same criteria, we tested for overlaps with novel lncRNAs that have recently been reported (Supplemental Table S4; Eser et al. 2016).
Using these overlap criteria, we classified all known and novel lncRNAs based on their proximity to nearby mRNAs. We defined lncRNAs as Intergenic if they do not overlap with any nearby mRNA, Sense-overlapping lncRNAs if they overlapped with any mRNA on the same strand, Antisense if they overlap ≥1 nt with a mRNA on the opposite strand, and Bidirectional if their TSS was <300 nt up-or downstream of a TSS of a mRNA on the opposite strand. Naturally, given the compact fission yeast genome, there was some overlap between these classes. We reassigned lncRNAs present in two classes using the following criteria. The lncRNAs classified as both Bidirectional and Intergenic (482 lncRNAs) or as Bidirectional and Antisense (1068 lncRNAs) were assigned to Bidirectional lncRNAs only. The 135 lncRNAs classified as both Sense-overlapping and either 86 Antisense or 49 Bidirectional lncRNAs were, respectively, assigned to Antisense or Bidirectional lncRNAs only.

Translation analysis of lncRNAs
For the ribosome profiling analysis, we systematically looked for overlaps between the translated regions defined before  and all annotated and novel lncRNAs. The significance of enrichments among different lncRNA classes was determined using the prop.test function in R.

Growth phenotypes of mutant cells
For a semi-quantitative analysis of cell growth, we used spot assays. After overnight preculture, yeast cells were adjusted to the same OD value (∼0.6). For each strain, 5 µL of six serial (fivefold) dilutions were spotted onto YE plates and grown at 32°C.

Nucleosome profiling
Mononucleosomal DNA (MNase digested) from exponentially growing wild-type (972 h − ) cells in EMM2 was generated as reported (Lantermann et al. 2009). Two independent biological repeats were performed. Sequencing libraries from MNase-digested DNA were prepared using the NEBNext ChIP-Seq Library Prep Master Mix Set for Illumina (E6240S). Pair-end 50 bp reads were obtained with an Illumina MiSeq sequencer at the Genomics and Genome Engineering Facility at the UCL Cancer Institute. MNAse sequencing data was mapped using Bowtie2 (Langmead and Salzberg 2012). Nucleosome maps for visualization were performed with nucwave (Quintales et al. 2015), following the web recommendations (http:// nucleosome.usal.es/nucwave/). Data were analyzed using the "computeMatrix reference-point" and "heatmapper" functions from the deeptools package with transcription start sites as reference points (Ramírez et al. 2014).

SUPPLEMENTAL MATERIAL
Supplemental material is available for this article.