extfieldsofstudy: Medicine Biology

text:
Loss of SMYD1 Results in Perinatal Lethality via Selective Defects within Myotonic Muscle Descendants

SET and MYND Domain 1 (SMYD1) is a cardiac and skeletal muscle-specific, histone methyl transferase that is critical for both embryonic and adult heart development and function in both mice and men. We report here that skeletal muscle-specific, myogenin (myoG)-Cre-mediated conditional knockout (CKO) of Smyd1 results in perinatal death. As early as embryonic day 12.5, Smyd1 CKOs exhibit multiple skeletal muscle defects in proliferation, morphology, and gene expression. However, all myotonic descendants are not afflicted equally. Trunk muscles are virtually ablated with excessive accumulation of brown adipose tissue (BAT), forelimb muscles are disorganized and improperly differentiated, but other muscles, such as the masseter, are normal. While expression of major myogenic regulators went unscathed, adaptive and innate immune transcription factors critical for BAT development/physiology were downregulated. Whereas classical mitochondrial BAT accumulation went unscathed following loss of SMYD1, key transcription factors, including PRDM16, UCP-1, and CIDE-a that control skeletal muscle vs. adipose fate, were downregulated. Finally, in rare adults that survive perinatal lethality, SMYD1 controls specification of some, but not all, skeletal muscle fiber-types.

Introduction
Skeletal myogenesis is a critical biological process that is essential for development, regeneration, and homeostasis [1,2]. A number of key transcription factors (TFs) control skeletal muscle development, including regulators of progenitors (e.g., PAX3 and PAX7), regulators of myogenic commitment (e.g., MYF5, MYOD1, and MYF6), and regulators of embryonic myocyte fusion (e.g., MYOG). These and other TFs act as a complex and often compensatory overlapping network to regulate both spatial and temporal expression.
SET and MYND Domain 1 (SMYD1) is a cardiac and skeletal muscle-specific, chromatin-binding factor that regulates growth and morphogenesis during fetal and neonatal development in mice [3,4], humans [5,6], and zebrafish [7,8]. SMYD1 is expressed within the heart and within fast-twitch skeletal muscle cells [3][4][5][6]. SMYD1 family members catalyze methylation of histone H3 lysine K4 in vitro, and as initially suggested by Gottlieb et al. [3], act within the developing heart as transcriptional repressors in vivo [4][5][6]. In some cases [4], SMYD1 inhibits the transcription of target genes by binding directly to target promoters. Alternatively, SMYD1 can serve as a direct target of TFs, including serum response factor (SRF) and MYOG, to direct sarcomere function in vivo [9]. Within sarcomeres, SMYD1 localizes to the M-line, where it physically associates with myosin [10], and within the nucleus, where

Mating Scheme
MyoG-Cre mice were provided by Dr. Eric Olson as F1 first generation hybrids of C3H and C57BL/6J. Recombination induced by MyoG-Cre is initiated in anterior somites at E9.0, and by E9.5, expression is detected in all somites and is confined and maintained in muscle tissue throughout embryogenesis and into adulthood [18]. The same breeding scheme as described previously [3] was employed to create Smyd1 Flox/KO ; MyoG-Cre +/− mice (BL6.FVB.129.C3H hybrids) as well as the following breeders of the same background: Smyd1 Flox/WT ; MyoG-Cre +/− and Smyd1 KO/WT ; MyoG-Cre +/− . However, in order to avoid breeding Cre to homozygosity, we began crossing the above genotypes with Cre-negative mice harboring the Smyd1 Flox and Smyd1 KO alleles. It later became apparent that Smyd1 Flox/Flox ; MyoG-Cre +/− and Smyd1 Flox/KO ; MyoG-Cre +/− displayed the same phenotype. In order to simplify the number of genotypes created and to create higher numbers of interesting genotypes, we mated Smyd1 Flox/Flox to Smyd1 Flox/WT ; MyoG-Cre +/− . These mice were only crossed to one another, thereby maintaining an inbred, but non-isogenic, hybrid population that was similar, but isolated, from the other Cre-induced CKO lines.

In Situ Hybridization
Embryos at 15.5 days post coitum (dpc) or younger were extracted in DEPC PBS (1 mL DEPC/L, autoclaved) and tail pieces were removed for genotyping. The remainder of the embryo was fixed in at 10-fold volumes of 4% paraformaldehyde at 4 • C on a rotator overnight. Samples were then washed and either stored at 4 • C or embedded in paraffin for tissue sectioning as previously described [19].

Protein/Immunohistochemistry
Embryos and adult tissues were extracted in PBS and fixed in at least 10-fold volumes of 10% formalin or 4% paraformaldehyde at room temperature on a rotator for 24-72 h. Embryos older than 15.5 dpc were skinned before fixation. Tail segments prior to fixation were genotyped. Samples were washed three times for at least 1 h in PBS and either stored at 4 • C or prepared for paraffin embedding and tissue sectioning as previously described [19].

End-Point RT-PCR
Trizol homogenates were thawed and RNA was isolated according to the Invitrogen Trizol protocol. For embryonic samples, isopropanol precipitation was supplemented with 5-10 µg of RNAse free glycogen. After drying the RNA pellet, samples were resuspended in 25 µL deionized DEPC water and immediately DNase1 treated with 15 µL mix D (2.7× DNase I buffer (Invitrogen, Carlsbad, CA, USA); 1 U/µL RNaseOut (Invitrogen, Carlsbad, CA, USA); 2 U/µL DNase I (Roche, Indianapolis, IN, USA)) to a final volume of 40 µL and incubated at 37 • C for 1 h. The enzyme was heat inactivated in 4 µL of 25 mM EDTA and incubation at 70 • C for 10 min. RNA quality and concentration were analyzed using a nanodrop spectrophotometer. Samples were boiled for 5 min and iced for 5 min before adding 13 µL of mix 2 (3.8× 1st Strand Buffer, 20 mM DTT, 1 U/µL RNaseOut (Invitrogen, Carlsbad, CA, USA)). After pre-warming the reaction mixture for 2 min at 37 • C, 2 µL of 200 U/µL MMLV (Invitrogen, Carlsbad, CA, USA) was added to each reaction except the negative control. Samples were reverse transcribed for 1.5 h at 37 • C. The enzyme was heat inactivated at 85 • C for 5 min and the samples were stored at 4 • C. Next, 1-2 µL of the cDNA reaction were used for RT-PCR with various gene primer sets ( Figure S1) or with primers for Smyd1: 5 -ctgcgggcagtgcaagtt-3 and 5 -gctcacaggagcagtcaaagtag-3 ; β-actin: 5 -tacgagggctatgctctc-3 and 5 -cgcagctcagtaacagtc-3 .

RT-qPCR
Total cellular RNA was isolated using an RNeasy Mini Kit (QIAGEN, Santa Clara, CA, USA). cDNA synthesis was performed with qScript cDNA supermix (Quanta Biosciences, Inc., Denver, CO, USA). RT-qPCRs were performed using PerfeCTa SYBR Green FastMix (Quanta) with 1 µL of 20×-diluted cDNA generated from 500 ng of total RNA. RT-qPCR primers were designed to amplify the junction between two exons. Primer sequences for qPCR are listed in Figure S2. Ct values were normalized against GAPDH. qPCR was performed with an initial denaturation step of 10 min at 95 • C, followed by 15 s at 95 • C and 60 s at 60 • C for 30 cycles using an MX3000P™ QPCR system (Stratagene, La Jolla, CA, USA). The normalized level of a mRNA was determined as 2−Ct(GOI)/2−Ct(CTL), where Ct is the threshold cycle, GOI is the transcript of interest, and CTL is the housekeeping control (assuming that Ct is inversely proportional to the initial concentration of mRNA and that the amount of product doubles with every cycle). The statistical significance of RT-qPCR results was determined by two-tailed t-test in Microsoft Excel. All sample values were normalized to GAPDH levels (primer sequences noted in Figure S2). GAPDH levels were not significantly altered among different groups (data not shown).

RNA In Situ Hybridization
Embryos were harvested at various pre-and postnatal gestation times, fixed in 4% paraformaldehyde overnight, and in situ hybridization was performed by whole-mount hybridization as described previously [3,19]. Following pre-hybridization procedures, tissue sections were hybridized at 55 • C with a sense and antisense riboprobe for murine Myogenin. Briefly, RNA antisense MyoD was obtained by RT-PCR after denaturing RNA for 10 min at 70 • C in PCR buffer (2.5 mM MgCl2), 1 mM dNTPs, 50 ng of oligo primer (dT) (Gibco, Carlsbad, CA, USA) and 200 units of the reverse transcriptase enzyme Superscript II (Gibco, Carlsbad, CA, USA). Primers used for the amplifications designed from a combination of the complementary MyoD gene sequence with the T7 bacteriophage promoter sequence in its 5 end were 5 -agcactacagtggcgactca-3 and 5 -ggccgctgtaatccatca-3 . Hybridization was conducted overnight with antisense and sense probes ( Figure S2) using~7.5 × 10 5 disintegrations per minute(dpm) per slide. The unhybridized probe was removed through stringent washes and treatment with RNase A. Slides were subsequently coated with K.5 nuclear emulsion (Ilford, UK) and exposed at 4 • C for 12 days. Slides were developed, counter-stained with hematoxylin, and examined using bright and dark field optics.

