Cell and animal cultures, sample preparation and crosslinking
S. arctica coenocytic culture was grown in marine broth (Difco, 3704 g l−1) at 12 °C in 25 cm2 flasks. Cells were passaged every 7 days using a 1:100 dilution. To synchronize cells in the G1/early S phase, an 8-day old culture was treated with 200 mM hydroxyurea (Sigma-Aldrich, catalogue no. H8627) for 18 h in the presence of 0.3% dimethylsulfoxide (DMSO). Synchronized cells were pelleted at 2,000g for 5 min at 12 °C, washed twice with Ca2+/Mg2+-free artificial sea water (CMFSW) (10 mM HEPES (pH 7.4), 450 mM NaCl, 9 mM KCl, 33 mM Na2SO4, 2.5 mM NaHCO3) and flash-frozen in liquid nitrogen. Frozen cells were then reconstituted in CMFSW and crosslinked with 1% formaldehyde (Thermo Scientific, catalogue no. 28906) for 10 min under vacuum. The crosslinking reaction was quenched with 128 mM glycine for 5 min in the vacuum desiccator, followed by a 15 min incubation on ice. Cells were pelleted at 4 °C for 10 min at 2,000g, washed once with CMFSW, reconstituted in CMFSW to the concentration of 2 M ml−1 and crosslinked with 3 mM DSG (Thermo Scientific, catalogue no. A35392) for 40 min at room temperature on a rotating wheel. The reaction was quenched with 400 mM glycine for 5 min. Double-crosslinked cells were pelleted at 4 °C for 15 min at 2,000g and flash-frozen in liquid nitrogen.
C. owczarzaki strain ATCC30864 was maintained in axenic culture at 23 °C in the ATCC (American Type Culture Collection) medium 1034 (modified PYNFH medium) in 25 cm2 flasks. For subculture, filopodial cells were passaged every 2–3 days using a dilution of 1:100. Before collection, filopodial cells were synchronized in G1 or early S phase by treating a filopodial culture of 70–80% confluency with 100 mM hydroxyurea (Sigma-Aldrich, catalogue no. H8627) for 18 h (ref. 61). Synchronized cells were scraped off the surface and pelleted at 2,200g for 5 min at room temperature. Collected cells were crosslinked as described in ref. 62. Briefly, cells were crosslinked with 1% formaldehyde (Thermo Scientific, catalogue no. 28906) in PBS for 10 min on a rotating wheel at room temperature. The crosslinking reaction was quenched with 128 mM glycine for 5 min at room temperature followed by extra incubation on ice for 15 min. The crosslinked cells were pelleted at 4 °C for 10 min at 2,000g and washed once with ice-cold PBS. Cells were diluted in PBS to the concentration of 2 M ml−1 and also crosslinked with 3 mM DSG (Thermo Scientific, catalogue no. A35392) for 40 min at room temperature on a rotating wheel. The crosslinking was quenched with 400 mM glycine for 5 min. Cells were pelleted at 4 °C for 15 min at 2,000g and flash-frozen in aliquots of 2 million cells.
S. rosetta was cocultured with Echinicola pacifica bacteria in artificial sea water supplemented with 20% cereal grass media (CGM3) at 23 °C. To synchronize the cell culture in the G1 or early S phase, cells from a 3-day-old culture were pelleted at 2,000g for 10 min and diluted in 4% CGM3 in artificial sea water to the concentration of 300,000 cells per ml. Cells were treated with 0.05 mM aphidicolin (Sigma-Aldrich, catalogue no. 178273) in the presence of 0.3% DMSO. After 18 h of incubation, cells, including chain colonies, fast and slow swimmers, were pelleted at 2,000g for 15 min. To remove bacteria from the choanoflagellate culture, collected cells reconstituted in 1 ml of culture media were passed through a Ficoll layer (1.6% Ficoll (Sigma-Aldrich, catalogue no. F5415), 0.5 M sorbitol, 50 mM Tris-HCl (pH 8.8), 15 mM MgCl2, 1% artificial sea water) by centrifugation at 1,000g for 10 min at 4 °C. Pelleted choanoflagellate cells were then double-crosslinked with 1% formaldehyde in CMFSW and 3 mM DSG in CMFSW as described above for C. owczarzaki. The crosslinked cells were pelleted at 4 °C for 15 min at 2,000g and flash-frozen in liquid nitrogen.
E. muelleri sponges gemmules were hatched and grown for 1 week in Strekal’s media63 in 150 × 25 mm culture dishes (Corning, catalogue no. 353025). To isolate phagocytic choanocyte cell population, specimens were fed for 10 min with 0.5 µm fluorescent carboxylate-modified FluoSpheres (Invitrogen, catalogue no. F8813) added to Strekal’s media to final 0.02% concentration (1:100 dilution of stock 2% FluoSpheres slurry)64. Sponges were washed once with Strekal’s media, and 1% formaldehyde solution in Strekal’s media was added to crosslink specimens for 10 min at room temperature with occasional mixing. To quench formaldehyde, 128 mM glycine was added and incubated for 5 min at room temperature and 15 min on ice. Crosslinked sponge specimens were washed twice with ice-cold Strekal’s media. Roughly 80 specimens were transferred in 5 ml of the Strekal’s media and dissociated by trituration until all tissue was removed from the gemmule husks (roughly ten trituration passages). The dissociated cell suspension was filtered through a 40-µm cell strainer, and cells were diluted to 2 M ml−1 concentration. The second crosslinking was performed with 3 mM DSG (Thermo Scientific, catalogue no. 20593) in Strekal’s media for 40 min at room temperature on a rotating wheel. The reaction was quenched with 400 mM glycine for 5 min at room temperature. Crosslinked cells were pelleted at 4 °C for 15 min at 2,000g, and then resuspended in 2 ml of ice-cold Strekal’s media with 2 µg ml−1 Hoechst 33342 (Thermo Scientific, catalogue no. 62249). Choanocytes were isolated using a BD FACS Aria II sorter with BD FACSDiva v.6.1.3 (BD Biosciences) as cells showing both FluoSphere fluorescence and Hoechst nuclei staining. Fluorescence-activated cell sorting (FACS) profiles were analysed with FlowJo v.10.7 (Extended Data Fig. 1b).
M. leidyi specimens were kept in 300-ml glass beakers with 5–10 individuals at 21 °C in artificial sea water (Red Sea, catalogue no. R11055) with a salinity of 27 ppt. Ctenophores were fed daily with a mixture of living rotifers (Brachionus sp.) and brine shrimps (Artemia salina). The water was exchanged once a week. For all experiments, adult lobate animals were starved for 2 days before collection. To dissociate animal tissue, roughly five adult animals (10 mm long) were transferred into CMFSW and washed twice to exchange the buffer. Animal tissue was dissociated into single cells in 5 ml of fresh CMFSW by triturating every 2 min for a total of 10 min. The efficiency of tissue dissociation was monitored under the microscope. Dissociated cells were filtered through a 40-µm cell strainer and diluted to 2 M ml−1 for the subsequent formaldehyde crosslinking. Cells were crosslinked in 1% formaldehyde in CMFSW for 10 min at room temperature. The reaction was stopped with 128 mM glycine for 5 min at room temperature and 15 min on ice. Crosslinked cells were pelleted at 4 °C for 10 min at 2,000g, washed once with CMFSW and resuspended to 2 M ml−1 for a second crosslinking with 3 mM DSG in CMFSW. The crosslinking reaction was stopped after 40 min of incubation at room temperature on a rotating wheel with 400 mM glycine for 5 min. The crosslinked cells were pelleted at 4 °C for 15 min at 2,000g.
H. californensis specimens from the first generation (F1) of a laboratory-reared culture at the Monterey Bay Aquarium (USA) were flash-frozen and pulverized in liquid nitrogen. Extracted cells and nuclei were filtered through a 40-µm cell strainer and pelleted by centrifugation at 4 °C for 10 min at 2,000g. Cells were double-crosslinked with 1% formaldehyde in CMFSW and 3 mM DSG in CMFSW as described for M. leidyi.
T. adhaerens and C. collaboinventa colonies were grown in 200 × 30 mm glass Petri dishes at 21 °C in artificial sea water (Red Sea, catalogue no. R11055) with a salinity of 33 ppt. Placozoans were fed once a week with unicellular algae (Pyrenomonas sp.), the water was exchanged every second week. To prepare single-cell suspension, roughly 500 animals were collected, washed twice with CMFSW and resuspended in 1 ml of CMFSW supplemented with 2 mM EDTA. Animal tissue was triturated every 2 min for a total of 10 min at room temperature. The efficiency of dissociation was monitored under the microscope. Dissociated cells were filtered through a 40-µm cell strainer, diluted to 2 M ml−1 and crosslinked with 1% formaldehyde in CMFSW for 10 min at room temperature on a rotating wheel. The reaction was quenched with 128 mM glycine for 5 min at room temperature and 15 min on ice. Cells were pelleted at 4 °C for 10 min at 2,000g, washed once with CMFSW and resuspended in 3 mM DSG in CMFSW for a second crosslinking. After 40 min of incubation at room temperature on a rotating wheel, 400 mM glycine was added to stop the reaction and cells were pelleted at 4 °C for 15 min at 2,000g.
N. vectensis NvElav1::mOrange transgenic line65 was maintained in one-third artificial sea water (Red Sea, catalogue no. R11055) with salinity of 14 ppt. To isolate NvElav1::mOrange positive cells, 1.5–2-month-old animals starved for 1 day before the experiment were crosslinked with 1% formaldehyde in Ca2+/Mg2+-free one-third sea water (one-third CMF: 17 mM HEPES (pH 7.4), 167 mM NaCl, 9 mM NaHCO3, 3.3 mM KCl) for 10 min under vacuum. The crosslinking reaction was stopped by adding 128 mM glycine and incubating the tissue under vacuum for 5 min, followed by a 15 min incubation on ice. The crosslinked tissue was dissociated into single cells by incubating the tissue with 10 mg ml−1 of Protease XIV (Sigma-Aldrich, catalogue no. P5147) in one-third CMF and 1 mM CaCl2 for 5 min at 24 °C triturating the tissue every 1 min. The digested tissue was pelleted at 800g for 5 min, reconstituted in one-third CMF supplemented with 2 mM EDTA and 2 µg ml−1 Hoechst 33342 (Thermo Scientific, catalogue no. 62249), and the trituration continued for another 5–10 min. Dissociated cells were filtered through a 40-µm cell strainer, and neurons were isolated using a BD FACS Aria II as cells showing both the mOrange signal and Hoechst nuclei staining (Extended Data Fig. 1b). Isolated NvElav1::mOrange positive cells were also crosslinked with 3 mM DSG for 40 min at room temperature.
Micro-C library preparation
Micro-C libraries were prepared as previously described11,12 with the following modification. Double-crosslinked cells (2 million cells per sample) with 1% formaldehyde and 3 mM DSG were permeabilized with 500 µl of MB1 buffer (10 mM Tris-HCl (pH 7.4), 50 mM NaCl, 5 mM MgCl2, 1 mM CaCl2, 0.2% NP-40, protease inhibitor cocktail) for 20 min on ice with occasional trituration. Cells were pelleted at 4,500g for 5 min at 4 °C and washed once with MB1 buffer. To digest chromatin to a 80% monomers to 20% dimer and oligomers nucleosome ratio, an appropriate amount of MNase (Takara Bio, catalogue no. 2910a) was added (Extended Data Fig. 1a), and samples were incubated for 10 min at 37 °C with mixing at 850 rpm. The digestion reaction was stopped with 4 mM EGTA (pH 8.0) followed by incubation at 65 °C for 10 min without agitation. Cells were washed twice with ice-cold MB2 buffer (10 mM Tris-HCl (pH 7.4), 50 mM NaCl, 10 mM MgCl2, 0.1% BSA) and pelleted at 4,500g for 5 min at 4 °C. Next, to repair the fragment ends after MNase digestion, pelleted cells were resuspended in the repair reaction mix (5 µl of 10× NEBuffer 2.1, 34 µl of nuclease-free water, 1 µl of 100 mM ATP, 2.5 µl of 100 mM DTT) supplemented with 2.5 µl of 10 U µl−1 T4 PNK (NEB, catalogue no. M0201). After 15 min of incubation at 37 °C with 850 rpm agitation, 5 µl of 5 U µl−1 Klenow Fragment (NEB, catalogue no. M0210) was added to generate 3′–5′ overhangs in the absence of dNTPs for a subsequent incorporation of biotin-labelled dNTPs. The reaction mixture was incubated for another 15 min at 37 °C at 850 rpm. To biotinylate DNA fragment ends, the mixture of dNTPs was added to the reaction mix (2.5 µl of 10× T4 DNA Ligase buffer, 11.875 µl of nuclease-free water, 5 µl of 1 mM Biotin-dATP (Jena Bioscience, catalogue no. NU-835-BIO14), 5 µl of 1 mM Biotin-dCTP (Jena Bioscience, catalogue no. NU-809-BIOX), 0.5 µl of a mixture of 10 mM dTTP and dGTP, 0.125 µl of 20 mg ml−1 BSA). After 45 min of incubation at room temperature with interval mixing at 850 rpm, the reaction was stopped with 30 mM EDTA (pH 8.0) followed by incubation at 65 °C for 20 min without agitation. The chromatin from lysed cells and nuclei was pelleted at 10,000g for 10 min at 4 °C and washed twice with MB3 buffer (50 mM Tris-HCl (pH 7.5), 10 mM MgCl2). Finally, the chromatin was resuspended in 1,200 µl of proximity ligation mix (920 µl of nuclease-free water, 120 µl of 10× T4 DNA Ligase buffer, 100 µl of 10% Triton X-100, 12 µl of 20 mg ml−1 BSA, 36 µl of 50% PEG 4000, 12 µl of 5 U µl−1 T4 DNA ligase (Thermo Scientific, catalogue no. EL0012)) and incubated at room temperature for at least 2.5 h. To remove biotin from unligated ends, pelleted chromatin was treated with 2 µl of 100 U µl−1 Exonuclease III (NEB, catalogue no. M0206) for 5 min at 37 °C and agitation 850 rpm. Then, chromatin was decrosslinked and deproteinased overnight at 65 °C at 850 rpm in the presence of 350 mM NaCl, 1% SDS and 1 mg ml−1 proteinase K (Roche, catalogue no. 3115879001). The DNA was purified using DNA Clean & Concentrator-5 kit (Zymo Research, catalogue no. D4014) and eluted in 50 µl of 10 mM Tris-HCl (pH 8.0) (Extended Data Fig. 1c). Next, biotinylated proximity ligated DNA fragments were captured with Dynabeads MyOne Streptavidin (Life Technologies, catalogue no. 65602). DNA ends were prepared for adapter ligation and dA-tailed using NEBNext End repair/dA-tailing mix (NEB, catalogue no. E7546). The Y-shaped Illumina adapters were ligated with NEBNext Ultra II Ligation Module (NEB, catalogue no. E7595S), and the final library was amplified using NEBNext High-Fidelity 2× PCR Master Mix (NEB, catalogue no. M0541). The final libraries were double-size selected with Ampure XP (Beckman Coulter, catalogue no. A63881) resulting in libraries ranging from 350 to 750 bp in length. The detailed Micro-C stepwise protocol is reported in Supplementary Text 1.
High molecular weight gDNA extraction for genome sequencing
Genomic DNA (gDNA) from C. owczarzaki (Cowc) strain ATCC30864 was extracted with Blood & Cell Culture DNA Mini Kit (Qiagen, catalogue no. 13323). The library was constructed by the use of Ligation Sequencing Kit (Oxford Nanopore, catalogue no. SQK-LSK109) and NEBNext Companion Module (NEB, catalogue no. E7180), and sequenced with the R9.4.1 Flow Cell set on a MinION device (Oxford Nanopore). We obtained 4.3 M reads with an estimated Oxford Nanopore N50 of 5.4 kb.
E. muelleri gDNA was isolated using the Nanobind Tissue (Circulomics, catalogue no. NB-900-701-01) from 177 mg of frozen tissue of clonal juvenile sponges hatched from overwintering cysts (gemmules). Gemmules were obtained from the head tank of the Kapoor Tunnel (Sooke Reservoir), part of the drinking water system of the city of Victoria, British Columbia, Canada21. Short DNA fragments of less than 10 kb were removed with Short Read Eliminator Kit (Circulomics, catalogue no. SS-100-101-01). gDNA was quantified with a Qubit fluorometer and sequenced on an Oxford Nanopore using a PromethION flow cell (R9.4), producing 5.31 million reads with an estimated Oxford Nanopore N50 of 18.97 kb.
To reduce the level of heterozygosity during the assembly of M. leidyi genome (below), an animal culture was established from a single individual through self-fertilization. High molecular weight DNA was isolated from 5–8 animals (3–5 cm) starved for 24 h before flash-freezing. Frozen tissues were powdered with mortar and pestle, dissolved in 10 ml of urea extraction buffer (50 mM Tris-HCl (pH 8.0), 7 M Urea, 312.5 mM NaCl, 20 mM EDTA (pH 8.0), 1% w/v N-lauroylsarcosine sodium salt) as described in ref. 66 and incubated for 10 min at room temperature on a rocking platform 20 rpm. gDNA was then purified twice with a phenol-chloroform-isoamyl alcohol mixture pH 7.7–8.3 (Sigma-Aldrich, catalogue no. 77617), precipitated with 0.7 volume of 100% isopropanol and subsequently washed twice with 70% ethanol. Finally, the isolated DNA was subjected to another round of purification with Nanobind Tissue kit (Circulomics, catalogue no. NB-900-701-01), followed by short-read elimination with the Short Read Eliminator Kit (Circulomics, catalogue no. SS-100-101-01). Sequencing was performed on Oxford Nanopore using PromethION flow cell (R9.4). We obtained 4.54 million reads with an estimated Oxford Nanopore N50 of 36.84 kb.
ATAC-seq library preparation
ATAC-seq libraries from M. leidyi and from sorted choanocytes of E. muelleri were prepared using Omni-ATAC protocol as described previously67. Briefly, two M. leidyi adult specimens were dissociated using CMFSW with 0.25% α-Chymotrypsin (Sigma-Aldrich, catalogue no. C8946). To isolate nuclei, dissociated cells were transferred into cold hypotonic ATAC lysis buffer adjusted for marine animals (10 mM Tris-HCl (pH 7.5), 35 mM NaCl, 3 mM MgCl2, 0.1% Tween-20, 0.01% NP-40, 0.01% digitonin, 70 µM Pitstop (Abcam, AB1206875MG)). Cell lysis was stopped after 2 min by adding marine ATAC wash buffer (10 mM Tris-HCl (pH 7.5), 35 mM NaCl, 3 mM MgCl2, 0.1% Tween-20, 1% BSA). Nuclei were then pelleted and resuspended in cold PBS buffer with 0.8 M Sorbitol. We used 50,000 nuclei per each tagmentation reaction.
To sort choanocytes of E. muelleri, 7 days posthatching sponges were fed with 0.5 µm fluorescent carboxylate-modified FluoSpheres (Invitrogen, catalogue no. F8813). After 10 min of incubation, sponges were washed twice with Strekal’s media, collected and dissociated for 15 min at 28 °C using Protease XIV (Sigma-Aldrich, catalogue no. P5147) in Strekal’s media and 1 mM CaCl2. Cells were pelleted at 800g for 5 min at room temperature and resuspended in Strekal’s media with 2 mM EDTA (pH 8.0). Further dissociation and trituration of sponge tissue continued for another 15 min at room temperature. Cells were filtered through a 40-µm cell strainer, stained with 2 µg ml−1 Hoechst 33342 and sorted using FACS. Sorted cells were lysed for 3 min in ATAC lysis buffer (10 mM Tris-HCl (pH 7.5), 10 mM NaCl, 3 mM MgCl2, 0.1% NP-40, 0.1% Tween-40, 0.01% digitonin). For each tagmentation reaction we used 100,000 nuclei. ATAC-seq libraries were prepared as described previously67 and sequenced on Illumina NextSeq 500 using High-Output 75 cycles.
iChIP–seq library preparation
For S. arctica, C. owczarzaki, S. rosetta, M. leidyi, E. muellleri, T. adhaerens, C. collaboinventa and N. vectensis, double-crosslinked cells, as above, were washed with PBS, resuspended in 500 µl of cell lysis buffer (20 mM HEPES (pH 7.5), 10 mM NaCl, 0.2% IGEPAL CA-630, 5 mM EDTA, protease inhibitors cocktail) and incubated on ice for 10 min. Samples were centrifuged at 16,000g for 10 min at 4 °C. The resulting pellets were resuspended in bead beating buffer (20 mM HEPES (pH 7.5), 10 mM NaCl, 5 mM EDTA, protease inhibitors cocktail), and then transferred to 0.2 ml tubes containing acid-washed glass beads (Sigma-Aldrich, G8772). Cells were lysed by vortexing five times for 30 s. The supernatant was transferred to a 1.5-ml sonication tube, SDS was added to 0.6% and samples were sonicated 3–5 cycles of 30 s on, 30 s off in a Bioruptor Pico (Diagenode) to generate 200–300 bp fragments. Chromatin was diluted with 5 volumes of dilution buffer (20 mM HEPES (pH 7.5), 140 mM NaCl), centrifuged at 16,000g for 10 min at 4 °C and stored at −80 °C before use.
Chromatin immunoprecipitation was performed as previously described68 with the following modifications. Briefly, for each species 100 ng of chromatin was used for immunoprecipitation. The pool of chromatin was incubated for 14–16 h at 4 °C with 5 µl (1:50 dilution) of anti-H3K4me1 (Cell Signaling, catalogue no. 5326), 6 µl (3.4 µg) of anti-H3K4me2 (Abcam, catalogue no. ab32356), 2.5 µl (1:100 dilution) of anti-H3K4me3 (Millipore, catalogue no. 07-473), 5 µl (5 µg) of anti-SMC1 (Thermo Fisher, A300-055A) or 2 µl (2 µg) of anti-H3 (Abcam, catalogue no. ab1791) and recovered using a 1:1 mix of Protein A (Sigma-Aldrich, catalogue no. 16-661) and Protein G (Sigma-Aldrich, catalogue no. 16-662) magnetic beads. Immunoprecipitated complexes were washed, reverse crosslinked for 3 h at 68 °C, deproteinased and then purified using Ampure XP beads (Beckman Coulter, catalogue no. A63881). Final libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit (New England BioLabs) according to the manufacturer’s protocol. ChIP–seq libraries were sequenced on Illumina NextSeq 500 sequencer using High-Output 75 cycles.
MARS-seq library preparation
Single-cell libraries were prepared from freshly dissociated and sorted choanocytes of E. muelleri as previously described8. To collect cells for MARS-seq libraries, 7 days posthatching, sponges were fed with 0.5 µm fluorescent carboxylate-modified FluoSpheres (Invitrogen, catalogue no. F8813). Animal tissues were dissociated and prepared for sorting as described above for ATAC-seq. Dissociated cells were sorted through FACS into four 384-well MARS-seq plates. In total, 1,536 single-cell libraries were prepared and sequenced on an Illumina NextSeq 500 using High-Output 75 cycles.
Chromatin proteomics
Chromatin proteomics samples were prepared as previously descibed69 with minor modifications. Briefly, double-crosslinked cells of M. leidyi (1 million per replicate) were solubilized in 1 ml of lysis buffer (4 M guanidine thiocyanate, 100 mM Tris-HCl (pH 8.0), 10 mM EDTA, 2% N-lauroylsarcosine sodium salt) and incubated for 10 min. Next, before adding DNA-binding beads (Invitrogen, catalogue no. 37002D), cell lysate was mixed with 1 ml of 2-propanol. The beads were separated on a magnet and the supernatant was saved as the unbound control. The beads were washed using 1 ml of wash buffer (1:1 lysis buffer to 2-propanol ratio), transferred to a 1.5-ml sonication tube and washed again with 1 ml of 80% ethanol. The chromatin was then eluted in 200 µl of 10 mM Tris-HCl (pH 8.0) containing proteinase inhibitors (Roche, catalogue no. 04693132001) and sonicated using a Bioruptor Pico at 4 °C for 3 cycles (30 s ON, 30 s OFF). To remove RNA-binding proteins, RNase A (Roche, catalogue no. 10109142001) was added to the sonicated samples, which were then incubated at 37 °C with agitation in a thermomixer. Afterwards, chromatin was re-bound to the beads by adding 250 µl of lysis buffer, vortexing and then sequentially adding 300 µl of 2-propanol. The beads were washed twice with 1 ml of 80% ethanol, and proteins were digested on the beads using trypsin (Promega, catalogue no. V5111) and LysC (NEB, catalogue no. P8109S).
Chromatographic and mass spectrometric analysis
Samples were analysed using an Orbitrap Eclipse mass spectrometer (Thermo Fisher Scientific) coupled to an EASY-nLC 1200 (Thermo Fisher Scientific (Proxeon)). Peptides were loaded directly onto the analytical column and were separated by reversed-phase chromatography using a 50-cm column with an inner diameter of 75 μm, packed with 2-μm C18 particles (Thermo Fisher Scientific, catalogue no. ES903).
Chromatographic gradients started at 95% buffer A (0.1% formic acid in water) and 5% buffer B (0.1% formic acid in 80% acetonitrile) with a flow rate of 300 nl min−1 and gradually increased to 25% buffer B and 75% A in 52 min and then to 40% buffer B and 60% A in 8 min. After each analysis, the column was washed for 10 min with 100% buffer B.
The mass spectrometer was operated in positive ionization mode with nanospray voltage set at 2.4 kV and source temperature at 305 °C. The acquisition was performed in data-dependent acquisition mode and full mass spectrometry scans with one micro-scan at resolution of 120,000 were used over a mass range of m/z 350–1,400 with detection in the Orbitrap mass analyser. Automatic gain control was set to ‘standard’ and injection time to ‘auto’. In each cycle of data-dependent acquisition analysis, following each survey scan, the most intense ions above a threshold ion count of 10,000 were selected for fragmentation. The number of selected precursor ions for fragmentation was determined by the ‘Top Speed’ acquisition algorithm and a dynamic exclusion of 60 s. Fragment ion spectra were produced by means of high-energy collision dissociation at normalized collision energy of 28% and they were acquired in the ion trap mass analyser. Automatic gain control and injection time were set to ‘Standard’ and ‘Dynamic’, respectively, and an isolation window of 1.4 m/z was used. Digested bovine serum albumin (NEB, catalogue no. P8108S) was analysed between each sample to avoid sample carryover and to assure stability of the instrument, and Qcloud70 was used to control instrument longitudinal performance during the project.
Acquired spectra were analysed using the Proteome Discoverer software suite (v.2.5, Thermo Fisher Scientific) and the Mascot search engine (v.2.6, Matrix Science71). The data were searches against M. leidyi database and a list of common contaminants (16,042 entries)72 as well as all the corresponding decoy entries. For the peptide identification a precursor ion mass tolerance of 7 ppm was used for the MS1 level, trypsin was chosen as enzyme and up to three missed cleavages were allowed. The fragment ion mass tolerance was set to 0.5 Da for MS2 spectra. Oxidation of methionine and N-terminal protein acetylation were used as variable modifications whereas carbamidomethylation on cysteines was set as a fixed modification. False discovery rate in peptide identification was set to a maximum of 1%.
Peptide quantification data were retrieved from the ‘Precursor ions quantifier’ node from Proteome Discoverer (v.2.5) using 2-ppm mass tolerance for the peptide extracted ion current. The obtained values were used to calculate protein fold-changes and their corresponding P value and adjusted P values.
DAP-seq (DNA affinity purification sequencing) library preparation
The DNA-binding domains of candidate zf-C2H2 proteins were cloned from complementary DNA (cDNA) library of M. leidyi into the pIX-HALO vector using NEBuilder HiFi DNA Assembly Master Mix (NEB, catalogue no. E2621). The obtained HALO-fusion constructs were translated using the TnT SP6High-Yield Wheat Germ Protein Expression System (Promega, catalogue no. L3260). Next, an adapter-ligated DNA library was prepared from native gDNA of M. leidyi using NEBNext Ultra II FS DNA library prep kit (NEB, catalogue no. E7805) or PCR amplified gDNA. The binding to HALO-zf-C2H2 fusion proteins and recovery of adapter-ligated gDNA libraries was performed as described in ref. 73. The generated DAP-seq libraries were sequenced in paired-end mode on an Illumina NextSeq 500 using High-Output 75 cycles.
De novo genome assembly and scaffolding
We made preliminary genome assemblies of C. owczarzaki from Oxford Nanopore reads basecalled by Guppy v.6.0.1 using NextDenovo v.2.5.0 (ref. 74), Flye v.2.9.0 (ref. 75) and NECAT v.0.0.1 (ref. 76), which produced 20, 141 and 56 contigs including the mitochondrial genome, respectively. For the Flye assembly, we only used 5,000 bp or longer reads. We then integrated the three assemblies by manually comparing them to each other, with a help of reciprocal large-scale alignments generated with minimap2 (ref. 77). The integrated assembly was polished with the Nanopore reads using Flye ten times, and with Illumina reads78 using HyPo79 twice. A chromosome-scale duplication, which was in the end included in chromosome 15 after the 3D assembly (Extended Data Fig. 2d), was temporarily removed before annotating the genome. Finally, we manually inspected the whole assembly sequence together with the mapped Illumina data, Nanopore data and the previous Sanger sequence data25, and navigated them on the Integrative Genomic Viewer (IGV)80 to find and fix errors occurred during the consensus calling. We also manually phased chimeric haplotypes for some genes using the long reads. In total, 7,937 nucleotides were manually inserted or deleted at 430 sites and 1,193 nucleotides were substituted at 1,081 sites.
We produced two new genome assemblies for E. muelleri and M. leidyi. In both cases, we used Oxford Nanopore reads after base call correction using Guppy v.5.0.17 (using the dna_r9.4.1_450bps_sup_prom.cfg configuration the super-accurate base calling model, and a filtering reads with min_qscore=10). Then, we used two different long-read assemblers (Flye v.2.9-b1768, ref. 75 and Shasta v.0.8.0, ref. 81) and various assembly strategies (filtering by read length at 0, 10 and 50 kb), and selected the best resulting draft assemblies for each species. To that end, we evaluated the contiguity (measured using the contig N50), completeness and occurrence of uncollapsed haplotypes for each draft (Extended Data Fig. 2a,b). Contiguity was evaluated using total assembly length and contig N50. Completeness was measured with the fraction of conserved orthologues recovered by BUSCO v.5.1.2 (ref. 82) (using the genome mode and the metazoa_odb10) and the fraction of mappable genes from the original assemblies (mapped using Liftoff v.1.6.1, ref. 83). The presence of uncollapsed haplotypes was assessed with the distribution of per-base sequencing depths, calculated using the pbcstat utility in purge_dups v.1.2.5 (ref. 84) (for which we remapped the input reads to the assembly with minimap2 2.18-r1015 (ref. 77), using the -x map-ont preset for long-read mapping) (Extended Data Fig. 2e,f).
The best drafts for each species were produced using the following parameter combinations: (1) for E. muelleri, we used the Shasta assembler with the Nanopore configuration (–config Nanopore-Oct2021 flag), without filtering by read length (estimated sequencing depth roughly 100×) and (2) for M. leidyi, we used Flye with reads filtered at 50 kb (estimated sequencing depth roughly 150×), the raw Nanopore read configuration (–nano-raw flag) and an estimated total assembly size of 200 Mb.
Then, we used purge_dups to collapse putative uncollapsed haplotypes in each assembly, in the following manner: (1) we split the assembly into contigs with the split_fa utility; (2) we aligned the genome to itself with minimap2 and the -x asm5 preset; (3) we used the read alignments to the unsplit assembly (produced with minimap2 -x map-ont) to obtain the sequencing depth histogram and calculate coverage cutoffs with pbcstat and calcuts, respectively; (4) we used these cutoffs and the mapped reads to remove haplotigs and overlaps for the draft, with purge_dups proper and using two rounds of alignment chaining (-2 flag) and finally (5) we reevaluated the assembly quality using per-base sequencing depth distributions (above) and reductions in the fraction of duplicated BUSCO orthologues.
Chromosome-level assembly
To obtain chromosome-level genome assemblies, generated Micro-C libraries were mapped to de novo draft genome assemblies (C. owczarzaki, E. muelleri and M. leidyi) or current genome assemblies (T. adhaerens ASM15027v1, ref. 22, S. arctica24, S. rosetta GCA_000188695.1, ref. 26, C. collaboinventa85) using Juicer v.1.6 (ref. 86) with an option -p assembly. Proximity ligation alignments were used by 3D de novo assembly pipeline87 to order and orient available contigs into chromosomes with the following parameters: S. arctica -r 3 –editor-repeat-coverage 10, C. owczarzaki -r 0 –editor-repeat-coverage 4, S. rosetta -r 3 –editor-repeat-coverage 2, E. muelleri -r 2 –editor-repeat-coverage 10, M. leidyi -r 2 -i 1000 –editor-repeat-coverage 2, T. adhaerens -r 3 –editor-repeat-coverage 2 and C. collaboinventa -r 3 –editor-repeat-coverage 2. The resulting assemblies were manually reviewed and corrected with Juicebox Assembly Tools88 (Extended Data Fig. 2c–f). Finally, chromosome-level genome assemblies were polished with Medaka (v.1.5.0) to correct possible sequence errors such as indels and mismatches, as follows: (1) first, we mapped the Nanopore reads to the chromosome-level assembly using the minimap2-based mini_align utility; (2) we then used Medaka consensus to obtain consensus sequences, specifying a batch size of 200 (–batch 200 flag) and the r941_prom_sup_g507 configuration (–model flag) and (3) we merged the consensus and variant calls for all chromosomes into a polished assembly using Medaka stitch.
Genome annotation
To annotate the C. owczarzaki genome, we did not mask the repeats because the intergenic regions are very small25 and, thus, masking only increased annotation failure on duplicated genes. We used BRAKER2 (ref. 89) with OrthoDB90 protein sequence collections as hint data, as well as with RNA-seq data from a previous study61. The three preliminary annotations, evidenced by metazoan proteins, protozoan proteins and RNA-seq data, were combined with TSEBRA91, giving rise to 9,069 annotated transcripts. Finally, we manually searched and fixed wrong annotations by navigating the assembly on IGV80, comparing the combined annotation with the three preliminary annotations together with the mapped RNA-seq data. By this careful inspection, we modified or newly annotated 1,871 transcripts including alternatively spliced ones. Compared to the previously published proteome25 (v.2), only 4,076 out of 8,792 proteins (including alternatively spliced ones) had completely matched sequences to the those predicted in this study, allowing simple amino acid mismatches probably accounting for polymorphisms.
To annotate M. leidyi genome we first downloaded developmental Illumina RNA-seq samples (GSE93977), trimmed them with fastp and built a de novo Trinity assembly, which was mapped to the genome using gmap92. The RNA-seq was also directly mapped to genome using HISAT2 (ref. 93) with the –dta parameter, and genome-based transcriptomes were built for each sample using StringTie94. Merged mapped RNA-seq samples were then used to find high-quality intron junctions using Portcullis. The combination of Trinity, StringTie and Portcullis intron junctions were then fed to Mikado for transcript selection. The best resulting gene models based on mapping to UniProt were then used to train an Augustus model for M. leidyi. Augustus was used for an ab initio gene prediction, using exonic hints from Mikado, intron hints from Portcullis and coding sequence hints from a MetaEuk95 run with query fasta files combining proteins from H. californensis and UniProt. Mikado transcripts and Augustus gene models were then merged using EVidenceModeler (scores of 10 for Mikado transcripts and 2 for Augustus gene models). The resulting gene models were updated with PASA96 to incorporate the untranslated regions from the Mikado transcripts.
To annotate S. arctica, S. rosetta, E. muelleri, T. adhaerens and C. collaboinventa genome assemblies, gene models from previous assemblies were mapped onto new coordinates using Liftoff (v.1.6.1)83 with -overlap 1 -flank 1 options.
Repeat annotation
Repetitive sequences and transposable elements were annotated using EDTA (v.2.1.0)97 with the following parameters: –sensitive 1 –anno 1 (Extended Data Fig. 2g). For H. sapiens, we used RepeatMasker (v.open-4-0-3) annotation of GRCh38 genome released by UCSC.
DNA methylation calling from Nanopore long-read sequencing data
The fast5 files obtained from the PromethION were used as input for Megalodon (v.2.5), with the Remora model dna_r9.4.1_e8 sup for 5hmc_5mc modification only on CG dinucleotides. We then built bigwig files using the bedGraphToBigWig tool from UCSC. The Megalodon CG methylation calls were compared to previously published Whole-Genome Bisulfite Sequencing remapped to the new reference genomes using Bismark (SRR8346013 and SRR10356110)21,48. Both data sources were congruent, yet Nanopore had deeper and broader coverage, we used Megalodon methylation data for subsequent analysis.
Micro-C data processing
Micro-C data were processed using the 4D Nucleome processing pipeline98. Briefly, raw reads were mapped to the reference genome using bwa mem (v.0.7.17-r1188) with the -SP5M option. The mapped reads were sorted and filtered with pairtools (v.0.3.0)99. Pairs that mapped within a 2-bp distance from each other were considered duplicates. We also discarded reads mapping within the distance of 200 bp, which eliminates self-ligated pairs and reads mapping to adjacent nucleosomes. Only uniquely mapping pairs and 5′ most unique alignments of multiple ligations pairs were aggregated into 200-bp bin contact matrices and multiresolution .cool or .hic files (Extended Data Fig. 1d). Contact matrices were normalized with cooler (v.0.8.11)100 using the iterative correction and eigenvector (ICE) balancing method101 for .cool files or with Juicer tools86 using Knight–Ruiz balancing102 for .hic files. All contact heatmaps were visualized with either Cooltools (v.0.5.1)103 or Coolbox (v.0.3.8)104 and genome assembly heatmaps were visualized using HiGlass105.
Reproducibility between replicates was estimated using the stratum-adjusted correlation coefficient (SCC) implemented in HiCRep106 at resolutions of 1, 2, 5, 10, 25 and 50 kb (Extended Data Fig. 1f). The SCC scores were averaged across chromosomes. Biological replicates with SCC score estimated above 0.7 at resolutions equivalent to roughly 20,000 bins per species genome (resolution of 10 kb for S. arctica, 1 kb for C. owczarzaki, 2 kb for S. rosetta, 10 kb for E. muelleri, 10 kb for M. leidyi, 5 kb for H. californensis, 5 kb for T. adhaerens, 5 kb for C. collaboinventa and 10 kb for N. vectensis) were pooled to obtain final chromatin interaction matrices. Technical replicates were first merged, deduplicated and only then combined into the final contact maps.
The decay of the average contact frequency over genomic distance from 1 kb to 100 Mb was calculated using Cooltools (v.0.5.1)103. The decay curves were calculated for each chromosome separately, and then averaged across chromosomes (Extended Data Fig. 1g).
Compartment analysis
Compartment analysis was performed on observed-over-expected contact maps at resolutions equivalent to 5,000, 10,000, 20,000 and 50,000 bins per species genome (Extended Data Fig. 3a) using Cooltools eigs-cis103. We visually examined calculated eigenvectors, and, for each organism, the E1 vector corresponded to the compartmentalization pattern of contact maps. Active (A) and inactive (B) compartment types were assigned by GC content (for all species except C. owczarzaki) or H3K4me3 chromatin signal for C. owczarzaki, such that higher GC regions or positively correlated with H3K4me3 signal regions correspond to A compartment. Saddle plots were generated using the Cooltools saddle module. Specifically, the eigenvectors were sorted from lowest to highest value and combined into 40 groups according to their eigenvector value. The first (bottom 2.5% E1 values) and last (top 2.5% E1 values) groups were ignored to exclude potential outliers. The observed-over-expected value of the remaining 38 groups was averaged across all bins and chromosomes and visualized as saddle plots.
Compartment strength was calculated as the ratio of homotypic (AA + BB) over heterotypic (AB + BA) compartment contacts. We choose the top 20% of observed-over-expected values for both homotypic and heterotypic interactions. To assess the error in estimating compartment strength, we compared the compartments strength across different resolutions as well as performed visual inspection of the contact maps (Extended Data Fig. 3b). The latter showed varying degrees of accuracy in identifying compartment types between species, with the algorithm performing particularly poorly on M. leidyi due to the lack of well-defined chromatin compartmentalization in this species in our sample. Therefore, we assigned an extra intermediate compartment I to the intermediate eigenvalues close to zero. To that end, we modelled the genome-wide eigenvalues distribution as a Gaussian mixture with three components using the normalmixEM function from the mixtools R package (v.2.0.0) as described107. The B–I and I–A thresholds were defined as intersection points between components (Extended Data Fig. 3c).
To characterize the distribution of genomic features in the A, I and B compartments, we calculated cumulative H3K4me3 chromatin signals and RNA-seq expression values for each compartment region. Furthermore, we estimated the percentage of bases annotated as transposable elements or coding gene regions within these compartments. All the values are presented as −log2(1 − the value’s quantile). Thus, a normalized value of six means that the coverage is in the upper 1–2−6 quantile, that is, in the upper 1/64th of the distribution (Extended Data Fig. 3d).
Insulation profiles and boundaries classification
To compute the insulation profiles, we first determined the optimal resolution and window sizes for a target genome. To that end, we calculated insulation scores using Cooltools insulation module103 at resolutions roughly equivalent to 50,000, 100,000, 200,000 and 400,000 genomic bins per species genome with a sliding window for each resolution that is ×5, ×10 and ×25 the applied resolution (Extended Data Fig. 4a). The resolution and two window sizes with maximum average insulation scores were considered optimal because they reflected the strongest partitioning of genomes into isolated domains. Insulation boundaries located within two bins of unmappable genomic region were removed.
Identified insulation boundaries were categorized into strong and weak using the peak prominence of their boundary strength distributions (Li threshold) as implemented in the Cooltools insulation score module. Strong boundaries were further annotated with overlapping genomic features that fall within one bin of the annotated feature from the insulation boundary. For example, if compartment boundaries were called at the resolution 5 kb, then the maximum distance to the closest insulation boundary is ±10 kb.
To estimate internal interactions within contact domains, rescaled pile-ups were generated using coolpup.py108. Contact domains were defined as valleys between two strongly insulated regions, which were not further from each other than 100 kb.
Loop calling and annotation
Chromatin loops were identified using SIP v.1.6.1 (ref. 109) on KR-normalized contact matrices. In M. leidyi, the SIP peak caller was applied with the following parameters: -norm KR -g 3.0 -min 2.0 -max 2.0 -mat 5000 -d 10 -res 400 -sat 0.01 -t 2000 -nbZero 6 -factor 4 -fdr 0.05 -isDroso false. For T. adhaerens, chromatin loops were called with the following parameters: -norm KR -g 5.0 -min 4.0 -max 4.0 -mat 5000 -d 20 -res 100 -sat 0.01 -t 2000 -nbZero 6 -factor 4 -fdr 0.05 -isDroso false; for C. collaboinventa: -norm VC_SQRT -g 1.5 -min 3.0 -max 3.0 -mat 5000 -d 20 -res 500 -sat 0.01 -t 2000 -nbZero 6 -factor 2 -fdr 0.05 -isDroso false; for H. californensis: -norm KR -g 2.5 -min 3.0 -max 3.0 -mat 5000 -d 10 -res 1000 -sat 0.01 -t 2000 -nbZero 6 -factor 4 -fdr 0.05 -isDroso false. Identified loops were then filtered based on APSscore, removing high-intensity signals outside the normal distribution of APSscore values. This threshold ensured accurate removal of false positive regions that corresponded to structural genomic rearrangements, such as inversions or assembly artefacts. For H. californensis, we kept only annotated loops with values greater than ten. Chromatin loops in N. vectensis and focal chromatin contacts in E. muelleri were annotated manually.
Each loop anchor was assigned a promoter or enhancer identity based on their epigenetic signature. We calculated quantile normalized counts per million (CPM) coverage of H3K4me3, H3K4me2 and H3K4me1 ChIP signals in 1-kb (T. adhaerens, C. collaboinventa, N. vectensis), 2-kb (M. leidyi, E. muelleri, D. melanogaster) or 10-kb (H. sapiens) windows from a centre of a loop anchor.
METALoci autocorrelation analysis
METALoci110 (v.0.3.0) analysis was applied to explore the spatial distribution and autocorrelation of epigenetic signal in T. adhaerens contact maps. For each region of interest at 800 bp resolution, the top 20% pairs of contacts were used to create a two-dimensional graph layout by means of the Kamada–Kawai algorithm111 (Fig. 3b, top left panel) using the ‘metaloci layout’ with default parameters. Next, the signal of interest (H3K4me3 ChIP, ATAC, genic exon annotation) measured in the 800-bp genomic bin was mapped onto the built graph layout using ‘metaloci lm’ with default parameters. Spatial and epigenetic signals were embedded into Voronoi diagrams for enhanced visualization as a Gaudí plot (Extended Data Fig. 8c), and the local Moran’s index (LMI) analysis112,113 was applied for each bin of the Gaudí plot.
According to LMI analysis, each bin is assigned to one of the four distinct groups, called LMI quadrants, based on the signal value in a bin and average signal value in its neighbourhood. If a bin and its neighbourhood have similar amounts of signal (low or high), then this bin is assigned to a low–low (blue) or high–high (red) quadrant. If a bin and its neighbourhood have different amounts of signal, then the bin is assigned to a low–high (cyan) or high–low (orange) quadrants, respectively. Significantly colocalized bins according to LMI, in which a P value is obtained using a permutation test, are highlighted by colour in the LMI scatterplots (Extended Data Fig. 8c,j, left panels). An analogous colouring scheme is applied to the Voronoi diagrams of the Gaudí plots. Hence, the highlighted blue and red bins on a Gaudí plot represent bins in which the signal is significantly colocalized in the space. Thus, ATAC-seq, H3K4me3 and motif score signals are significantly enriched inside the nested focal contacts (Fig. 3b and Extended Data Fig. 8c), whereas exons are significantly enriched outside loop contacts (Fig. 3b). METALoci code is available at the GitHub repository (https://github.com/3DGenomes/METALoci).
Motif analysis
Loop anchor regions of M. leidyi (n = 8,523) and H. californensis (n = 478) were scanned for enriched motifs with HOMER114 in de novo motif discovery mode. As background sequences, we used random genomic regions of equivalent size and GC content (n = 38,810 in M. leidyi and n = 49,097 in H. californensis). For motif enrichment analysis in T. adhaerens (n = 3,557) and C. collaboinventa (n = 4,037), we scanned loop anchor regions using random genomic sequences of equivalent size (n = 32,004 in T. adhaerens and n = 36,178 in C. collaboinventa) as background. Loop anchor regions of N. vectensis (n = 327) were scanned for enriched motif using random genomic sequences (n = 45,268) of equivalent size and GC content as background. In addition, we used ATAC-seq accessible neuronal promoter regions (n = 22,961) as background to scan for enriched motifs in genomic regions that overlap ATAC-seq peaks located at the non-loop insulation boundaries (n = 9,016). To annotate genomic regions with identified motifs, we used the monaLisa package115, selecting percentile threshold of motif scores by comparing the motif score distributions in the target regions with genome-wide motif score distributions (Extended Data Figs. 7d,f, 8i and 9f).
Whole-genome alignment and sequence conservation analyses
We evaluated the degree of sequence conservation of the M. leidyi genome by comparing it to other ctenophores (B. microptera, P. bachei and Hormiphora californiensis). To that end, we first aligned all genomes to each other using Cactus v.2.6.4 (ref. 116), following a progressive approach guided by the species trees of ctenophores, namely: ((M. leidyi, B. microptera), (H. californiensis, P. bachei)). Second, we used the hal2maf utility from the HAL toolkit v.2.2 (ref. 117) to create MAF (multiple alignment format) alignments of each chromosome, using M. leidyi as reference. To identify conserved regions in these genomes, we used the rphast v.1.6.1 implementation of the Phast toolkit118, as follows: (1) we used phyloFit119 to create an initial null model of neutral change based on the fourfold degenerate codon positions of each genome’s coding regions, using a general reversible nucleotide transition matrix and the predefined species tree; (2) we used phastCons to optimize this model using the expectation–maximization procedure, re-estimating transition probabilities and tree parameters at each step (the optimization step was performed using only the longest chromosome in each genome).
Loop synteny analysis in M. leidyi
We evaluated the degree of syntenic conservation of the loop regions in M. leidyi compared to the other ctenophore genomes, and compared it to that of length-matched regions not involved in loops. To that end, we first identified orthologous genes across the four ctenophore species (M. leidyi, B. microptera, P. bachei and H. californensis) using Broccoli v.1.2 (ref. 120) to obtain orthologous gene pairs (step 4), using predicted peptide sequences as input (longest isoform per gene only). Within Broccoli, we used the maximum-likelihood gene tree inference algorithm (based on IQ-TREE121) and set a k-mer length of 10,000 to avoid the removal of paralogous sequences from the analysis. Second, we mapped pairs of loop anchor regions from M. leidyi (2,353 pairs of promoter–enhancer and 99 promoter–promoter loops, n = 2,452 in total) to their closest overlapping genes, and used these genes and their orthologs as anchors to map these regions to the other ctenophore genomes. In parallel, we randomly selected length-matched, non-loop overlapping regions from the M. leidyi genome to compare their synteny conservation with that of loop regions (using the randomizeRegions function in the regioneR R package122 to select 3× background regions, n = 7,356). Then, for each pair of species, we evaluated the synteny conservation of the foreground (loop) and background regions (random non-loops) from the point of view of the flanking synteny-anchoring genes, using two different metrics: (1) the fraction of shared orthologous genes between the flanking genes across all regions in the foreground and background sets (testing the significance of the difference using a χ2 test for given probabilities) and (2) the distributions of per-region shared orthologs (tested using the one-sided Wilcoxon rank sum test with continuity correction).
ATAC-seq analysis
We used previously published datasets of ATAC-seq for C. owczarzaki17, S. rosetta39, T. adhaerens123, C. collaboinventa123, D. melanogaster124 and H. sapiens125 as well as newly generated datasets for N. vectensis, M. leidyi and E. muelleri. Sequenced reads were demultiplexed and converted to fastq files using bcl2fastq v.2.20 Illumina. Raw reads were filtered and trimmed with Trimmomatic (v.0.39)126 before mapping to the reference genome with bwa mem (v.0.7.17-r1188) and duplicates were marked with bamsormadup from biobambam2 (https://github.com/gt1/biobambam2). Using deeptools alignmentSieve aligned reads were filtered and shifted with -ATACshift, which corresponds to mate reads being shifted +4 and −5 bp for positive and negative strands, respectively. To generate nucleosome-position data tracks, nucleosome-free and nucleosome-bound regions were defined using the following length thresholds 0–120 and 150–240 bp, respectively. ATAC peaks were called with MACS2 (ref. 127) on shifted nucleosome-free regions. Footrpint ATAC score was calculated using TOBIAS v.0.13.3 (ref. 128).
ChIP–seq analysis
We analysed publicly available dataset for D. melanogaster129 and H. sapiens130 and 34 newly generated ChIP–seq datasets as described below. Raw reads after removal of 3′-adapters and quality filtering with Trimmomatic (v.0.39)126 were aligned to the reference genome with bwa mem (v.0.7.17-r1188). Duplicated reads were marked with bamsormadup (https://github.com/gt1/biobambam2), and peaks were called using MACS2 (v.2.2.6)127. Aggregated density plots were visualized with deeptools (v.3.1.3)131.
DAP-seq analysis
Raw reads from amplified and native gDNA fragments bound by HALO-zf-C2H2 protein fusions were analysed as described for ChIP–seq. Motif enrichment analysis was performed using HOMER114 in de novo motif discovery mode for MACS2 identified narrow peaks resized to 300 bp (for CTEP1 n = 14,638; for CTEP2 n = 10,615). GC- and size-normalized random genomic regions were used as background (for CTEP1 n = 25,964; for CTEP2 n = 30,744).
RNA-seq analysis
We used previously published datasets of bulk poly-A enriched RNA-seq for S. arctica132, C. owczarzaki61, S. rosetta39, D. melanogaster124 and H. sapiens133 (Supplementary Table 2). To process data, raw reads were aligned to the reference genome using STAR (v.020201)134 in –quantMode to estimate the number of read counts per gene. In downstream analysis, gene counts were reported as −log2(1 − gene counts quantile).
MARS-seq analysis and single-cell expression atlases
Single-cell MARS-seq libraries generated previously8,135 were aligned to new reference genomes of E. muelleri (GCA_049114765.1), M. leidyi (GCA_048537945.1) and N. vectensis (GCA_932526225.1) using Liftoff or de novo annotated gene models. To improve single-cell RNA-seq quantification, gene annotations for E. muelleri and M. leidyi have been extended using GeneExt136. Briefly, MARS-seq alignment files have been subsampled to 100 M reads and MACS2 (ref. 127) was used to call peaks using default parameters. Intergenic peaks were filtered based on the 20th percentile of the genic peak coverage and each gene was extended to the most distant peak within 5,000 nucleotides. Metacell and clustering analyses were performed as previously described8. The single-cell expression atlas for T. adhaerens was obtained from a previously published dataset123.
Public datasets used in this study
All public datasets used in this study are listed in Supplementary Table 2. ATAC-seq, ChIP–seq and RNA-seq datasets were analysed as described above.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
