SAFB associates with nascent RNAs and can promote gene expression in mouse embryonic stem cells

Cell culture

Male mouse ESCs that express doxycycline-inducible Xist from the Rosa26 locus (derivation described in Trotman et al. (2020) were grown in DMEM (Gibco) supplemented with 15% Qualified Fetal Bovine Serum (Gibco), 1% Pen/Strep (Gibco), 1% l-Glutamine (Gibco), 1% Non-Essential Amino Acids (Gibco), 100 µM betamercaptoethanol (Sigma), and 0.2% LIF. Cells were maintained in incubators set at 37°C and 5% CO₂. Media was replaced daily.

Generation of WT parent ESCs used for this study

To generate the WT parent ESC line from which we ultimately deleted SAFB and SAFB2, full-length Xist-expressing ESCs from Trotman et al. (2020) were deleted of their hygromycin B resistance gene via Lipofectamine transfection of a plasmid expressing FlpE (Addgene #20733) (Beard et al. 2006). An amount of 5 µg of FlpE recombinase was mixed with 5 µL of P3000 reagent, 7.5 µL of Lipofectamine 3000 reagent, and Opti-MEM media (Gibco #31985-070) to a total volume of 250 µL. The reagents were incubated for 5 min at room temperature before being added to cells with fresh media. After 24 h, cells were pulsed with puromycin (2 µg/mL) for 72 h. Ninety-six hours after transfection, ESCs were trypsinized to single-cell suspension and plated onto irradiated fibroblast feeder cells (500–2000 cells/10 cm plate) until individual colonies were visible by eye (4–5 d). Individual colonies were then selected and grown in individual wells for genotyping. After genotyping, candidate clonal colonies underwent hygromycin B (50 µg/mL) selection to verify loss of resistance. Genotyping primers used are in Supplemental Table S10.

Generation of SAFB/2 knockout ESCs

sgRNAs to delete SAFB and SAFB2 were designed to the mm10 genome using CRISPOR with the specifications 20 bp-NGG—Sp Cas9, Sp Cas9-Hf1, eSp Cas9 1.1 (Concordet and Haeussler 2018). sgRNA sequences are found in Supplemental Table S10. Guides were cloned into the pX330 plasmid as specified in Cong et al. (2013) (Addgene plasmid #42230). To delete SAFB and SAFB2, Parent ESCs were seeded at 0.5 × 10⁶ cells per well in a six-well plate. The following day, the cells were transiently transfected using Lipofectamine 3000 (Invitrogen L3000-015): 800 ng of sgRNA plasmid pool and 200 ng puro resistant GFP plasmid (1 µg total) were mixed with 2 µL P3000 reagent, 7.5 µL Lipofectamine 3000 reagent and Opti-MEM media (Gibco #31985-070) to a final volume of 250 µL. The reagents were incubated for 5 min at room temperature before being added to cells with fresh media. After 24 h, cells were pulsed with Puromycin (2 µg/mL) for 48 h. After puromycin selection, cells were trypsinized to single cells and plated onto irradiated fibroblast feeder cells (500–2000 cells/10 cm plate) until individual colonies were visible by eye (4–5 d). Individual colonies were then selected and grown in individual wells for genotyping.

The two WT and two DKO lines that were selected for further study, along with the WT Parent line, were then rendered dox-inducible by transfection of the rtTA-expressing plasmid described in Kirk et al. (2018). One day prior to transfection, parent and DKO ESCs were seeded at 0.5 × 10⁶ cells per six-welled well. The following day, 500 ng of rtTA plasmid and 500 ng of transposase (1 µg total DNA) were mixed with 2 µL P3000 reagent, 7.5 µL Lipofectamine 3000 reagent and Opti-MEM media (Gibco #31985-070) to 250 µL. The reagents incubated for 5 min at room temperature before being added to cells with fresh media. After 24 h, cells underwent G418 selection (50 µg/mL) for 12 d.

cDNA expression plasmids

Plasmids expressing full-length or truncated versions of SAFB, GFP, and HNRNPU were designed in silico based on existing vector backbones from Schertzer et al. (2019b) and synthesized by Genewiz. All plasmids have been deposited into Addgene.

Generation of cDNA-expressing ESCs

One day prior to transfection, ESCs were seeded at 0.5 × 10⁶ cells per well of a six-well plate. The following day, 850 ng of cDNA plasmid and 150 ng of piggyBac transposase from Kirk et al. (2018) were mixed with 2 µL P3000 reagent, 7.5 µL Lipofectamine 3000 reagent, and Opti-MEM media (Gibco #31985-070) to a final volume of 250 µL. The reagents were incubated for 5 min at room temperature before being added to cells with fresh media. After 24 h, cells underwent hygromycin B selection (50 µg/mL) for 1 wk.

HNRNPU-expressing ESCs were then transfected with the rtTA from Kirk et al. (2018) as described above. The DKO cells in which GFP and SAFB cDNA constructs were introduced had previously been stably transfected with rtTA.

PCR

Genomic DNA was collected from 0.8 × 10⁶ cells with 500 µL lysis buffer (100 mM Tris-HCl pH 8.1, 0.5 mM EDTA pH 8.0, 200 nM NaCl, 0.2% SDS) + 80 µL Proteinase K (Denville) + 8 µL linear acrylamide (Thermo Fisher) and incubated at 55°C overnight. Twice the volume of ice cold 100% ethanol was added. Samples were then vortexed and rotated end-over-end at 4°C for 15 min. Samples were spun at max speed for 5 min at 4°C. The lysis buffer/ethanol mixture was then removed, and the DNA pellet was washed with 70% ethanol, after which the DNA pellet was resuspended in 1× TE (10 mM Tris pH 8.0, 1 mM EDTA) and incubated overnight at 56°C. DNA concentration was measured via Nanodrop and diluted to 50 ng/µL. PCR was performed with ChoiceTAQ (Denville CB4050) as follows: 25 µL PCR reaction mixture (2.5 µL 10× PCR reaction buffer, 0.2 µL 10 mM dNTPs, 0.25 µL 100 µM primers, 0.25 µL Choice TAQ polymerase, 3 µL DNA template [50 ng/µL], and ddH₂O to 25 µL) ran in Bio-Rad C1000 Touch or T100 thermocycler (initial denaturation at 95°C for 3 min; 25 cycles of 95°C for 30 sec, 58°C–62.5°C annealing for 30 sec, and 72°C for 30–45 sec extension time). PCR primers and conditions are in Supplemental Table S10.

RT-qPCR

Equal amounts of RNA (0.5–1 µg) were reverse transcribed using the High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific #4368813) with the random primers provided, and then diluted with 30 µL 1× TE. For RIP RT-qPCR, 2 µL of the eluted sample was used in RT reactions. An amount of 10 µL qPCR reactions were performed using iTaq Universal SYBR Green (Bio-Rad) and custom primers on a Bio-Rad CFX96 system with the following thermocycling parameters: initial denaturation at 95°C for 10 min; 40 cycles of 95°C for 15 sec, 60°C for 30 sec, and 72°C for 30 sec followed by a plate read. The primer concentration used for all qPCR reactions in this study was 0.5 µM. Standard curves were used in all qPCR analyses and were prepared by RT of equal volume of WT sample to other samples. After RT, five fivefold serial dilutions were made (six total standards including undiluted RT reaction) and added in duplicate to qPCR plates. After the qPCR run, samples were normalized to standard curve read using the Bio-Rad CFX Manager Software. See Supplemental Table S10 for all primer sequences used.

Antibodies

All antibodies used for this study are listed in Supplemental Table S11.

Western blot

To isolate protein for Western blotting, 0.8 × 10⁶ cells were washed with 1× PBS and then lysed with 500 µL RIPA buffer (10 mM Tris-Cl [pH 7.5], 1 mM EDTA, 0.5 mM EGTA, 1% NP40, 0.1% sodium-deoxycholate, 0.1% sodium dodecyl sulfate, 140 mM sodium chloride) supplemented with 1mM PMSF (Thermo Fisher #36978) and 1× Protease Inhibitor Cocktail (PIC; Sigma Product #P8340). Cell suspensions were rotated for 15 min at 4°C, then spun down at high speed at 4°C for 15 min and the supernatant was collected. Prior to western blotting, protein levels were quantified using the DC assay from Bio-Rad (product #5000006). 4× SDS loading buffer (Sigma Aldrich recipe: 0.2 M Tris-HCl pH 6.8, 0.4 M DTT, 8% [w/v] SDS, 6 mM Bromophenol blue, 4.3 M Glycerol) was added to samples to 1× final concentration. Samples were then boiled for 5 min at 95°C, and equal microgram amounts were loaded onto Bio-Rad TGX Stain Free Gels. Samples were run at 50 V until past stacking gel, then at 150 V for 1–2 h. Gels were transferred to PVDF (Immobulon #IPVH00010) membrane either for 1 h at 125 V at 4°C or overnight at 25 V at 4°C. Membranes were blocked for 45 min in 1× TBST + 5% milk. Membranes were then incubated with primary antibody either overnight at 4°C or for 1–3 h at RT. Membranes were washed 3× for 5 min each in 1× TBST. Secondary antibodies were diluted in 1× TBST + 5% milk and incubated with membranes for 45 min (1:100,000; Invitrogen). Membranes were then washed 3× in 1× TBST washes for 10 min, before being imaged with ECL (Thermo Fisher #34096). Antibodies used were FLAG (Sigma F1804, 1:1000), GAPDH (Abcam, ab9484, 1:1000), Total H3 (Proteintech 17168-1-AP, 1:1000), hnRNPU (Santa Cruz sc-32315 1:500), SAFB (Bethyl A300-812A, 1:3000), SAFB2 (Proteintech 11642-1-AP, 1:10,000), Goat anti-Mouse (Thermo Scientific A16072, 1:100,000), and Goat anti-Rabbit (Thermo Scientific G21234, 1:100,000).

Formaldehyde crosslinking of ESCs

For RIP and IP-MS, cells were grown to 75%–85% confluency, trypsinized and counted. Cells were washed twice in cold 1× PBS then rotated for 30 min in 10 mL of 0.3% formaldehyde (1 mL 16% methanol-free formaldehyde [Pierce #28906] in 49 mL 1× PBS) at 4°C. Formaldehyde was quenched with 1 mL of 2 M glycine for 5 min at room temperature. Cells were washed 3× in cold 1× PBS, then resuspended in 1× PBS at 10 × 10⁶ cells per mL and spun down. PBS was aspirated and pellets were snap frozen in a liquid nitrogen bath and immediately transferred to −80°C.

For ChIP, cells were grown to 75%–85% confluency and counted. Cells were washed once with 1× PBS and crosslinked with 0.6% formaldehyde for 10 min at room temperature. Formaldehyde was quenched with 557 µL of 2.5 M glycine and washed twice with cold 1× PBS. Cells were then scraped in 2 mL cold 1× PBS supplemented with 1× PIC (PIC; Sigma product #P8340). An amount of 10 mL 1× PBS + 0.05% tween-20 was added to collect the cells. Cells were spun down, resuspended in 1× PBS at 10 × 10⁶ cells per mL and spun down. PBS was aspirated and pellets were snap frozen in a liquid nitrogen bath and immediately transferred to −80°C.

RNA-IPs (RIPs)

RIPs were performed similar to Schertzer et al. (2019a), which is a protocol originally adapted from Hendrickson et al. (2016) and Raab et al. (2019). An amount of 25 µL protein A/G agarose beads (Santa Cruz sc-2003) were washed three times in blocking buffer (0.5% BSA in 1× PBS) and incubated overnight at 4°C with 10 µg antibody (anti-SAFB; Bethyl 812-300A; FLAG, Sigma F1804; V5, Sigma V8012; mouse IgG, Invitrogen 02-6502). A total of 10 × 10⁶ cells were resuspended in 500 µL RIPA Buffer (50 mM Tris-HCl, pH 8, 1% Triton X-100, 0.5% sodium-deoxycholate, 0.1% SDS, 5 mM EDTA, 150 mM KCl) supplemented with 1× Protease Inhibitor Cocktail (PIC; Sigma product #P8340), 2.5 µL SuperaseIN (Thermo Fisher Scientific AM2696) and 0.5 mM DTT and sonicated twice for 30 sec on and 1 min off at 30% output using the Sonics Vibracell Sonicator (model VCX130, serial #52223R). Samples were spun down at high speed and 50 µL total lysate was saved for input. Beads were washed three times in 1 mL fRIP buffer (25 mM Tris-HCl pH 7.5, 5 mM EDTA, 0.5% NP-40, 150 mM KCl) and resuspended in 450 µL fRIP buffer supplemented as above with PIC, SuperaseIN and DTT, then mixed with sonicated samples. Samples were rotated overnight at 4°C, then washed once with 1 mL fRIP buffer and resuspended in 1 mL PolII ChIP Buffer (50 mM Tris-HCl pH 7.5, 140 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 0.1% Sodium-deoxycholate, 0.1% Sodium dodecyl sulfide) before transferring to a new 1.7 mL tube. Samples were rotated at 4°C for 5 min, spun down at 1200g, and the supernatant aspirated. Samples were washed twice more with 1 mL PolII ChIP Buffer, once with 1 mL High Salt ChIP Buffer (50 mM Tris-HCl pH 7.5, 500 mM NaCl, 1 mM EDTA, 1 mM EGTA, 0.1% sodium-deoxycholate, 0.1% sodium dodecyl sulfide, 1% Triton X-100), and once in 1 mL LiCl buffer (20 mM Tris pH 8.0, 1 mM EDTA, 250 mM LiCl, 0.5% NP-40, 0.5% sodium-deoxycholate); each wash included a 5-min rotation at 4°C. At the final wash, samples were transferred to a new 1.7 mL tube. After the final wash, inputs were thawed on ice and bead samples were resuspended in 56 µL water, 33 µL of 3× reverse-crosslinking buffer (3× PBS, 6% N-lauroyl sarcosine and 30 mM EDTA), 5 µL 100 mM DTT, 20 µL Proteinase K, and 1 µL of SuperaseIN. Samples were incubated for 1 h at 42°C, then 1 h at 55°C, then 65°C for 30 min, and mixed by pipetting every 15 min. Afterward, 1 mL TRIzol was added, samples were vortexed, 200 µL CHCl₃ was added, samples were vortexed, and finally spun at 12,000g for 15 min at 4°C. The aqueous phase was then extracted and to that one volume of 100% ethanol was added. Samples were vortexed and applied to Zymo-Spin IC Columns (from #R1013) and spun for 30 sec at top speed on a benchtop microcentrifuge. An amount of 400 µL of RNA Wash Buffer (Zymo #R1013) was added and samples were spun at top speed for 30 sec. For each sample, 5 µL DNase I and 35 µL of DNA Digestion Buffer (Zymo #R1013) was added directly to the column matrix and incubated at room temp for 20 min. An amount of 400 µL of RNA Prep Buffer was then added (Zymo #R1013), and columns were spun at top speed for 30 sec. An amount of 700 µL RNA Wash Buffer (Zymo #R1013) was then added, and columns were spun at top speed for 30 sec. An amount of 400 µL RNA Wash Buffer was then added, and columns were spun at top speed for 30 sec. The flow through was discarded and columns spun again for 2 min to remove all traces of wash buffer. Columns were transferred to a clean 1.7 mL tube, 15 µL of ddH₂O was added to each column, and after a 5-min incubation, samples were spun at top speed to elute.

Total RNA isolation

ESCs were grown in six-well plates to ∼80% confluency. Cells were washed twice with 1× PBS and 1 mL of TRIzol was added per well. Samples were pipetted up and down at least 10 times, transferred to a microcentrifuge tube and briefly vortexed. Samples were incubated at RT for 5 min, then 200 µL of chloroform was added. Afterwards, samples were vortexed and incubated for 3 min at RT. Samples were spun down at 12,000g for 15 min at 4°C. The upper aqueous phase was moved to a new tube and 8 µL of linear acrylamide (Thermo Fisher, AM9520) was added. Then 500 µL of 100% isopropanol was added, and samples were vortexed and incubated at RT for 10 min. Tubes were spun down 12,000g for 10 min at 4°C. Supernatant was removed and pellets washed with 1 mL of cold 75% ethanol. Samples were briefly vortexed and spun down at 7500g for 5 min at 4°C. Supernatant was discarded and pellet was dried by repeated spin down and aspiration. Final pellets were resuspended in 100 µL water by pipetting up and down.

RNA sequencing

For RIP-seq inputs, 100 ng of RNA prepared from RIP input samples was used for library preparations. For RIP-seq RIP samples, 9 µL of RIP sample (from 15 µL total) were used. For total RNA-seq, 900 ng of total RNA was used. Each library preparation included 1 µL of 1:250 dilution of ERCC Spike-In RNAs (Ambion #4456653). An amount of 10 µL total was prepped using the KAPA RNA HyperPrep Kit with RiboErase (Kapa Biosystems; product #KR1351). Sequencing was performed on an Illumina Next-seq 500, using high-output, 75-cycle kits.

H3K27me3 and total H3 ChIP-seq

The day before sonication, 25 µL of protein A/G agarose beads (Santa Cruz sc-2003) were washed three times in block solution (0.5% BSA in 1× PBS) before being resuspended in 300 µL blocking solution. An amount of 10 µL per 10 million cells of antibody (Abcam mouse monoclonal ab6002) was added, then beads and antibody conjugated via rotation overnight at 4°C.

On the day of sonication, 10 million ESCs crosslinked with 0.6% formaldehyde were resuspended in lysis buffer 1 (50 mM HEPES pH 7.3, 140 mM NaCl, 1 mM EDTA, 10% glycerol, 0.5% NP-40, 0.25% Triton X-100, and 1× PIC [PIC; Sigma product #P8340]) incubated for 10 min at 4°C, and then incubated with lysis buffer 2 (10 mM Tris-HCl pH 8.0, 200 mM NaCl, 1 mM EDTA pH 8.0, 0.5 mM EGTA pH 8.0, and 1× PIC) for 10 min at room temperature. For H3K27me3 ChIPs, cells were resuspended in lysis buffer 3 (10 mM Tris-HCl pH 8.0, 100 mM NaCl, 1 mM EDTA pH 8.0, 0.5 mM EGTA pH 8.0, 0.1% Na-deoxycholate, 0.5% N-lauroyl sarcosine, and 1× PIC) and then sonicated.

ChIPs were then performed by incubating sonicated cell lysates at a concentration of 20 million cells/1 mL of lysis buffer 3 containing 1% Triton X-100 with preconjugated H3K27me3 antibody/agarose beads overnight at 4°C. After overnight H3K27me3 ChIP, beads were washed 5× in RIPA buffer (50 mM HEPES pH 7.3, 500 mM LiCl, 1 mM EDTA, 1% NP-40, and 0.7% Na-deoxycholate) for 5 min each and then once in TE. To elute the DNA, beads were resuspended in Elution buffer (50 mM Tris pH 8.0, 10 mM EDTA, and 1% SDS) and placed on a 65°C heat block for 17 min with frequent vortexing. Crosslinks were reversed overnight at 65°C, eluates were incubated with Proteinase K and RNase A, and DNA was extracted with phenol/chloroform and precipitated with ethanol. DNA was prepared for sequencing on the Illumina platform using Next Reagents (NEB) and Agencourt AMPure XP beads (Beckman Coulter).

IP-mass spectrometry sample preparation

An amount of 40 µL protein A/G agarose beads (Santa Cruz sc-2003) were washed three times in blocking buffer (0.5% BSA in 1× PBS) and incubated overnight at 4°C with 20 µg antibody (anti-FLAG; Sigma F1804). 30 × 10⁶ of ESCs crosslinked with formaldehyde as described above were resuspended in 500 µL RIPA Buffer (50 mM Tris-HCl, pH 8, 1% Triton X-100, 0.5% sodium-deoxycholate, 0.1% SDS, 5 mM EDTA, 150 mM KCl) supplemented with (1× PIC [PIC; Sigma product #P8340], 2.5 µL SuperaseIN [Thermo Fisher Scientific product #AM2696], and 0.5 mM DTT) and sonicated twice for 30 sec on and 1 min off at 30% output using Sonics Vibracell Sonicator (Model VCX130, serial #52223R). Samples were spun down at high speed and 50 µL total lysate was saved for input. Beads were washed three times in fRIP buffer (25 mM Tris-HCl pH 7.5, 5 mM EDTA, 0.5% NP-40, 150 mM KCl) and resuspended in 450 µL fRIP buffer supplemented with 1× PIC (PIC; Sigma product #P8340), 2.5 µL SuperaseIN (Thermo Fisher Scientific product #AM2696), and 0.5 mM DTT to bring samples to a 1:1 ratio of RIPA/fRIP buffer. Samples rotated overnight at 4°C. At 4°C, samples were then washed once with 1 mL fRIP buffer and then resuspended in 1 mL PolII ChIP Buffer (50 mM Tris-HCl pH 7.5, 140 mM NaCl, 1 mM EDTA, 1 mM EGTA, 1% Triton X-100, 0.1% Na-deoxycholate, 0.1% SDS) and transferred to a new Eppendorf tube. Samples were washed two more times with PolII ChIP Buffer, twice with High Salt CLIP Buffer (50 mM Tris-HCl pH 7.4, 1 M NaCl, 1 mM EDTA, 1% NP-40, 0.5% Na-deoxycholate, 0.1% SDS), and resuspended in 1 mL LiCl buffer (20 mM Tris pH 8.0, 1 mM EDTA, 250 mM LiCl, 0.5% NP-40, 0.5% Na-deoxycholate) and moved to a new microcentrifuge tube. Samples were then resuspended in cold 1× PBS and moved to a new Eppendorf tube. Samples were washed 3× with 1 mL cold 1× PBS. Twenty-five percent of samples were saved for western blot, and the remaining 75% were subjected to on-bead trypsin digestion as previously described (Rank et al. 2021). Briefly, after the last wash buffer step during affinity purification, beads were resuspended in 50 µL of 50 mM ammonium bicarbonate (pH 8). On-bead digestion was performed by adding 50 µL of 50 mM ammonium bicarbonate (pH 8) and 1 µg trypsin and incubated, shaking, overnight at 37°C. Beads were pelleted and supernatants were transferred to fresh tubes. The beads were washed twice with 100 µL LC–MS grade water, and washes added to the original supernatants. Samples were acidified by adding formic acid to a final concentration of 2%, to pH ∼ 2. Peptides were desalted using peptide desalting spin columns (Thermo), lyophilized, and stored at −80°C until further analysis.

LC/MS/MS analysis

The peptide samples were analyzed by LC/MS/MS using an Easy nLC 1200 coupled to a QExactive HF Biopharma mass spectrometer (Thermo Scientific). Samples were injected onto an Easy Spray PepMap C18 column (75 µm id × 25 cm, 2 µm particle size) (Thermo Scientific) and separated over a 2-h method. The gradient for separation consisted of 5%–45% mobile phase B at a 250 nL/min flow rate, where mobile phase A was 0.1% formic acid in water and mobile phase B consisted of 0.1% formic acid in acetonitrile (ACN). The QExactive HF was operated in a data-dependent mode where the 15 most intense precursors were selected for subsequent fragmentation. Resolution for the precursor scan (m/z 350–1700) was set to 60,000, while MS/MS scans resolution was set to 15,000. The normalized collision energy was set to 27% for HCD. Peptide match was set to preferred, and precursors with unknown charge or a charge state of 1 and ≥7 were excluded.

RIP- and RNA-seq alignment

RIP and RNA-seq data were aligned to the mm10 mouse genome using STAR with default parameters (Dobin et al. 2013). Alignments with a quality score ≥30 were retained for subsequent analyses (Li et al. 2009). For analysis of multimapping reads with TElocal (Jin et al. 2015), additional parameters (‐‐winAnchorMultimapNmax 100 and ‐‐outFilterMultimapNmax 100) were specified in STAR alignments.

RIP peak calling

After alignment and filtering, all SAFB RIP-seq data from WT ESC replicates were concatenated, and using SAMtools, were split into two files, corresponding to alignments that mapped to the positive and negative strands of the genome, respectively. Using a custom perl script, the strand information within the positive and negative strand alignment files was randomized so as to better match the criteria of the MACS peak caller, which uses the average distance between positive and negative strand alignments to estimate the fragment length (Zhang et al. 2008). Putative peaks were called on strand-randomized positive and negative strand alignment files, respectively, using default MACS parameters and not providing a background file (Zhang et al. 2008). Peak bed files were converted to SAF format and reads under each putative peak were counted from SAFB RIP-seq alignments performed in WT and DKO ESCs using featureCounts (Liao et al. 2014). We retained putative peaks that were represented by at least five reads in at least two of the three WT ESC lines profiled. We then used DESeq2 under default parameters to identify those putative peaks that were ascribed a P-value of <0.05 by DESeq2 when comparing signal between SAFB/2 DKO ESCs and WT ESCs. Lastly, we retained only those putative peaks that harbored an average aligned-reads-per-million-total-reads (RPM) signal of at least twofold less in SAFB/2 DKO ESCs compared to WT ESCs. This yielded 32,354 regions that were potentially enriched in their association with SAFB in WT ESCs. As a final filtering step, we used featureCounts to count the number of reads under these 32,354 regions that aligned to the mm10 genome with quality scores of ≥30 from within the SAFB-FLAG and GFP-FLAG RIP-seq data sets. From the initial set of 32,354 regions, 23,853 had a total of at least five reads distributed between the SAFB-FLAG and GFP-FLAG; of these, 1356 regions had higher signal in the GFP-FLAG compared to the SAFB-FLAG RIP-seq data set and were dropped from further analysis, yielding a total of 30,998 regions that we defined as SAFB-associated peaks (Supplemental Table S1).

RIP scatter plots

Scatter plots in Figure 2 were constructed using featureCounts to count the reads under each of the 30,998 SAFB peaks in each data set. Read counts were then plotted using R (R Core Team 2021).

SAFB motif analysis

To identify the motifs associated with SAFB peaks, we provided the sequences of the two thousand peaks with the greatest level of SAFB signal (top ∼10% of peaks) as input to the sensitive, thorough, rapid, enriched motif elicitation (STREME) tool from the MEME Suite (Bailey et al. 2015; Bailey 2021). Randomized control sequences with lengths the same as each of the peak sequences were developed with weighted nucleotide occurrence based on the mononucleotide content of the mm10 reference genome. The ‐‐rna flag was specified to account only for single-stranded analysis and motif width was restricted to between four and eight nucleotides; the motifs with the top three most significant P-values are shown in Figure 2D.

UCSC wiggle density plots

UCSC wiggle density plots were made from filtered sam files using custom perl scripts. Tracks of individual and pooled replicates are located here: https://genome.ucsc.edu/s/recherney/Cherney_Safb_2022.

Intersection of SAFB peaks and genic features

To identify the genic features under each SAFB peak, the GENCODE Basic vM25 GTF was downloaded and modified to include annotations of 5′- and 3′-UTRs (Frankish et al. 2021). From this file, features mapping to protein-coding genes (GENCODE transcript_type: “protein_coding”) and lncRNAs (GENCODE transcript_types: “bidirectional_promoter_lncRNA,” “macro_lncRNA,” “antisense,” “3prime_overlapping_ncRNA,” “lincRNA,” “processed_transcript,” “sense_intronic,” and “sense_overlapping”) were extracted and intersected with SAFB peaks using bedtools (Quinlan and Hall 2010). Peaks were classified as exon-overlapping if they fell within a gene and overlapped >50% of the exon in question, otherwise, they were classified as intron-overlapping. The classification of each peak can be found in Supplemental Table S1.

Intersection of SAFB peaks and repeat-masked elements

To determine whether SAFB peaks overlapped with repeat-masked genomic elements more than would be expected by random chance, repeat-masked elements were first extracted from all 20 mouse autosomes and the X chromosome using the UCSC genome browser MySQL relational database (Lee et al. 2022). Peaks were then intersected with repeat-masked elements using bedtools (Quinlan and Hall 2010). To estimate what level of intersection might be expected from random chance, we performed 1000 repetitions of the following process: the starting position of each peak in the set of SAFB peaks was shifted randomly to a new position between 2000 and 10,000 bases upstream or downstream; then, each complete set of randomized set of peaks was intersected with repeat-masked elements extracted from UCSC. The error bars in Figure 2G represent the standard deviation of the number of intersections with each class of repeat-masked element from each set of randomized peaks.

Intersection of multimapping RIP reads with repeat-derived elements

To determine the relative representation of repeat-derived elements in multimapping reads from SAFB RIPs, reads were aligned to mm10 with STAR using the parameters recommended by TElocal (‐‐winAnchorMultimapNmax 100 and ‐‐outFilterMultimapNmax 100; Jin et al. 2015). SAMtools was then used to extract alignments with MAPQ = 0 (i.e., multimapping reads; Li et al. 2009). Relative read abundance over repeat-derived elements was then calculated with TElocal, using the prebuilt “mm10_rmsk_TE.gtf.locInd” index and the “‐‐stranded reverse” option (Jin et al. 2015). Counts from TElocal were converted to RPM values (reads-per-million-total-reads). We retained only those elements that were represented by an RPM value of ≥1 summed across all data sets, were represented by an average RPM of ≥0.25 in the WT RIPs, and that had greater than or equal to twofold higher RPM in the WT compared to DKO RIPs. To determine the expected representation of each class of repeat in this final list of filtered elements, we summed the genomic space occupied by each class of repeat in the mm10_rmsk_TE.gtf from TElocal, and used this information to calculate the expected genomic space occupied by each class of repeat in our final list of filtered elements. We then used χ² tests to determine whether the actual genomic space occupied by each class of repeat in our final list of filtered elements differed significantly from what was expected. Classes of repeat that were represented by less than 25 elements in our final filtered list were not plotted in Figure 2H. Processed data used to generate Figure 2H are included in Supplemental Table S2.

Differential gene expression analyses

To detect genes that were differentially expressed between WT and SAFB/2 DKO ESCs, we performed RNA-seq on “Input” RNA extracted from the same sonicated extracts of formaldehyde-crosslinked WT and DKO ESCs samples that were used to perform SAFB RIP-seq described in Figure 2—the RNA extraction protocol is detailed in the “RNA-IPs” section of the methods. “Input” RNA Reads were aligned to mm10 using STAR and default parameters (Dobin et al. 2013), alignments were filtered to retain only those reads with a mapping-quality of ≥30 using SAMtools (Li et al. 2009), and then the number of filtered reads mapping to each GENCODE vM25 gene was counted using featureCounts (Liao et al. 2014): [-g gene_name -s 2 -a gencode.vM25.basic.annotation.gtf -o]. Genes that had less than 10 total reads summed across all five samples (three WT and two DKO) were excluded from downstream analyses. Counts were loaded into DESeq2, and the genes that had adjusted P-values for differential expression between WT and DKO samples of <0.05 were retained and reported as “significant” in Supplemental Table S3.

Differential splicing analyses

Fastq files were aligned to the mm10 genome using STAR v2.7.10b using default parameters (Dobin et al. 2013). Files were sorted and indexed using SAMtools v1.16 (Li et al. 2009) before using rMATS 4.1.1 to compare splicing events in WT and DKO ESCs, and in GFP- and SAFB-rescue DKO ESCs (Shen et al. 2014). Events were called significant if their PSI values in WT versus DKO ESCs, or GFP- versus SAFB-rescue DKO ESCs, were ascribed an adjusted P-value of <0.05. Significant events can be found in Supplemental Table S6.

Gene set enrichment analyses

Gene set enrichment analyses of differentially expressed genes were performed using the MSigDB webserver (Liberzon et al. 2015). The lists of significantly down- and up-regulated genes were input, orthology mapped onto the human genome, and queried for overlap with the Hallmark, Chemical and Genetic Perturbations, and GO Biological Process Gene Sets. The top 20 most significantly overlapping gene sets from each search are reported in Supplemental Table S4.

Assessing SAFB RIP signal over spliced and unspliced transcripts

To determine the extent to which SAFB was enriched over expressed spliced and unspliced transcripts, we took advantage of the kallisto algorithm, which was designed to enable probabilistic alignment of short-read RNA-seq data (Bray et al. 2016). To enable the detection of unspliced transcripts, we created a version of the GENCODE vM25 basic transcriptome that for each gene, included one representative unspliced transcript that began at the first annotated transcription start and the last annotated transcription end (vM25_basic_complete.fa).

In parallel, because like all RIP- or CLIP-seq data sets, our RIP-seq data sets contained reads that align to genomic regions that were not classified as peaks, we selected for our downstream analyses only the subset of RIP-seq reads that aligned under SAFB peaks. To do this, SAFB RIP-seq reads from WT and DKO data sets were aligned to mm10 using STAR and filtered for mapping quality ≥30 using SAMtools (Li et al. 2009; Dobin et al. 2013). We then used SAMtools to split the alignments by strand, and for each stranded alignment file, again using SAMtools, selected the subset of alignments from the WT and DKO data sets that aligned under each SAFB peak. Still using SAMtools, we converted the bam alignments back into fastq data. These final fastq files represent the subset of RIP-seq data that aligned under each classified peak and exclude most noise in the WT data set. Input RNA-seq and the subset RIP-seq data were then aligned with kallisto to an index made from vM25_basic_complete.fa, using the options [-l 200 -s 50 ‐‐rf-stranded]. For each transcript isoform in vM25_basic_complete.fa, TPM counts reported from the DKO RIP-seq data set were subtracted from TPM counts in the WT RIP-seq data set. The output file then underwent the following filtering parameters: (i) Transcript isoforms that we had previously filtered out prior to performing DESeq analyses were excluded; (ii) for each gene remaining, we retained the single spliced isoform with the highest expression level as representative; and (iii) transcript isoforms whose expression in the WT input total RNA-seq data were less than 0.125 TPM were excluded, including 25 and 14 genes that were originally called significantly up- and down-regulated by DESeq2, respectively. Finally, we split the transcripts of the genes that did not significantly change in expression upon DKO into three categories: those with low (0.125–1 TPM), medium (>1 but <16 TPM) and high (>16 TPM) levels of expression. Data used to make plots in Figure 3D are included in Supplemental Table S5.

IUPred2

Protein disorder plots for SAFB were constructed using the webserver for IUPred2A (Erdos and Dosztanyi 2020).

Kallisto analyses of rescue RNA-seq data

To determine the relative abundance of spliced and unspliced transcript isoforms between the SAFB-FL-WT, SAFB-ΔDD3, and GFP rescue data sets, RNA-seq data were aligned using kallisto to an index made from vM25_basic_complete.fa, using the options [-l 200 -s 50 ‐‐rf-stranded]. The output file then underwent the same filtering as in Figure 3D: (i) Transcripts isoforms that we had previously filtered out prior to performing DESeq analyses were excluded (62,934), (ii) for each gene remaining, we retained the single spliced isoform with the highest expression level as representative, and lastly, (iii) transcript isoforms whose averaged expression of the two SAFB-FL-WT replicates in the total RNA-seq data were less than 0.125 TPM were excluded, including 23 and 20 genes that were originally called significantly up- and down-regulated by DESeq2, respectively. We then normalized each individual replicate TPM value to the SAFB-FL-WT TPM average. P-values were determined using a paired t-test. Data used to make plots in Figure 5D and E are included in Supplemental Table S5.

Mass spectrometry data analysis

Raw data files were processed using MaxQuant version 1.6.15.0 and searched against the reviewed mouse database (containing 17,051 entries), appended with a contaminants database, using Andromeda within MaxQuant. Enzyme specificity was set to trypsin, up to two missed cleavage sites were allowed, and methionine oxidation and amino-terminus acetylation were set as variable modifications. A 1% FDR was used to filter all data. Match between runs was enabled (5 min match time window, 20 min alignment window), and a minimum of two unique peptides was required for label-free quantitation using the LFQ intensities. Perseus was used for further processing (Tyanova et al. 2016). Only proteins with >1 unique + razor peptide were used for LFQ analysis. Proteins with 50% missing values were removed and missing values were imputed from normal distribution within Perseus. Log₂ fold change (FC) ratios were calculated using the averaged log₂ LFQ intensities of the IP sample compared to the GFP control. Proteins with log₂ FC > 1 were considered biological interactors and analyzed further. All analyzed protein interaction data are present in Supplemental Table S7. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository (Perez-Riverol et al. 2022) with the data set identifier PXD038103.

DAVID GO analysis

GO analyses were conducted using DAVID (Huang da et al. 2009; Sherman et al. 2022). Genes were searched against UP_ KW_BIOLOGICAL_PROCESS, UP_KW_CELLULAR_COMPONENT, UP_KW_MOLECULAR_FUNCTION, GOTERM_BP_DIRECT, GOTERM_CC_DIRECT and GOTERM_MF_DIRECT. The list shown in Figure 6 represents the union of the top twelve most enriched CC and MF GO terms from the WT SAFB and HNRNPU pulldowns that also passed an FDR of <0.01. The bubbleplot in Figure 6 was created using tidyverse v 1.3.1 package in R version 4.0.4 (Wickham et al. 2019; R Core Team 2021).

Custom gene set enrichment analyses

To determine whether the SAFB and HNRNPU immunoprecipitates were enriched in proteins found in nuclear speckles (Saitoh et al. 2004), paraspeckles (Yamazaki and Hirose 2015), or proteins that harbor RS domains (Cascarina and Ross 2022), we used the lists of proteins reported in the aforementioned references and the 69 and 165 proteins that we classified as enriched over GFP control in the SAFB and HNRNPU IPs, respectively (Supplemental Table S7). We followed the gene ranking metric referred to as “log₂ ratio of classes” (LRC) and the GSEA framework described in Subramanian et al. (2005). However, we decided to use a custom version of GSEA, instead of the standard version, to account for (i) the limited number of genes in our two data sets, (ii) the number of replicates available (two for each sample), which is lower than the canonical threshold of at least seven recommended by the GSEA authors for phenotype permutations, and (iii) the few gene sets (only three), which were tested at the same time. Briefly, prior to calculating enrichments, the average LFQ values per protein per data set were calculated and divided by the average LFQ values of each protein in the GFP IP; this ratio was then log₂ transformed. Enrichment and depletion scores were then calculated separately for each data set. To calculate the enrichment/depletion score (EDS), for each gene set and data set of interest, we first converted log₂-transformed ratios into a ranked list. The highest rank was defined as the numerical value that corresponds to the total number of rows in the list in question and was assigned to the corresponding gene in the data set that had the highest log₂-transformed ratio in the list. The lowest rank was defined as a value of 1 and was assigned to the gene in the data set that had the lowest log₂-transformed ratio in the list. These ranks were then assigned to the genes of each gene set and their averages became the EDSs, which are specific for each data set (SAFB and HNRNPU), gene set (speckle, paraspeckle, and SR proteins), and condition (full protein vs. protein with a deleted domain). The neutral point (NP) for each data set in each gene set is equal to the [(# genes in the data set +1)/2]. Specifically, the EDS of each gene set was then defined as the average rank of genes in each data set that were present in the gene set. Gene sets whose EDS > NP were classified as “enriched” and those whose EDS < NP were classified as “depleted.” To assess statistical significance, we generated random EDS values by averaging the ranks of as many randomly selected data set genes as those present in each gene set and repeated this process, for each data set, each gene set and each protein form (full or with a deleted domain) 100,000 times. Then, we assessed the probability that each EDS was produced by chance following the approach outlined in Mielke and Berry (2007), in which the P-value for enrichment ≈ [(number of permutated cases with EDS ≥ NP)/(number of total permutations performed)] or, in the case of depletion, using as the numerator of this ratio [# permutated gene sets with EDS ≤ NP]. Resulting P-values were Bonferroni-corrected, thus controlling for the FWER as recommended in Olejnik et al. (1997) and Subramanian et al. (2005). The FWER was assessed at three levels: 0.10 (*), 0.05 (**), and 0.01 (***). We performed Bonferroni correction by keeping together enriched and depleted gene sets, when they were present at the same time (namely, in the HNRNPU data set of the full protein), thus producing more conservative statistical results than performing the correction after splitting enrichment and depletion cases, as done in standard GSEA.

Analysis of H3K27me3 ChIP-seq data

H3K27me3 data were merged by genotype (WT and DKO, respectively) and total H3 data were merged. Data were then aligned to mm10 using bowtie2 (Langmead and Salzberg 2012). Peaks were called using MACS2 with total H3 as a control under the following parameters: [macs2 callpeak -t -c bam -n -f BAM -g mm –broad –broad-cutoff 0.01] (Zhang et al. 2008). H3K27me3 peak locations are included in Supplemental Table S7. WT and DKO peak files were catenated and then merged using bedtools [bedtools merge -I in.file > out.file]. The independent WT and DKO data sets were intersected with the merged data file [bedtools intersect -wao -header -a union file -b wt or dko file > outfile] to identify WT and DKO-specific H3K27me3 peaks.

H3K27me3 analyses

To annotate gene promoters, we took the gencode.vM25.basic.annotation.gtf file and filtered features for transcripts only. We then took the coordinates of the transcription start site from each transcript and added 2 kb upstream and 1 kb downstream (3 kb total) and used this region as the transcript promotor region. We then used featureCounts (featureCounts -s 0 -F SAF -a promoter.file -o out.file) to align H3K27me3 reads to our annotated promoters. We then used bedtools (bedtools intersect -wao -header -a promoter.file -b k27peaksfile > out.file) to intersect our previously called H3K27me3 peaks in our WT and DKO files with our promoter file to find the levels of H3K27me3 at promoters.

To perform allelic H3K27me3 analyses, variant sequence data for the mm10 genome build was obtained from the Sanger Institute (http://www.sanger.ac.uk/resources/mouse/genomes/; Keane et al. 2011), and a CAST/EiJ (CAST) pseudogenome was created as in Calabrese et al. (2012, 2015). Reads were aligned to both the B6 and CAST version of mm10 using STAR, and those that had a mapping quality greater than or equal to 30 were extracted with SAMtools. Reads that overlapped B6 or CAST SNPs were detected using a custom perl script as in Calabrese et al. (2012, 2015). Allelic tiling density plots over chr6 were created as in Schertzer et al. (2019a), using a bin size of 4000 bp.

Allelic RNA-seq analysis

To determine the extent of gene silencing by Xist across genotypes, RNA-seq reads were aligned and processed as described in Trotman et al. (2020, 2023). Briefly, reads were aligned to B6 and CAST pseudogenomes using STAR and default parameters (Dobin et al. 2013), alignments were filtered to retain only those reads with a mapping-quality of ≥30 using SAMtools (Li et al. 2009), and uniquely aligning reads that overlapped B6 or CAST SNPs were summed under each gene using custom perl scripts as in Trotman et al. (2020, 2023). Read counts per gene were then normalized to reads-per-million total reads (RPM). Across SAFB genotypes, we examined the set of 242 genes whose B6 expression levels on chr6 differed between Xist-expressing and nonexpressing ESCs at an adjusted P-value threshold of <0.01 (Supplemental Table S9; Trotman et al. 2023).