MyoG-Cre Deletion is Specific to Skeletal Muscle
As described in the Materials and Methods section, we analyzed Cre excision fidelity for genomic DNA recombination of the Smyd1 locus in heterozygous mice. As shown in Figure 1A, there was no detectable deletion of the Flox compared to the wildtype (WT) allele, but efficient recombination was observed specifically within skeletal muscle tissue. We expect that 20% of pups to be homozygous for Smyd1-deficient alleles, 50% of pups should be heterozygous, and 30% should be wild type (WT). Litters were genotyped at weaning and the observed numbers were recorded. There was no significant reduction in litter size between birth and weaning.

MyoG-Cre-Mediated Smyd1 CKO Results in Perinatal Lethality
Unexpectedly, of the first 350 pups weaned and genotyped at 3-6 weeks, there were zero Smyd1 conditional knockouts (CKOs; Figure 1B). Dead pups were found during and shortly after birth, but significant numbers of pups did not disappear between postnatal day 1 and weaning. The birth of a few litters was observed. CKO pups were difficult births, probably because of their increased size due to extreme edema. These pups were often crushed during their birthing process. Those that were born intact showed no signs of movement nor breathing.
Genotypic ratios for these litters were analyzed from E9.5 to E17.5, and no significant variance from the expected Mendelian ratios were observed (Table 1). Fetuses were extracted at E18.5-an age at which they could be fostered and survive if viable. However, unlike their littermates, MyoG-Creinduced CKOs at this age made no attempts at locomotion nor breathing. Thus, CKOs were not dying simply due to birthing problems, but never showed signs of life during late stage embryogenesis. Whether or not the heart was pumping at this age was not determined, but without locomotion, these embryos could not survive. We expect that 20% of pups to be homozygous for Smyd1-deficient alleles, 50% of pups should be heterozygous, and 30% should be wild type (WT). Litters were genotyped at weaning and the observed numbers were recorded. There was no significant reduction in litter size between birth and weaning.

MyoG-Cre-Mediated Smyd1 CKO Results in Perinatal Lethality
Unexpectedly, of the first 350 pups weaned and genotyped at 3-6 weeks, there were zero Smyd1 conditional knockouts (CKOs; Figure 1B). Dead pups were found during and shortly after birth, but significant numbers of pups did not disappear between postnatal day 1 and weaning. The birth of a few litters was observed. CKO pups were difficult births, probably because of their increased size due to extreme edema. These pups were often crushed during their birthing process. Those that were born intact showed no signs of movement nor breathing.
Genotypic ratios for these litters were analyzed from E9.5 to E17.5, and no significant variance from the expected Mendelian ratios were observed (Table 1). Fetuses were extracted at E18.5-an age at which they could be fostered and survive if viable. However, unlike their littermates, MyoG-Creinduced CKOs at this age made no attempts at locomotion nor breathing. Thus, CKOs were not dying simply due to birthing problems, but never showed signs of life during late stage embryogenesis. Whether or not the heart was pumping at this age was not determined, but without locomotion, these embryos could not survive. Therefore, the CKOs are perinatal lethal-a result unanticipated given that only a single previous example of MyoG-Cre-induced lethality had been observed [17]. We return to this point in the present context below.

Skeletal Muscle Development in the MyoG-Cre-Induced Smyd1 CKO is Significantly Derailed
Viable embryos initially were inspected visually from E11.5 through E18.5 for gross defects. At E12.5, embryos appear normal after sectioning (data not shown). Some embryos were outwardly abnormal as early as E13.5, but by E15.5, CKOs were easily discerned from littermates by eye.
In order to analyze potential skeletal muscle disruption in detail, CKO pups at various embryonic ages were sectioned and compared to their WT littermates. At E14.5, Haematoxylin and Eosin (H&E)-stained transverse sections of CKOs displayed a marked subcutaneous edema most visible in the cervical region (arrows in Figure 2A and data not shown).
Next, we performed in situ hybridization on transverse sections employing as probe, Myog, which marks all myogenic lineages, on sections prepared at E14.5 ( Figure 2B). While the sections are not matched precisely in the rostro-caudal plane, there is a clear reduction of Myog-marked myofiber density in the CKO. However, muscle groups were affected differently. For example, at E15.5, CKO limbs displayed less cellularity than controls, in correspondence with the Myog in situ hybridization ( Figure 2B). Unanticipatedly, the pectoralis (white arrow) and sternohyoid muscles (yellow arrow) were no longer visible in sections; i.e., missing entirely ( Figure 2B). We also noted that the dorsal muscles (solid blue arrows) were significantly diminished, while other muscles, including the limbs, were clearly less affected. There also appeared to be an excess of brown adipose tissue (blue and white arrows) in the CKO ( Figure 2B and re-addressed below).
Failure of CKO muscular development also can be seen in the H&E-stained cross-sections of Figure 2C. For example, the scapula region revealed that the surrounding muscle cells were less dense and less organized (Figure 2(Ca-d)). Also, there was an emergence of quite large, multinucleated cells-myoblasts that had fused, but apparently had failed to elongate into myotubes (marked with arrows in Figure 2(Cd)). CKO tongues displayed the same phenotypes as scapula muscles, with evident multinucleated, unfused structures, particularly evident in the CKO tongue ( Figure 2(Ch), arrows). However, the masseter muscles appeared totally unscathed from observable morphologic defects (Figure 2(Ci-l)). Myogenin in situ hybridization of WT and CKO transverse sections reveal reduction of myocytes and general affliction of some, but not all, muscles. The pectoral (white arrow) and sternohyoid (yellow arrow) muscles did not develop in the CKO; the dorsal muscles (solid blue arrow) are significantly diminished, but other muscles (e.g., limbs) are clearly less afflicted, if at all. There is an excess of tissue tentatively classified as brown adipose tissue (blue and white arrows) in the CKO. (C) H&E stains (ah) of WT and CKO sections at E15.5. Scapulae (a-d) show asynchronous elongation and emergence of large, multinucleated cells-myoblasts that did not fuse-but failed to elongate into myotubes (black arrows in d). Stains of CKO tongues (e-h) displayed the same phenotypes as above with evident multinucleated, unfused structures (black errors in h). CKO masseter muscles (i-l) had no observable morphologic defects.
Failure of CKO muscular development also can be seen in the H&E-stained cross-sections of Figure 2C. For example, the scapula region revealed that the surrounding muscle cells were less dense and less organized ( Figure 2C(a-d)). Also, there was an emergence of quite large, multinucleated cells-myoblasts that had fused, but apparently had failed to elongate into myotubes (marked with arrows in Figure 2Cd). CKO tongues displayed the same phenotypes as scapula muscles, with evident multinucleated, unfused structures, particularly evident in the CKO tongue (Figure 2Ch, arrows). However, the masseter muscles appeared totally unscathed from observable morphologic defects ( Figure 2C(i-l)).
In summary, SMYD1 skeletal muscle deficiency does not affect all myotonic descendants equally. Trunk muscles are nearly ablated, forelimb muscles are disorganized and improperly differentiated, but some muscles, such as the masseter, are normal. We observed that muscles appear normal as late as E12.5, following post-migration from the somitic myotome (data not shown). However, it is possible that the molecular effect occurred earlier than the visible effect, and Smyd1 is involved in some essential event(s) prior to this migration. Finally, TUNEL staining on sections of WTs and CKOs detected no difference in the number of apoptotic cells (data not shown).
These results along with the increased numbers of large multinucleated cells near areas of reduced cell density suggest that the paucity of myofibers is secondary to degeneration. However, we cannot distinguish whether the myofibers did not form or if they instead underwent rapid degeneration. In summary, SMYD1 skeletal muscle deficiency does not affect all myotonic descendants equally. Trunk muscles are nearly ablated, forelimb muscles are disorganized and improperly differentiated, but some muscles, such as the masseter, are normal. We observed that muscles appear normal as late as E12.5, following post-migration from the somitic myotome (data not shown). However, it is possible that the molecular effect occurred earlier than the visible effect, and Smyd1 is involved in some essential event(s) prior to this migration. Finally, TUNEL staining on sections of WTs and CKOs detected no difference in the number of apoptotic cells (data not shown).
These results along with the increased numbers of large multinucleated cells near areas of reduced cell density suggest that the paucity of myofibers is secondary to degeneration. However, we cannot distinguish whether the myofibers did not form or if they instead underwent rapid degeneration.

MyoG-Cre-Induced Smyd1 CKOs Exhibit Proliferative Defects and Non-Autonomous Deposits of Brown Adipose Tissue
To determine whether the apparent reduction in cellular density is due to a decrease in proliferation, sections were analyzed by immunohistochemistry for Ki67, a protein present in all active, but not in quiescent, stages of the cell cycle [20]. CKO skeletal muscles analyzed at E12.5 showed normal proliferation (data not shown). However, at E14.5, we observed reduction in Ki67 staining in dorsal sections ( Figure 3A,B) but not in muscles surrounding the CKO humorous ( Figure 3E,F). Notably, the leading edge of proliferation within the CKO tongue epithelium is significantly lagging and appears to be anatomically misshaped (blue arrow in Figure 3D). proliferation, sections were analyzed by immunohistochemistry for Ki67, a protein present in all active, but not in quiescent, stages of the cell cycle [20]. CKO skeletal muscles analyzed at E12.5 showed normal proliferation (data not shown). However, at E14.5, we observed reduction in Ki67 staining in dorsal sections ( Figure 3A,B) but not in muscles surrounding the CKO humorous ( Figure  3E,F). Notably, the leading edge of proliferation within the CKO tongue epithelium is significantly lagging and appears to be anatomically misshaped (blue arrow in Figure 3D). We suspected that the Ki67 dense regions (particularly evident in CKO dorsal muscle (red arrows)) were fat deposits. Adipocytes can be divided into classes of white (WAT), brown (BAT), or beige cells [21]. We tentatively assigned these adipocytes as "standard" BAT based on: (1) the eosin component of H&E stain distinguishes WAT from BAT (i.e., eosin colors BAT, but not WAT); (2) unfixed CKO cryostat sections were enriched in deposits staining positively with Sudan Black-a dye specific for BAT (data not shown). We suspected that the Ki67 dense regions (particularly evident in CKO dorsal muscle (red arrows)) were fat deposits. Adipocytes can be divided into classes of white (WAT), brown (BAT), or beige cells [21]. We tentatively assigned these adipocytes as "standard" BAT based on: (1) the eosin component of H&E stain distinguishes WAT from BAT (i.e., eosin colors BAT, but not WAT); (2) unfixed CKO cryostat sections were enriched in deposits staining positively with Sudan Black-a dye specific for BAT (data not shown).
In CKOs, BAT deposits were both increased as well as more heavily stained with H&E (pink arrows in Figure 3A,B). The CKO BAT area appears to be surrounded by other smaller cells-most likely lymphocytes and other immune cells. BAT composition, function, and numbers undergo significant changes in response to nutritional and mitochondrial demands (readdressed below).
Excess BAT was confirmed in other tissue sections (e.g., Figure 2B,C). Perhaps relevant in this context was our observation that, following shRNA knockdown of SMYD1 in C2C12 myocytes, higher levels of triglycerides-a component of adipose tissue-were apparent (L. Zhu and H.O. Tucker, data not shown). However, a previously considered essential BAT regulator, Prdm16, was downregulated-the opposite effect of that which was predicted (Figure 4).

Skeletal Muscle SMYD1 Deficiency Does not Affect Transcript Levels of Myogenic Regulatory Factors
In an attempt to identify pathways disrupted by skeletal muscle loss of Smyd1, we analyzed the transcript pool at E12.5 in an attempt to preclude secondary and indirect effects. Because the myogenic program was significantly derailed, we first examined 14 muscle-specific TFs, including the helix-loop-helix regulatory factors (MRFs) MYOG, MYOD1, MYF5, and MYF6 [22]. As shown both by real time RT-qPCR ( Figure 4) and by endpoint RT-PCR at E13.5 ( Figure S3), no significant changes in expression were observed.
likely lymphocytes and other immune cells. BAT composition, function, and numbers undergo significant changes in response to nutritional and mitochondrial demands (readdressed below).
Excess BAT was confirmed in other tissue sections (e.g., Figure 2B,C). Perhaps relevant in this context was our observation that, following shRNA knockdown of SMYD1 in C2C12 myocytes, higher levels of triglycerides-a component of adipose tissue-were apparent (L. Zhu and H.O. Tucker, data not shown). However, a previously considered essential BAT regulator, Prdm16, was downregulated-the opposite effect of that which was predicted (Figure 4).

Skeletal Muscle SMYD1 Deficiency Does not Affect Transcript Levels of Myogenic Regulatory Factors.
In an attempt to identify pathways disrupted by skeletal muscle loss of Smyd1, we analyzed the transcript pool at E12.5 in an attempt to preclude secondary and indirect effects. Because the myogenic program was significantly derailed, we first examined 14 muscle-specific TFs, including the helix-loop-helix regulatory factors (MRFs) MYOG, MYOD1, MYF5, and MYF6 [22]. As shown both by real time RT-qPCR ( Figure 4) and by endpoint RT-PCR at E13.5 ( Figure S3), no significant changes in expression were observed. Factors critical to BAT formation within the immune system (Immune BAT) and factors critical to progenitor differentiation to skeletal muscle or adipose tissue (Muscle vs. Adipose Fate) are significantly altered in the CKO relative to WT. E12.5 embryo preparation, RNA preparation, RT-qPCR, and associated analyses were performed as detailed in the Materials and Methods section employing primer-pairs listed in Figure S2. Data from three biological replicates were averaged. Absolute quantification of each target was performed using a standard curve as a reference (Roche Light Cycler software version 1.5, Pleasanton, CA, USA). PCR primers were designed using Primer3 software [23,24]. Significance was determined by student's t-test as: * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001.
Given the severity of the phenotype, this result was unanticipated. However, MRFs also were not targeted in the only other reported case of Myog-Cre-induced lethality [17]. Instead, the target of Myog-Cre was the Mef2c muscle enhancer, whose lethality occurred perinatally or postnatally (P) depending on genetics [17]. In a C57BL/6 mixed background, 100% lethality was observed at P10, whereas in other backgrounds (e.g., 129/SvEv)~50% of the CKO mice survived. However, the survivors displayed a fiber-type switching phenotype; i.e., myoblasts appeared normal, but differentiation from myoblasts to myotubes was blocked.
The MyoG-Cre:Smyd1 and Mef2c CKOs share a number of similarities. For example, we observed an abundance of primary myogenic cells but also the appearance of large multinucleated (fused) myoblasts that had failed to elongate (Figure 2(B,Cf,j)). This suggested that CKO myoblasts did not properly differentiate into myotubes-an event controlled to a great extent by MYOD1 and MYOG [22]. However recent reports also have linked similar highly polarized fusion deficiencies to melanoma cell adhesion molecule (MCAM) [22,25,26]. MCAM asymmetrically distributes at the tip of elongating myotubes where it colocalizes with actin binding and cell polarity regulators to facilitate directional elongation.
A possible disruption unrelated to MRFs is SMYD1 enzymatic function. SMYD1 is a musclespecific histone methyltransferase [5,27] that is capable of globally altering the transcriptome post-transcriptionally, while yielding no apparent alteration in transcript abundances. This is further evidenced by the lack of gene expression differences in Smyd1 WT and CKOs as determined by microarrays of both E12.5 and E14.5 muscle (T.L.R. and H.O.T; data not shown). Perhaps SMYD1 functions here by posttranslational methylation-an event that has been documented in a number of studies [5,27]. However, there is no evidence that any of the MRFs are methylated nor otherwise posttranslationally modified in a manner which would be catalyzed directly by SMYD1.

Central Regulators of BAT are Unaffected by SMYD1 CKO
A striking phenotypic feature of the MyoG-Cre:Smyd1 CKO was loss of regional BAT (Figures 2B  and 3). Transcriptional control of BAT formation has been analyzed extensively [28]. The TFs peroxisome proliferator-activated receptor gamma (PPARγ) [29,30] and the CCAAT/enhancer-binding proteins (C/EBPα/β/δ) [31] are essential transcriptional regulatory components that precede the formation of BAT. A downstream consequence of these pathways is de novo mitochondrial biogenesis activation [28].
As shown in Figure 4, RT-qPCR analyses detected no alteration of any of these factors in Smyd1-MyoG CKO muscle relative to controls. Thus, classical mitochondrial BAT biology appears unscathed by loss of SMYD1.

SMYD1 Regulates Immune Factors Critical for BAT Development and Physiology
Alternatively, regulation of BAT metabolism by both the innate and the adaptive arms of the immune system has been well established [32][33][34][35]. Recently, the impact and the molecular characterization of BAT-associated T regulatory (Treg) cells have been defined [36]. Several proinflammatory markers are induced by Treg-deficient BAT, including the cytokine interleukin 6 (IL-6), the chemokine ligand 2 (CCL2), and tumor necrosis factor alpha (TNFα) [37].
As anticipated from those reports, several previously demonstrated Treg-BAT factors were upregulated by SMYD1 loss (Figure 4). These include IL-6, CCL2, TNFα, and several additional pro-inflammatory chemokines (CCL6, 7, and 9). Also, as shown previously [36], we observed modest upregulation of Foxp3, the central transcriptional regulator of Treg differentiation, and several of its direct target genes, including BCL11, RUNX1, and CBFβ. Both BCL11 and RUNX1 also have been shown to reciprocally regulate Foxp3 [38].
The relevance of TNFα and IL-6 modulation with respect to muscle biology is particularly striking. TNFα is essential for incorporating fatty acids into muscle tissue, whereas there is a direct correlation between IL6 and fatty acid oxidation in skeletal muscle [39]. Both are inflammatory cytokines required for the maintenance and differentiation of skeletal muscle. TNFα plays a bipolar role in muscle regulation. In its absence, myogenic differentiation is blocked by failure to activate p38 [40]. In contrast the absence of TNFα strongly contributes to cachexia, or wasting atrophy [41]. Overexpression and knockdown experiments have documented that the myogenic markers Myog, α-actin, and p21 directly correlate with IL-6 levels [42]. Although TNFα levels are extremely low in SMYD1 CKO muscle (Figure 4), they are not strictly dependent on SMYD1 KD levels ( Figure S4). We also observed a tight negative correlation of IL-6 levels with gender at E13.5 ( Figure S4).

SMYD1 as a Potential Determinant of Muscle vs. Adipose Fate
Skeletal muscles (slow and fast twitch) and fat tissues (white, beige, and brown) differ broadly in morphology, anatomical locations, transcriptional profiles, and physiologic function [43,44]. BAT and skeletal muscle derive from the same multi-potent progenitors that originate within the central dermomyotome [45,46]. Expression of PAX7, an essential multi-potent progenitor TF, is turned off in dermal and BAT but continues to be expressed in MYOD+ embryonic and postnatal muscle progenitors [47,48].
It was shown recently that two MRFs, MyoD and Myf5, act downstream of PAX7 to repress the BAT determining factor, PRDM16, thereby acting as a "molecular switch" to repress BAT fate [48].
To determine if SMYD1 might be a component of this fate-switch circuit, we compared the levels of the above three components in E12.5 WT and CKO embryos. While we observed no effect on PAX7, MYOD nor MYF5, CKO skeletal muscle was depleted of PRDM16-a quintessential requirement for BAT fate choice ( Figure 4).
These data suggest that SMYD1 may control muscle vs. adipose fate via modulation of Prdm16 via a pathway different from that described for PAX7. Loss of PRDM16 in adipose tissue results in significant upregulation of thermogenic genes, including the brown fat-uncoupling protein, UCP-1 [49], which previously was touted as a "master regulator" of conventional BAT [28][29][30], and cell death activator-A (CIDE-A), which controls metabolic rates and higher lipolysis in BAT [50]. Yet as shown in Figure 4, both of these downstream modulators of PRMD16 are downregulated in MyoG-Cre/Smyd1 muscle. This result was the opposite to what would be expected if SMYD1 directly repressed PRDM16, as levels of these direct PRDM16 targets would be upregulated. An alternative mechanism is therefore required.
One possibility is that SMYD1 independently regulates transcription of PRDM16, UCP-1, and CIDE-A. Or, perhaps SMYD1 executes Prdm16 activation via transcriptional or posttranscriptional control. Finally, SMYD1 may mediate cross-talk with adipogenic energy stores or absorption of nutrients. Alternatively, if Smyd1-deficient multipotent progenitors are incapable of signaling the break-down of BAT to release metabolites, BAT would accumulate more efficiently in CKOs than in control embryos.
None of these hypotheses explain why the proliferation marker, Ki67, stained BAT more densely in SMYD1-deficient muscle (Figure 3). One explanation is that while PAX7 is repressed in normal BAT, it continues to be expressed in proliferating muscle progenitor cells in embryonic and young postnatal muscle that give rise to BAT [51,52]. Enlarged BAT foci also have been attributed to increased proliferation and differentiation of brown pre-adipocytes [53]. Alternatively, the observation might have a physiologic explanation. For example, SMYD1-deficient muscle might be unable to efficiently absorb nutrients. In this scenario, excess nutrients within the system could be adsorbed by brown adipose, resulting in abnormal proliferation and expansion of BAT.

SMYD1 Controls Specification of Muscle Fiber-Type in Rare Adult Survivors
After mating and genotyping over 190 litters of MyoG-Cre/Smyd1 mice, we observed no postnatal CKOs on the non-isogenic background described in the Materials and Methods section. However, two of four litters produced by one particular mother-son cross yielded pups in which CKOs survived and appeared outwardly normal into adulthood-albeit not at normal Mendelian ratios ( Figure 5A).
Cre expression was lower than in age-matched control embryos ( Figure 5B). Because our mouse lines were not isogenic, this particular branch may over time have lost or gained a genetic or epigenetic modification(s) that delayed/reduced Cre expression. If so, the data indicate that SMYD1 is not required for skeletal muscle function per se but is essential for maintaining the development and differentiation of skeletal muscle. This is supported by the Myf6-Cre data reported previously [12,13] and noted below.  DNA analyses revealed that a muscle-specific recombination event had occurred within the Smyd1 locus. The mutation(s) did affect SMYD1 expression at the RNA level, perhaps because MyoG-Cre expression was lower than in age-matched control embryos ( Figure 5B). Because our mouse lines were not isogenic, this particular branch may over time have lost or gained a genetic or epigenetic modification(s) that delayed/reduced Cre expression. If so, the data indicate that SMYD1 is not required for skeletal muscle function per se but is essential for maintaining the development and differentiation of skeletal muscle. This is supported by the Myf6-Cre data reported previously [12,13] and noted below.
Previously unpublished data obtained from shRNA knockdowns of SMYD1 in C2C12 myocytes indicated that SMYD1 repressed Calcineurin (CaN) regulation of slow fiber formation (L. Zhu, personal communication). Consistent with this was our previous observation that SMYD1 forms a stable complex with CaN [54]. We further showed that knockdown of SMYD1 led to dissociation of the complex and activation of CaN-mediated dephosphorylation [54]. However, and in apparent conflict with these conclusions, ablation of Smyd1 postnatally at 6 wk via Myf6-Cre produced mutant mice that exhibited loss of fast-twitch muscle markers [13].
The opportunity afforded by the rare adult MyoG-Cre/Smyd1 skeletal muscle survivors allowed us to address this apparent paradox by providing adequate amounts of dissectible, regional-specific muscle groups of known fiber types. These include the soleus muscle, which is composed primarily of slow muscle fibers; the tibialis anterior, primarily fast muscle fibers; and the gastrocnemius, which is a mixture of the two fiber types [55,56]. In skeletal muscle, there are four types of myosin heavy chains: a slow oxidative MyHC1, a fast oxidative MyHC2A, a fast/intermediate MyHC2X, and a fast glycolytic MYHC2B [55,56].
As shown in Figure 5C, at postnatal day 12 (P12), SMYD1 knockdown was robust within gastrocnemius and within the tibialis anterior muscles, but less effective in the soleus (probably because Cre-expression there was less efficient ( Figure 5C). Surprisingly, based on previous observations, we detected no expression of slow muscle markers (MyHC-1 and TPN1S) within CKO fast tibialis anterior muscles, nor in these markers (including MyHCIIA) in either in the tibialis anterior or the soleus. However, glycolytic fast muscle myosin (MyHC-IIB) was upregulated in CKO slow soleus ( Figure 5C). Alternatively, intermediate oxidative-glycolytic myosin (MyHC-IIX) was decreased in the slow soleus but increased in the mixed gastrocnemius and fast tibialis anterior fibers of Myog-Cre:Smyd1 CKOs ( Figure 5C). These results were confirmed in SMYD1-deficient soleus at P24 both at the RNA level ( Figure S5) and at the protein level for a pan antibody against MyHC-IIX that recognizes all slow twitch fibers ( Figure S6).
We suspect these phenomena were due to epigenetic silencing, as such, phenomena typically are not "all or none". This is evidenced by the mechanism of reversible methylation, which plays a role in many adaptive processes [57]. In the present context, such a process could lead to incremental decreases in Cre transgene expression from generation to generation. One might expect such alterations to be both stochastic (i.e., a random percentage acquires a change) and inheritable (i.e., once acquired, reversion is unlikely in subsequent generations).
These data, as with our previous reports [12,13], suggest that loss of SMYD1 results in slow-to-fast transformation of muscle fiber-type in the soleus. Physiologically, loss of slow-twitch, mitochondrial-rich fibers would lead to reduction in oxidative metabolism and greater fatigue. If mitochondrial-rich slow fibers are compromised, Myog-Cre:Smyd1 CKO mice would rely on type II a fast fibers to balance this deficiency. Accordingly, MyHC 2a and MyHC 2b muscle function should remain unaffected or perhaps even slightly increased on loss of SMYD1. Muscle fiber composition is a critical determinant of the capacity for mitochondrial fuel oxidation, ATP synthesis, and muscle endurance [56,58]. Thus, it would be informative to subject Myog-Cre:Smyd1 adult survivors to treadmill and distance running tests.

Concluding Remarks
We find it remarkable that SMYD1 deficiency leads to a number of multiple phenotypes dependent on the developmental stage and/or the Cre-recombinase employed for targeting [5][6][7][8][9][10][11][12][13]. Is there a single fundamental mechanism underlying this apparent heterogeneity? We submit that, yes, the diverse consequences of SMYD1 loss depends upon its bi-functionality as both a lysine (K) and a histone (H) MTase. Further diversification most surely rests on its ability to catalyze multiple histone marks, including H3K4m3, H3K4m2, and H4K20m2, as well as a growing list (currently up to nine) of confirmed KMTas targets [4,5]. While no additional "direct" targets were identified in this report, we judge it very likely that much, if not all, of the transcriptional changes documented in Figure 4 result from protein stabilization or other post-transcription functions. We suggest that identification of additional methylation targets and systematically assigning them into specific developmental stages of heart and skeletal muscle development will be key in establishing the basis for SMYD1-deficient heterogeneity.
While we feel that our study opened another new avenue for SMYD1 research, additional work is required to better interpret the complex phenotype. For example, a systematic histologic time course across each stage of prenatal development needs to be linked to that for full gene expression analyses. At present, results are compromised as only certain stages were analyzed, and there may be some discrepancy between the histology (most penetrant at E14.5) and the gene expression (analyzed primarily at E12.5), since the authors suggest an increase in BAT instead of muscle (at E14.5), but the gene expression shows a decrease in expression of BAT markers, Prdm16 and UCP1, at E12.5. There are also lots of data mentioned but not included. A full time course analysis will make it possible to fully interpret what appears to be an interesting and unusual phenotype.


extfieldsofstudy: Biology Medicine

text:
Single-cell full-length total RNA sequencing uncovers dynamics of recursive splicing and enhancer RNAs

Total RNA sequencing has been used to reveal poly(A) and non-poly(A) RNA expression, RNA processing and enhancer activity. To date, no method for full-length total RNA sequencing of single cells has been developed despite the potential of this technology for single-cell biology. Here we describe random displacement amplification sequencing (RamDA-seq), the first full-length total RNA-sequencing method for single cells. Compared with other methods, RamDA-seq shows high sensitivity to non-poly(A) RNA and near-complete full-length transcript coverage. Using RamDA-seq with differentiation time course samples of mouse embryonic stem cells, we reveal hundreds of dynamically regulated non-poly(A) transcripts, including histone transcripts and long noncoding RNA Neat1. Moreover, RamDA-seq profiles recursive splicing in >300-kb introns. RamDA-seq also detects enhancer RNAs and their cell type-specific activity in single cells. Taken together, we demonstrate that RamDA-seq could help investigate the dynamics of gene expression, RNA-processing events and transcriptional regulation in single cells.

1. The authors use NSR primers to reduce the representation of rRNAs that would otherwise dominate the sequencing output. This seems to be a critical step for the protocol to be usefulhowever, I assume that this process is not perfect. What is the percentage of reads in the sequencing libraries that correspond to ribosomal RNA? How does this compare to ribosomal depletion strategies involving probes for capture or degradation? (Note that quantification from sequencing data requires some care, as the many rRNA copies may preclude alignment and counting if done naively.) 2. Some more elaboration on the sample sizes would be desirable. How many diluted mESC RNA samples were generated with each protocol to calculate the CV2 values in Figure 1E? I would expect around 100 samples at least, given that real plate-based scRNA-seq experiments would be expected to use the entire plate. A related concern is the number of plates/batches that were used, not just the number of wells on each plate. scRNA-seq data is subject to considerable plate effects (see https://doi.org/10.1101/025528, https://doi.org/10.1038/srep39921), and I would strongly recommend examining not just the well-to-well variability, but also the plate-to-plate variability, e.g., when plates are processed on different days, with different operators, etc.
3. With respect to the conclusion that RamDA-seq provides full-length coverage, Figure 2A may be misleading as the profile depends on the true length of transcripts. One can imagine that, even for a perfect technique that gave completely uniform coverage, the profile would not be flat if the underlying distribution of full transcript lengths was uneven. I would suggest an additional plot where the transcript lengths are scaled to a fixed interval. This ensures a protocol yielding fulllength uniform coverage would also yield a flat profile.
4. The numbers on page 10 do not match up with the numbers in Figure 3B.
5. The linear regression to identify the sawtooth pattern for recursive splicing does not seem appropriate. In particular, I do not understand why the authors are using a dummy binary variable (page 52). It seems more logical to use the distance from the closest upstream RS site candidate as a real-valued covariate to model the coverage at each bin, given that coverage drops linearly with distance. A related issue is the claim that "the sawtooth wave pattern... weakened as differentiation progressed". While the decrease in coverage is clear, is this a specific reduction in recursive splicing or simply a general downregulation of the entire gene?
6. For Figure 5B, I would like to see separate coverage profiles for the forward and reverse strands. This will confirm that they are indeed eRNAs, which are characterised by inverted bimodality around the enhancer. The "background" is also constructed from random genomic regions -how were these chosen? Were these matched for GC content, mappability/sequenceability, DNA accessibility, etc.?
Obviously, it would be trivial to obtain low "background" coverage if the background regions were defined in difficult-to-map regions such as telomeres.
7. The read processing steps are very inconsistent across the different analyses. For some analyses, alignment was performed using HISAT2 (page 42), while for others, Bowtie2 was used (page 49). For quantification of gene expression, some analyses used Sailfish (pages 42, 45), others used Cufflinks (page 43) and others used featureCounts on aligned data (page 46). It would be advisable to standardize the quantification pipeline for all analyses, to avoid confusion and simplify the interpretation of the results. The use of Cufflinks seems particularly unnecessary, as I do not see many references to newly-discovered unannotated transcripts in the manuscript. 8. The selection of genes on page 48 is not coherent. Why do the authors identify the top HVGs within each phase, and then select the genes that are DE between phases? This will artificially encourage the separation of cells from different phases on the diffusion map. If the aim is to create an unbiased diffusion map, they should just use highly variable genes detected across the entire population of cells (from all phases). Otherwise, the authors should justify their choice of genes to construct a biased pseudotime.

MINOR:
1. Genomic contamination is an important consideration for RNA-seq protocol development. There is no harm in showing the data as a Supplementary on the last line of page 6.
2. Figure 1D uses a rather unconventional metric to examine the similarity between each scRNAseq library and the bulk data (the number of genes with counts >= 10 and |FC| <= 2). A more natural statistic would be the Pearson's correlation between the counts, which can be directly interpreted as the proportion of variance explained. I suspect this would give similar results without the need for the arbitrary threshold on the fold-change. Otherwise, it would be necessary to show similar results at different thresholds. Figure 5C is misleading as the expression values have simply been copy-andpasted to form a continuum. I would suggest colouring the copies differently to avoid giving the idea that there is information on the cycle number for each cell.

Supplementary
4. Page 49 makes references to the cell cycle phases, while the section refers to time series data. I presume that this is a mistake. 5. Consider using TREAT (specifically, the glmTreat function) to test against a log-fold change threshold, rather than post-hoc filtering on the log-fold changes on page 44. See https://support.bioconductor.org/p/62286/#62287 for a discussion of this topic by some of the edgeR authors.
Reviewer #2 (Remarks to the Author): In this paper Hayashi et al. propose RamDA-seq, a full length total RNA-seq method that shows near complete full-length transcript coverage, and, as the authors claim, especially high sensitivity to non-poly(A) RNA. The method was tested to reproduce dynamical trajectories of differentiating embryonic stem cells, revealing hundreds of dynamically regulated non-poly(A) transcripts. They also demonstrate that the data resulting from RamDA-seq can detect recursive splicing, enhancer RNAs and their cell-type specific activity in single-cells. This is a well written paper, with the data well assembled and presented and the results clearly stated, and we support its publication as we think that the single-cell RNA-seq community will greatly benefit from this method. However, we have a number of major and minor comments that we would like the authors to address before accepting the paper fully.
Major comments 1. Important details of how the RT enzyme works are missing from this paper. While we understand that it is described in greater detail in the other paper (Hayashi et al., NAR, in review) by the same authors, it would be good to add a supplementary note describing details about the underlying biology (e.g. how does strand displacement amplification work?) 2. The authors claim that "RamDA-seq protocol did not produce a library of DNA derived from genomic DNA or environmental DNAs (data not shown)". Can the authors provide some evidence of this claim? Supplementary Fig 3c shows that the proportion of intronic/intergenic reads is higher in RamDA-seq ( >40%) than other methods. How do we reconcile these two statements? Some commentary is essential. 3. For Fig. 1d and 1e, I believe a useful (and easy) control would be to perform RamDA-seq on 1 microgram RNA and compare its performance to the gold standards (rdRNA-seq and paRNA-seq) 4. It is not clear if the data shown on Fig. 2b is an average across multiple single cells for RamDAseq or coverage within a single cell. It would be useful to see what the coverage within example single cells looks like. Also it is not clear why Mdb1/Mdn1 (the name of the gene is mismatched between the figure and the text) was chosen? Can the authors comment on its significance? 5. Are non-poly(A) transcripts enriched specifically in Clusters 1,2 and 5 in Fig. 3d. If so, why? 6. Related to Fig. 3, a discussion on the biological significance of Neat1 is missing. The authors should relate to what was published 7. The authors claim that "the sawtooth wave pattern in the first intron of Cadm1 weakened as differentiation progressed, whereas the pattern in the second intron of Robo2 strengthened". But this could merely be the result of changing expression levels. In other words, the data as presented are consistent with the fact that the recursive splicing stays exactly the same, it's only that it cannot be detected when global expression levels are low. Can the authors comment on the sensitivity of RS detection as a function of expression levels? 8. Also, can the authors detect heterogeneity across single cells in the presence/absence of RS within a specific gene? 9. Lastly, but importantly, given the high proportion of rRNA, intronic and intergenic reads in RamDA-seq, can the authors comment on the depth per cell needed to get useful information out?
Minor comments 1. The often used phrase "RamDA-seq is the single-cell RNA-seq method" very strange. This does not convey any meaning and I am sure is an error. Stated this way, it appears that there is no other single-cell RNA-seq method other than RamDA-seq. I would suggest a rephrasing as follows: "Altogether, based on the substantial resemblance …….., we conclude that RamDA-seq is a robust single-cell total RNA-seq method" 2. For comprehensiveness, it would be good if Fig. 2b involved SUPER-seq. Also can the authors share a few more representative examples like Fig. 2b (which is impressive) in the supplement? 3. The wording of Fig. 2c within the text is not clear. What does "longer range of transcripts" mean? I guess the authors are referring to a longer range of transcripts binned by lengthclarification would help.

RESPONSE TO REVIEWERS
Original reviewer comments appear below in italics with responses in blue.

RESPONSE TO REVIEWER 1:
We are thankful for the reviewer's insightful comments, which have helped us to significantly improve the paper.  with 10 pg and approximately 10-25% with cells. Despite having a higher rRNA read ratio than rdRNA-seq, RamDA-seq showed high sensitivity and reproducibility comparable to rdRNA-seq (please refer to the revised version of Supplementary Fig. 4).
Moreover, to our knowledge, there is no scRNA-seq method using probes for rRNA depletion. Developing scRNA-seq methods using rRNA depletion is very difficult because loss of RNA in the RNA purification step is extremely severe. Indeed, the minimum input RNA amount for standard rRNA depletion RNA-seq kits is approximately 25 ng of total RNA, which is far greater than the typical RNA amount in single cells (1-10 pg). Thus, we chose to use NSRs to reduce rRNA reads (Armour et al., Nat Methods, 2009;Ozsolak et al., Genome Research, 2010). Therefore, we have added the following text to the Introduction section (p.

4, line 3):
In addition, how to reduce the sequence derived from ribosomal RNAs (rRNA) that accounts for most of the total RNA is a major task for establishing single-cell total RNA-seq. This issue is encountered because single-cell RNA-seq uses a trace amount of total RNA as a template, and it is difficult to apply general rRNA-depletion methods that cause loss of RNA.
Finally, we have changed the sentence "the RamDA-seq library still contains a relatively high proportion (10-25%) of rRNA sequences ( Supplementary Fig. 3b)." in the Discussion section  Figure 1E? I would expect around 100 samples at least, given that real plate-based scRNA-seq experiments would be expected to use the entire plate. A related concern is the number of plates/batches that were used, not just the number of wells on each plate. Response: We agree that this point requires clarification. The number of diluted RNA samples was n= 3 or 4 for each scRNA-seq method. We showed the sample sizes in the revised version of Figure 1e.
In response to the concern raised by this reviewer, we also prepared ~96 diluted RNA samples (corresponding to the entire plate) with RamDA-seq, C1-RamDA-seq, and C1-SMART-seq v4 and calculated CV^2 values. RamDA-seq and C1-RamDA-seq showed high reproducibility (low CV^2 values) in most expression level bins ( Supplementary Fig. 4f).
Moreover, in accordance with the reviewer's comment, we evaluated the plate-to-plate variability of RamDA-seq. Thus, we have added new experimental data to Supplementary  Supplementary Fig. 16a). Plate-to-plate variability was evaluated using the correlation of expression levels and PCA. Although variabilities between plates were larger than variabilities within plates, plate-to-plate variability was much smaller than the variability between cell types (revised version of Supplementary Fig. 16d-f), and the proportion of variance explained in PCA was 0.6% (revised version of Supplementary Fig. 16g-h). We have added the following sentence to the Discussion section (p. 18, line 9): We also evaluated plate-to-plate variability (batch effect) of RamDA-seq ( Supplementary Fig. 16a). Although variabilities between plates were larger than variabilities within plates, plate-to-plate variability was much smaller than the variability between cell types, and the proportion of variance explained in principle component analysis (PCA) was 0.6% (Supplementary Fig. 16b-h). While reviewing our data analysis pipeline (as suggested by Reviewer #1's Major comment 7), we found that the plot in the previous version of Figure 2a was created using RefSeq annotation. To ensure the consistency with the other parts of the manuscript, we replaced the previous plot with a plot created using GENCODE (revised version of Fig. 2a and Supplementary Fig. 4i). This revision did not change the conclusion.
On the other hand, we agree that additional information on the read coverage as the reviewer suggested would be valuable. Thus, we have added the transcript coverage of scRNA-seq methods against relative transcript positions to the revised version of Supplementary Figure. 7d. RamDA-seq showed uniform transcript coverage, similar to rdRNA-seq (revised version of Supplementary Fig. 7d).
Comment 4: The numbers on page 10 do not match up with the numbers in Figure 3B.
Response: The reviewer's comment is correct. We corrected the numbers in the revised manuscript (p. 10, line 6). Response: We thank the reviewer for providing a suggestion on how to improve our analysis.
In response to the first point in the reviewer's comment, we implemented the proposed model, which uses the distance from the closest upstream RS site candidate as a realvalued covariate and fitted the model against mapping data on Robo2, Cadm1, and Magi1.
Additionally, we performed the Cox test for a non-nested model to evaluate the significance of fitting (briefly, the proposed model is selected against the baseline model if the difference in AIC (dAIC) is positive and p < 1e-5). Unfortunately, the proposed model failed to detect the sawtooth pattern in Robo2 and Magi1 (Letter Figure 1). The proposed model is forced to predict the same values at an RS site and its upstream annotated splicing junction; thus, it is likely to lack the flexibility to detect the sawtooth pattern of recursive splicing. Based on the new results, we would like to keep the previous version of the results.
Regarding the second point in this reviewer's comment and Reviewer #2's major comment 7, we reanalyzed the read coverages in detail. We found that the sawtooth wave pattern itself was robustly observed (revised version of Supplementary Fig. 13, middle panels) and that the 'height' of the pattern was associated with the expression levels of host genes (revised version of Supplementary Fig. 13, left panels). Therefore, we have changed the sentence "Interestingly, the sawtooth wave pattern in the first intron of Cadm1 weakened as differentiation progressed (Fig. 4b), whereas the pattern in the second intron of Robo2 strengthened (Fig. 4c). The sawtooth patterns were also confirmed using averaged RamDAseq read coverage for each time point (Fig. 4d-e). These results demonstrate that RamDAseq can detect RS at the single-cell level and exemplify the potential of RamDA-seq to detect dynamics in RNA processing events." in the revised Results section (p. 11, line 27) as follows: The Comment 6: For Figure 5B, I would like to see separate coverage profiles for the forward and reverse strands. This will confirm that they are indeed eRNAs, which are characterised by inverted bimodality around the enhancer. The "background" is also constructed from random genomic regions -how were these chosen?
Obviously, it would be trivial to obtain low "background" coverage if the background regions were defined in difficult-to-map regions such as telomeres.
Response: We agree that additional information on the coverage profiles for the forward and reverse strands as the reviewer suggested would be valuable. Regrettably, because RamDA-seq is a non-stranded sequencing method as mentioned in the Discussion section "It is also important to achieve strand-specific sequencing in RamDA-seq to distinguish overlapping transcripts." in the previous version of the manuscript, we are unable to perform the analysis. Thus, we have retained the original version of Figure 5b in the manuscript.
The background regions in the previous version of manuscript were selected randomly from intergenic regions, but we ignored GC content and other characteristics. We agree with the reviewer's concern that low-mappability regions would yield low read coverage. To eliminate this unwanted possibility, we prepared several other sets of background regions that have enhancer-like sequence characteristics and/or chromatin accessibility and made aggregation plots again. The read coverage in all background region sets was steadily and remarkably lower than that in ES-active enhancer regions (revised version of Supplementary Fig. 14c and 15c).  For time-series data, (1) we have changed the sentence "We identified 5,342 such transcripts, including 423 non-poly(A) transcripts (Fig. 3b), and divided the 5,342 transcripts into seven clusters based on expression patterns (Fig. 3c)." (p. 10, line 6) as follows: We identified 7,580 such transcripts, including 458 non-poly(A) transcripts (Fig. 3b), and divided the 7,580 transcripts into seven clusters based on expression patterns (Fig. 3c).
(3) The clusters in which Pou5f1 and PrE maker genes were included were changed from (4) The number of dynamically regulated eRNAs decreased (Fig. 5d) and, for interpretability, we set the number of clusters to five in the subsequent clustering (Fig. 5d,e). Accordingly, the clusters in which GATA4 motif was enriched were changed from 4 and 6 to 3 and 5, respectively. Thus, we have made the corresponding changes in the revised manuscript (p. 14, line 17) and have changed the sentence "We found 1,441 dynamically regulated eRNAs (FDR < 0.01) in seven clusters: three down-regulated clusters, three transiently up-regulated clusters, and one late up-regulated cluster (Fig. 5d-e)." (p. 14, line 11) and as follows: We found 1,338 dynamically regulated eRNAs (FDR < 0.01) in five clusters: two down-regulated clusters, two transiently up-regulated clusters, and one late upregulated cluster (Fig. 5d,e).
For analyses of recursive splicing and eRNA, we have established specialized pipelines based on the mapping data (by HISAT2) because, unlike expression level quantification, these analyses need read coverage data along genomic coordinates. To avoid confusion, we added sentences to the revised Methods section (p. 36, line 16; p. 38, line 13).
Cufflinks was applied to the bulk RNA-seq data (not single-cell RamDA-seq data) for transcriptome assembly. In the previous version of manuscript, the resulting gene models were used to demonstrate the potential of RamDA-seq for detecting annotated and newly discovered non-poly(A) transcripts and to define a set of enhancer RNAs with enhancer-like histone modification. Thus, we disagree that the use of Cufflinks is unnecessary and would rather retain the previous version of the analysis. We hope that this issue does not cause any unnecessary trouble. Response: We appreciate the reviewer's concerns regarding this point. In response to the reviewer's comment, we conducted an 'unsupervised' diffusion map analysis using highly variable genes selected regardless of the derived cell-cycle phases. RamDA-seq could detect differences among cells in different cell-cycle phases (revised version of Supplementary Fig. 8a).
Our aim in the previous version of manuscript was to reconstruct a pseudotime specifically reflecting a subjective time along the cell cycle; thus, we chose the gene selection strategy in the previous version of manuscript. We would prefer to retain the previous version of manuscript. To clarify our aim and avoid confusion, we added descriptions of diffusion maps using both the 'unsupervised' and 'supervised' gene selection methods (revised version of Supplementary Fig. 8a,b), and we have changed the sentences "This analysis divided the cells according to their derived cell-cycle phases, although the cells in each phase showed variance ( Supplementary Fig. 5a). Based on this observation, we hypothesized that the RamDA-seq data could be used to reconstruct a 'subjective time' for cells along the cell cycle. To address this possibility, we used DC1 as a proxy of subjective time for each cell along the cell cycle and searched for transcripts oscillating along DC1 by fitting a sine function (false discovery rate (FDR) < 0.01; Materials and Methods)." in the revised manuscript (Supplementary Note 5, p.6 line 12) as follows: The cells were largely separated according to the derived cell-cycle phase ( Supplementary Fig. 8a). This result indicates that RamDA-seq could detect variability among different cell-cycle phases within a single cell type.
Based on this observation, we hypothesized that the RamDA-seq data could be used to reconstruct a 'subjective time' for cells along the cell cycle. To address this possibility with use of information from a cell sorter, we first selected transcripts showing high variance among G1, S, and G2M, and performed diffusion map analysis. This procedure likely reconstructed a subjective time in each cell; the cells formed different clusters according to their derived cell, while the cells in each phase showed variance (Supplementary Fig. 8b). We then used DC1 as a proxy of subjective time for each cell along the cell cycle and searched for transcripts oscillating along DC1 by fitting a sine function (false discovery rate (FDR) < 0.01; see Methods).
We have noticed that the description of how to select highly variable genes in the previous version of manuscript was misleading. To avoid confusion, we have changed the sentence "For each transcript among the expressed genes, we first calculated the high CV of the averaged TPM within the G1, S, and G2M samples" in the revised manuscript (p. 34, line 21) as follows: For each transcript among the expressed genes, we first calculated averaged TPM across cells in each cell cycle phase and then calculated the CV using the averaged TPM values for G1, S, and G2M.

Minor Comments Comment 1: Genomic contamination is an important consideration for RNA-seq protocol development. There is no harm in showing the data as a Supplementary on the last line of page 6.
Response: In accordance with the reviewer's comment, we added quality control data for RamDA-seq library DNA, including the data from the RT enzyme minus and non-template control experiments in the revised manuscript (revised version of Supplementary Fig. 3).
Contamination of genomic or environmental DNA did not occur.
Comment 2:. Figure 1D uses a rather unconventional metric to examine the similarity between each scRNA-seq library and the bulk data (the number of genes with counts >= 10

and |FC| <= 2). A more natural statistic would be the Pearson's correlation between the counts, which can be directly interpreted as the proportion of variance explained. I suspect this would give similar results without the need for the arbitrary threshold on the fold-change.
Otherwise, it would be necessary to show similar results at different thresholds.
Response: As suggested by the reviewer, we showed the Pearson correlation coefficients between each scRNA-seq method and bulk rdRNA-seq (revised version of Supplementary   Fig. 4h). RamDA-seq and C1-RamDA-seq were most similar to bulk rdRNA-seq. We also showed the number of detected transcripts for each scRNA-seq method without any threshold for fold change (revised version of Supplementary Fig. 4d), which provided similar results to those in Figure 1d.

Comment 3:. Supplementary Figure 5C is misleading as the expression values have simply been copy-and-pasted to form a continuum. I would suggest colouring the copies differently to avoid giving the idea that there is information on the cycle number for each cell.
Response: Thank you for your helpful suggestions for data visualization. To avoid confusion, we showed the data points once, colored the copies differently, changed the x-axis scale so that ticks corresponded to the boundary of cycles, and added arrows (revised version of Supplementary Fig. 5d,e).

Comment 4: Page 49 makes references to the cell cycle phases, while the section refers to time series data. I presume that this is a mistake.
Response: The reviewer's comment is correct. To correct for this mistake, we have changed the sentences "We conducted a diffusion map analysis as follows. First, we selected expressed transcripts with a TPM of at least 1 in at least 10% of the cells. For each transcript among the expressed genes, we first calculated the high CV of the averaged TPM within the G1, S, and G2M samples. For each cell cycle phase, we next selected the top 5000 high-CV transcripts with an averaged TPM that was higher than those for the other phases. For the selected transcripts, we performed diffusion map analysis using the 'destiny' R package. We used DC1 as the pseudotime, i.e., we treated the expression levels as a function of DC1." in the revised manuscript (p. 35, line 18) as follows: We conducted a diffusion map analysis as follows: (1) (Fig. 2e, Fig. 3b,c, Supplementary Fig. 8c, Supplementary Fig. 10

RESPONSE TO REVIEWER 2:
We wish to express our appreciation for the reviewer's insightful comments, which have helped us to significantly improve the paper.

Major Comments
Comment 1: Important details of how the RT enzyme works are missing from this paper.
While we understand that it is described in greater detail in the other paper (Hayashi et al.,

NAR, in review) by the same authors, it would be good to add a supplementary note describing details about the underlying biology (e.g. how does strand displacement amplification work?)
Response: As suggested by the reviewer, to describe the details of the RT-RamDA reaction, we added a section "The principle of RT with random displacement amplification (RT-RamDA)" to the revised Supplementary Note 1.
Comment 2: The authors claim that "RamDA-seq protocol did not produce a library of DNA derived from genomic DNA or environmental DNAs (data not shown)". Can the authors provide some evidence of this claim? Supplementary Fig 3c shows that the proportion of intronic/intergenic reads is higher in RamDA-seq ( >40%) than other methods. How do we reconcile these two statements? Some commentary is essential.
Response: We confirmed that the RT enzyme minus and non-template control experiment did not produce library DNA (revised version of Supplementary Fig. 3), indicating that contamination of genomic or environmental DNA did not occur. Therefore, the QC metrics of scRNA-seq methods ( Fig. 1 and revised version of Supplementary Fig. 4) indicate that the observed high proportion of intronic and intergenic reads, which is also a hallmark of rdRNAseq, reflect the resemblance of RamDA-seq to rdRNA-seq. We believe that these reads were derived from transcripts, not endogenous genomic DNA.
In accordance with Reviewer #1's major comment 7, we changed the data set of transcript annotation from Refseq to GENCODE to ensure consistency through the manuscript. Thus,  Fig. 1d and 1e, I believe a useful (and easy) control would be to perform RamDA-seq on 1 microgram RNA and compare its performance to the gold standards (rdRNA-seq and paRNA-seq) Response: As suggested by the reviewer, we compared RamDA-seq and rdRNA-seq prepared with equal amounts of input RNA. Since RT-RamDA and a Nextera XT library preparation kit have limitations regarding the input amount of template, it is not easy to compare its performance with standard RNA-seq using the same amount of input RNA. Therefore, considering the maximum amount of input RNA, we prepared RamDA-seq libraries using 1 ng of RNA and compared their performance with that of rdRNA-seq libraries using 1 or 10 ng of RNA. Compared with rdRNA-seq (1 ng), RamDA-seq (1 ng) showed a higher proportion of rRNA reads (revised version of Supplementary Fig. 5 a), which was the case when comparing RamDA-seq (10 pg) with rdRNA-seq (1 μg) (revised version of Supplementary Fig. 4b,c and revised version of Supplementary Fig. 12c). Despite the difference, RamDA-seq (1 ng) and rdRNA-seq (1 ng) showed comparable sensitivity (number of detected transcripts) and reproducibility (SCC, PCC and CV^2) (revised version of Supplementary Fig. 5b-d). Moreover, the correlation between RamDA-seq (1 ng) and rdRNA-seq (1 ng or 10 ng) was as high as that between RamDA-seq (10 pg) and rdRNAseq (1 μg) (revised version of Supplementary Fig. 4h and revised version of Supplementary   Fig. 5c). These results support our previous conclusion that RamDA-seq is a single-cell total RNA-seq method that resembles rdRNA-seq. We have added the following text to the Results section (p. 7, line 23): We also compared large-volume-inputted RamDA-seq with rdRNA-seq using 1 ng of RNA. When we compared RamDA-seq and rdRNA-seq with the same input amount, there was hardly any difference in sensitivity and reproducibility except for the contamination rate of rRNA ( Supplementary Fig. 5).

Comment 4: It is not clear if the data shown on Fig. 2b is an average across multiple single cells for RamDA-seq or coverage within a single cell. It would be useful to see what the coverage within example single cells looks like. Also it is not clear why Mdb1/Mdn1 (the name of the gene is mismatched between the figure and the text) was chosen? Can the authors comment on its significance?
Response: We completely agree and apologize for this omission. We originally showed the read coverage in 10 pg data (n= 1) (which corresponds to a single cell) (Fig. 2b). To clarify this point, we added the numbers of samples of each scRNA-seq (n= 1) to the revised version of Figure 2b.
We first selected 25 expressed long transcripts as transcripts whose TPM >= 5 in rdRNAseq in mESC and whose transcript length >= 10 kb. We further selected Mdn1, which has the highest number of exons (102 exons) among the expressed long transcripts, for visualization. In response to Reviewer #2's minor comment 4, we showed the read coverage of 25 other expressed long transcripts (revised version of Supplementary Data 2 and 3 for 10 pg RNA and mESC data, respectively). We also added text to the legend of Figure 2b in the revised manuscript as follows: We selected Mdn1 as the gene with the highest number of exons (102 exons) in the 25 genes with length >= 10 kb and TPM >= 5 in rdRNA-seq results. Fig. 3d.

If so, why?
Response: In response to Reviewer #1's major comment 7, we changed the expression quantification method and reanalyzed the data. Although the proportion of non-poly(A) transcripts in each cluster changed, non-poly(A) transcripts were included in all clusters irrespective of the quantification methods (revised version of Fig. 3c). This result suggests that non-poly(A) transcripts might contribute to various cell functions. Accordingly, we have changed the sentence "Closer inspection of clusters 1, 2, and 5 revealed that these three clusters contained many non-poly(A) transcripts, as exemplified by two unannotated intergenic non-poly(A) transcripts and Hist1h1a (Fig. 3d)." in the revised manuscript (p. 10, line 9) as follows: The dynamically regulated non-poly(A) transcripts were spread in all clusters with various expression patterns, suggesting that non-poly(A) transcripts are involved in various cell functions. Using single-cell RT-qPCR, we validated the observed expression changes in several of the non-poly(A) transcripts, including two unannotated intergenic non-poly(A) transcripts (clusters 1 and 2) and Hist1h1a (cluster 5) ( Fig. 3d and Supplementary Fig. 12g). Fig. 3, a discussion on the biological significance of Neat1 is missing.

The authors should relate to what was published
Response: In accordance with the reviewer's comment, we added a new paragraph (p. 17, line 1) to the revised version of the Discussion as follows: Neat1 is an architectural component of paraspeckle nuclear bodies38, which regulate gene expression via capture of A-to-I edited mRNAs39 and transcription factors40, and is required for corpus luteum formation and establishment of pregnancy in mice41. The long non-poly(A) isoform Neat1-001, not the short poly(A) isoform Neat1-002, is essential for the formation of paraspeckles42. Although the two isoforms are transcribed from the same promoter, they show different expression patterns, and Neat1-001 is expressed only in a small subpopulation of cells in adult mouse tissues26. Therefore, distinguishing the expression of the two isoforms of Neat1 at the single-cell level is critical for studying their functions. In this study, RamDA-seq's sensitivity to non-poly(A) transcripts and full-length transcript coverage distinguished the expression of the two isoforms and revealed dynamic and differential regulation of the long non-poly(A) isoform (Fig. 3). These results suggest that RamDA-seq could be beneficial for investigation of temporal and spatial expression patterns of long non-poly(A) RNAs in single cells. Response: In response to the reviewer's comment and Reviewer #1's major comment 5, we reanalyzed the read coverages in detail. The sawtooth wave pattern itself was robustly observed (revised version of Supplementary Fig. 13, middle panels), and the 'height' of the pattern was associated with host gene expression levels (revised version of Supplementary   Fig. 13, left panels).
To estimate the sensitivity of RS detection as a function of host gene expression levels, we then performed a simulation in which intronic read coverage data were repeatedly subsampled (100 times per subsampling fraction), and the success rate of RS detection was calculated (Methods section). Then, the aggregated expression levels of host genes at each time point were mapped to the estimated sensitivity curve (revised version of Supplementary   Fig. 13, right panel). Note that the fraction of intronic reads was converted to host gene expression levels (TPM). The result shows that the sensitivity of RS detection is a sigmoidlike monotonically increasing function of host gene expression (revised version of Supplementary Fig. 13, right panel) and coincides with RS detection at each time point for Robo2, Cadm1, and Magi1 (revised version of Supplementary Fig. 13, middle panel). The detection limit varied from gene to gene. Altogether, we were able to estimate the sensitivity of RS detection. Based on these results, we have changed the sentences "Interestingly, the sawtooth wave pattern in the first intron of Cadm1 weakened as differentiation progressed (Fig. 4b), whereas the pattern in the second intron of Robo2 strengthened (Fig. 4c). The sawtooth patterns were also confirmed using averaged RamDA-seq read coverage for each time point (Fig. 4d-e). These results demonstrate that RamDA-seq can detect RS at the single-cell level and exemplify the potential of RamDA-seq to detect dynamics in RNA processing events." in the revised manuscript (page 11, line 27) as follows: The averaged RamDA-seq read coverage for each time point showed that the sawtooth wave patterns persisted at all time points for Cadm1 and Magi1, whereas the pattern was observed at only 48 and 72 h for Robo2 (Fig. 4d,e,   Supplementary Fig. 13). The height of the sawtooth wave pattern was associated with the expression level of host genes at all time points ( Supplementary Fig. 13), while the pattern was not observed when gene expression was hardly detected ( Supplementary Fig. 13). A simulation-based estimation of the sensitivity of RS detection showed that RS was robustly detected when host genes were sufficiently expressed (transcript per million (TPM) > 1), corroborating the above observations (Methods section; Supplementary Fig. 13, right panels). To report these results, we have added a new paragraph to the Results (p. 12, line 9) section as follows: Next, we attempted to address whether RamDA-seq can detect RS even in each single cell. Based on the RS detection sensitivity estimated above, we selected cells with sufficient intronic reads (RS detection probability > 0.95) (Supplementary Data 5, Methods section). We fitted linear regression models against the RamDA-seq read coverages of each single cell in Cadm1, Robo2, and Magi1. RS was detected in a subpopulation of cells (71 of 149 cells for Cadm1, 12 of 54 cells for Magi1, 1 of 1 cell for Robo2) although many cells in which RS was not detected also appeared to show the sawtooth pattern. However, interestingly, some other cells showed monotonically decreasing patterns, which correspond to 'normal' splicing (Fig. 4a). The monotonically decreasing pattern was also observed even when we filtered cells with a more stringent threshold of intronic reads (Supplementary Data 5). These observations raise the possibility of cell-to-cell heterogeneity in RS. Therefore, further investigation is needed to reveal the mechanisms and significance of the observed heterogeneity in RS. Collectively, these results demonstrate that RamDA-seq can detect RS in single cells.
We have also added a new paragraph to the Discussion section (p. 17, line 15) in the revised manuscript as follows: Unexpectedly, we observed cell-to-cell heterogeneity in read coverage patterns around RS sites, suggesting that some cells showed recursive splicing, and other cells showed normal splicing (Supplementary Data 5). These results indicate that RamDA-seq can detect cell-to-cell heterogeneity in RS and could help to address the mechanisms and relationship between transcription and splicing. Toward these goals, several important challenges remain. Given that some cells in which RS was not detected showed weak sawtooth patterns, RamDA-seq highlights the limitation of the current linear regression model used to detect RS in this study and the need for further improvement in computational methods to robustly detect RS using single-cell data. Another challenge is to experimentally and computationally distinguish biological and technical variabilities in RS at the single-cell level. We will address these challenges in the future.

Minor Comments
Comment 1: The often used phrase "RamDA-seq is the single-cell RNA-seq method" very strange. This does not convey any meaning and I am sure is an error. Stated this way, it appears that there is no other single-cell RNA-seq method other than RamDA-seq. I would suggest a rephrasing as follows: "Altogether, based on the substantial resemblance …….., we conclude that RamDA-seq is a robust single-cell total RNA-seq method" Response: In accordance with the reviewer's comment, we changed the description throughout the revised manuscript.
Comment 2: For comprehensiveness, it would be good if Fig. 2b involved SUPER-seq. Also can the authors share a few more representative exampleslike Fig. 2b (which is impressive) in the supplement?
Response: In accordance with the reviewer's comment, we added SUPeR-seq to the comparison in the revised version of Figure 2b. We also showed other examples in the revised version of Supplementary Data 2 and 3 (for 10 pg RNA and mESC data, respectively). Fig. 2c within the text is not clear. What does "longer range of transcripts" mean? I guess the authors are referring to a longer range of transcripts binned by lengthclarification would help.

Comment 3: The wording of
Response: In accordance with the reviewer's comment, we have changed the sentence "In addition, the fraction of exonic regions covered by the reads also indicated that RamDA-seq covered a wider range of transcripts than did the other methods ( Fig. 2c and Supplementary   Fig. 4)." (p. 8, line 20) as follows: In addition, the fraction of exonic regions covered by the reads indicated that RamDA-seq covered a higher fraction of exonic regions than did the other methods in all length bins (Fig. 2c, Supplementary Fig. 7a-c  Supplementary Figure 4: What is the difference between % of rRNA reads in FASTQ and % of mapped reads in rRNA sequences? This was not apparent from reading the caption or methods. I would have expected the latter to be larger than the former, given that rRNA reads must map (otherwise they would not be identified as mapped reads) but other reads may or may not map.
Supplementary Figure 16: This looks good, but I'm not clear on how the authors concluded that "variabilities between plates were larger than variabilities within plates". The within-plate correlations in D and F look very similar to the between-plate correlations, which suggestsencouragingly -that the plate effects are modest at worst.
Reviewer #2 (Remarks to the Author): I have reviewed the author's responses carefully and am satisfied. I recommend this paper for publication.

RESPONSE TO REVIEWER 1:
Original reviewer comments appear below in italics with responses in blue.
We wish to express our appreciation for the reviewer's comments, which have helped us to significantly improve the paper. Response: We thank the reviewer for the encouraging comment. We changed the sentences "Although variabilities between plates were larger than variabilities within plates, plate-toplate variability was much smaller than the variability between cell types, and the proportion of variance explained in principle component analysis (PCA) was 0.6% (Supplementary Fig.   16b-h)." in the Discussion section as follows:

Supplementary
The variabilities between plates were similar to the variabilities within plates, suggesting that the plate effects are modest at worst. Plate-to-plate variability was much smaller than the variability between cell types, and the proportion of variance explained in principal component analysis (PCA) was 0.6% ( Supplementary Fig. 16bh).