Lineage dynamics of the endosymbiotic cell type in the soft coral Xenia

Many corals harbour symbiotic dinoflagellate algae. The algae live inside coral cells in a specialized membrane compartment known as the symbiosome, which shares the photosynthetically fixed carbon with coral host cells while host cells provide inorganic carbon to the algae for photosynthesis1. This endosymbiosis—which is critical for the maintenance of coral reef ecosystems—is increasingly threatened by environmental stressors that lead to coral bleaching (that is, the disruption of endosymbiosis), which in turn leads to coral death and the degradation of marine ecosystems2. The molecular pathways that orchestrate the recognition, uptake and maintenance of algae in coral cells remain poorly understood. Here we report the chromosome-level genome assembly of a Xenia species of fast-growing soft coral3, and use this species as a model to investigate coral–alga endosymbiosis. Single-cell RNA sequencing identified 16 cell clusters, including gastrodermal cells and cnidocytes, in Xenia sp. We identified the endosymbiotic cell type, which expresses a distinct set of genes that are implicated in the recognition, phagocytosis and/or endocytosis, and maintenance of algae, as well as in the immune modulation of host coral cells. By coupling Xenia sp. regeneration and single-cell RNA sequencing, we observed a dynamic lineage progression of the endosymbiotic cells. The conserved genes associated with endosymbiosis that are reported here may help to reveal common principles by which different corals take up or lose their endosymbionts.

Article looked for two previously characterized cnidarian cells: the cnidocytes, which are used for prey capture and/or defence, and gastrodermal cells. The cells of cluster 11 express minicollagen and nematogalectin genes, which are markers of cnidocytes [18][19][20] (Fig. 2c). Further analysis revealed that cluster 11 contained two subclusters (Fig. 2d, Extended Data Fig. 3b). Minicollagen genes are expressed in both subclusters, whereas nematogalectin genes are preferentially expressed in one (Fig. 2e, Extended Data Fig. 3c). RNA in situ hybridization (ISH) confirmed the expression of a nemetogalectin gene to be more spatially restricted than that of Minicollagen 1 in Xenia pinnules (Fig. 2f, g, Extended Data Fig. 3d, e). Clusters 2, 12 and 16 express genes that encode collagens and proteases (Fig. 2h) that are known to be enriched in gastrodermis of Nematostella 18 . RNA ISH for Collagen 6, Astacin-like metalloendopeptidase 2 (both expressed by clusters 2 and 12) and the uncharacterized Xe_003623 gene (expressed by clusters 2, 12 and 16) confirmed the high expression of these genes in the gastrodermis ( Fig. 2h-j, Extended Data Fig. 3f-i). Thus, the clustering analyses and ISH identified cnidocytes and cells in the gastrodermis in Xenia.

Endosymbiotic cell type in Xenia sp.
To identify the cells that perform endosymbiosis, we took advantage of the autofluorescence of the member of the Symbiodiniaceae (Durusdinium) in our Xenia sp. (Methods). Using fluorescence-activated cell sorting (FACS), we separated alga-containing and alga-free Xenia cells (Fig. 3a, b) and performed bulk RNA-seq (Supplementary Table 3). By comparing these bulk transcriptomes with genes expressed in each cluster, we found that cells of cluster 16 exhibited the highest overall similarity to the alga-containing cells and most of the marker genes for cluster 16 (Supplementary Table 4) have a higher level of expression in alga-containing Xenia cells than that in alga-free Xenia cells (Fig. 3c, d, Supplementary Table 5). RNAscope ISH for two of the cluster-16 marker genes-one of which encodes a protein with lectin and kazal protease inhibitor domains (abbreviated LePin, encoded by a gene that we name LePin), and the other of which encodes Granulin 1-showed that these genes were expressed in alga-containing gastrodermal cells (Fig. 3e, f, Extended Data Fig. 4a, b). Additionally, on average 95% and 98% of alga-containing Xenia cells were positive for expression of LePin and Granulin 1, respectively (Extended Data Fig. 4c). On the basis of microscopy of cryopreserved tissue sections or FACS analyses, we estimated that on average 2-6% of Xenia cells contained algae and that tentacles have a higher percentage of alga-containing cells than do stalks (Extended Data Fig. 4d, Methods). This is consistent with the cluster-16 endosymbiotic cells being identified by scRNA-seq as a small fraction (382 cells, 1.4% of the total). Of the three gastrodermal cell clusters, cluster-16 cells therefore have a high likelihood of being a major cell type involved in endosymbiosis.
Among the top 89 marker genes enriched in the cluster-16 endosymbiotic cells, 67 encode proteins with domains of known or predicted functions, including receptors, extracellular matrix proteins, immune response proteins, phagocytosis and/or endocytosis proteins, or nutrient transporters (Extended Data Fig. 4e, Supplementary Table 4). Three proteins-encoded by CD36, DMTB1 and CUZD1-contain CD36 or scavenger receptor domains that are known to recognize a wide range of microbial surface ligands and mediate their phagocytosis, and that also modulate the innate immune response of the host 11,21,22 (Fig. 3g, Extended Data Fig. 5a). CUZD1 is the least understood, and is similar to DMBT1 in domain organization. DMBT1 functions in pattern recognition of microorganisms. In mammals, it is expressed on the surface of the gastrointestinal tract, where it recognizes polysulfated and polyphosphorylated ligands on microorganisms, represses the   inflammatory response and regulates the differentiation of gastrointestinal cells 23 . LePin and Granulin 1, which we used for ISH, have homologues in Exaiptasia, as well as stony and soft corals. Because LePin has an N-terminal signal peptide followed by multiple domains (including H-and C-type lectins and a Kazal-type serine protease inhibitor) (Extended Data Fig. 5b), it may confer selectivity for the Symbiodiniaceae. On the basis of previous studies of granulins in mammals 24 , Granulin 1 may modulate the immune response in Xenia endosymbiotic cells.
Phagocytosis of the Symbiodiniaceae by gastrodermal cells (which are of a similar size to these algal cells) requires substantial expansion of the host cell, but the genes that regulate this size expansion are unknown. Among the endosymbitoic marker genes that we found, Plekhg5 encodes a highly conserved RhoGEF (Fig. 3g, Extended Data Fig. 5c). In Xenopus, Plekhg5 localizes to the apical membrane of epithelial cells and recruits actomyosin to induce cell elongation and apical constriction 25 . Thus, Plekhg5 is a prime candidate for regulating the extension of the apical membrane to engulf algae of the Symbiodiniaceae during the early stages of phagocytosis in Xenia. Upon phagocytosis, algae of the Symbiodiniaceae are enclosed by the host membrane to form symbiosomes 26  to lysosomes 27 , the genes that are involved in the formation of the symbiosome remain unclear. Xenia sp. has two genes that encode lysosome-associated membrane glycoproteins, which are more similar to the previously characterized LAMP1 than to LAMP2 28 . In Xenia, Lamp1-L encodes a larger protein and is an endosymbiotic marker gene, whereas Lamp1-S encodes a smaller protein and is expressed across all cell clusters (Extended Data Fig. 5d, e). Because lysosome-associated membrane glycoproteins are known to regulate phagocytosis, endocytosis, lipid transport and autophagy 28 , Lamp1-L may regulate symbiosome formation and/or function (Fig. 3g). Several endosymbiotic marker genes encode enzymes that may promote the establishment of endosymbiosis or facilitate nutrient exchanges between alga and the host cell. For example, there are 17 genes that potentially participate in nutrient exchanges, as they encode transporters for sugar, amino acids, ammonium, water, cholesterol and choline (Fig. 3g

Lineage dynamics of endosymbiotic cells
To better understand the temporal dynamics of cluster-16 cells, we developed a Xenia regeneration model. We surgically cut away all tentacles from Xenia polyps and found that the stalks regenerated all tentacles in several days when cultured in the seawater from our aquarium that houses stock corals (Fig. 4a). Individual tentacles also regenerated into full polyps, but required a longer time (data not shown). BrdU labelling showed that some proliferated (BrdU + ) gastrodermal cells began to take up algae that were present either in the gastrodermis or in the seawater at day 4 of regeneration (Extended Data Fig. 6a, b). We performed scRNA-seq of the regenerating stalks and pooled the data with the scRNA-seq of non-regenerating samples (Methods). We used Monocle 2 to perform pseudotemporal ordering of all of the endosymbiotic Xenia cells 29 (Fig. 4b); Monocle 2 uses reversed graph embedding to construct a principal curve that passes through the middle of the cells in the t-SNE space. Because this trajectory analysis does not provide a direction of cell-state progression, we used velocyto 30 to determine the directionality of lineage progression of all cells and focused on the endosymbiotic cells in the regenerating sample. Velocyto calculates RNA velocity by comparing the number of unspliced and spliced reads, which measures the expected change in gene expression in the near future-thereby providing the directionality of cell-state change. This enabled the identification of early and late stages of endosymbiotic cells (green and red, respectively, in Extended Data Fig. 6c). The cell trajectory showed that the early and late-stage cells are mapped to relatively early and late pseudotime, respectively (Extended Data Fig. 6d). Thus, the pseudotime represents actual lineage progression. Modelling of gene expression revealed substantial changes along pseudotime. Further hierarchical clustering showed distinct gene-expression patterns, which helped to define five putative endosymbiotic cell states (Fig. 4c, Supplementary Table 6).
To further explore the cell dynamics in these five states, we compared single-cell transcriptomes to transcriptomes from the bulk RNA-seq of alga-containing or alga-free cells isolated by FACS, and plotted the expression correlation along pseudotime. State-3 cells showed the strongest correlation with the alga-containing cells, followed by state 2 and then state 1; state-4 and state-5 cells showed the least correlation (Fig. 4d). This suggests that state 3 represents mature, alga-containing cells. State-1 and state-5 cells showed correlations with alga-free cells (Fig. 4d). Given that these five states are present in our identified endosymbiotic cell type with a linear pseudotime progression, we hypothesize that state-1 cells are pre-endosymbiotic progenitors that can transition through state 2 to become state-3 mature alga-containing cells, and that state-3 cells could further transit through state 4 into state-5 post-endosymbiotic cells (Fig. 4e). In support of this, we found that the regenerating samples have higher percentages of state-1 (pre-endosymbiotic) and state-2 (transition 1) cells, and that the non-regeneration sample has more state-3 mature, state-4 (transition 2) and state-5 (post-endosymbiotic) cells (Fig. 4f).
We further verified our hypothetical endosymbiotic cell states in the regeneration paradigm by pulse-chase experiments (Methods). After cutting, Xenia sp. stalks were pulsed with EdU at day 3 and day 4 of regeneration. EdU was washed out, corals were allowed to continue regenerating and samples were collected at days 7, 9, 11, 13, 15, 17 and 19 (Extended Data Fig. 7a). Using FACS (Extended Data Fig. 7b-h), we calculated the percentages of EdU + alga-containing cells out of all alga-containing Xenia cells, and the percentages of all alga-containing Xenia cells out of all Xenia cells. We found an increase of EdU + alga-containing Xenia cells up to regeneration day 13, which may account for the increase in uptake of algae during tentacle growth (as tentacles have more alga-containing cells than the stalk) (Extended Data Figs. 4d, 7i). Thereafter, the percentage of total alga-containing cells remained constant, but the percentage of EdU + alga-containing cells gradually decreased (Extended Data Fig. 7i, j). Thus, these results support our hypothesis that the endosymbiotic cells progress from a progenitor state through an alga-uptake state and a mature alga-containing state, followed by loss of their algae.
Analysis of differentially expressed genes suggests the roles of each state in endosymbiotic cell lineage development and function. For example, the state-1 pre-endosymbiotic cells express WNT7b and WNT11, which may regulate progenitor-cell proliferation and differentiation through the Wnt signalling pathway 31,32 . Among 24 genes   preferentially expressed in state 3, 13 are endosymbiotic markers that are expressed at higher levels in the FACS-isolated alga-containing Xenia cells than that in the alga-free cells. By contrast, none of the genes preferentially expressed in state 5 is an endosymbiotic marker. Instead, state-5 cells preferentially express several oxidative-stress-response genes (see Supplementary Table 6 for detailed descriptions). Because increased oxidative stress is observed upon cellular ageing and during coral bleaching [33][34][35] , state-5 cells may represent a natural ageing state of endosymbiotic cells that are no longer able to hold on to their algae. Additional molecular studies exploring the function of the differentially expressed genes in each state are needed to further validate our five-state hypothesis.

Summary and outlook
Here we demonstrate the power of genomic and bioinformatic tools in studying coral biology. The Xenia sp. genome encodes essential components of RNA interference, such as Dicer and Ago, and DNA repair pathway proteins, which should enable the development of gene-manipulation tools to determine the mechanism of endosymbiosis. Although we focused on studying the endosymbiotic cell lineage, the regenerative processes for the other cell clusters can be similarly investigated in future analyses. Our studies suggest that Xenia endosymbiotic cells exist in five progressive states that are dynamic between homeostatic conditions and the regeneration process ( Fig. 4f). It will be important to further understand the endosymbiotic lineage progression under different environmental stressors and to test whether efficient recovery from bleaching relies on state-1 pre-endosymbiotic cells. It is also feasible to test whether forced regeneration by fragmenting bleached corals can stimulate the expansion of state-1 pre-endosymbiotic cells and the restoration of endosymbiosis.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-020-2385-7.

Genomic DNA isolation from Xenia sp.
To enable Nanopore DNA sequencing, we modified a protocol 37 that allowed the isolation of long DNA fragments. For each DNA preparation, one or two Xenia sp. colonies containing about 30 polyps were collected from the aquatic tank and washed 3 times for 5 min each with Ca 2+ -and Mg 2+ -free artificial seawater (449 mM NaCl, 9 mM KCl, 33 mM Na 2 SO 4 , 2.15 mM NaHCO 3 , 10 mM Tris-HCl, 2.5 mM EGTA, pH 8.0). Tentacles were cut away, as they secrete a lot of mucus (which affected the quality of the isolated DNA). The remaining stalks and the bases of individual Xenia colonies were placed in 100 μl DNAzol (Invitrogen) in a 1.5-ml microcentrifuge tube. The tissues were cut into small pieces by a scissor to make fragment sizes of about 1/10th of the original size. These fragments were further minced by a small pestle made for 1.5-ml microcentrifuge tubes (Fisher Scientific, 12-141-364). Then, 900 μl DNAzol was added, followed by vortexing the sample and then transferred to a 15-ml conical tube. Four millilitres of DNAzol and 50 μl of 10 mg/ml RNase A were then added to the tube and mixed, followed by incubation at 37 °C for 10 min. Then, 25 μl of 20 mg/ml proteinase K was then added, mixed and the tube was incubated at 37 °C for another 10 min. The sample was centrifuged at 5,000g for 10 min. The supernatant was transferred to another 15-ml tube. After adding 2.5 ml ethanol, the tube was gently mixed by inverting several times. The tube was left to stand at room temperature for 3 min followed by centrifugation at 1,000g for 10 min to pellet the DNA. The supernatant was discarded and the DNA pellet was resuspended in 500 μl TE (10 mM Tris-HCl, 1 mM disodium EDTA, pH 8.0). After the DNA had dissolved, 500 μl of phenol:chloroform:isoamyl alcohol (25:24:1) was added, and the tube was placed on the Intelli-Mixer RM-2S for mixing using programme C1 at 35 rpm for 10 min. The mixture was then transferred to a 2-ml phase-lock gel (QuantaBio, Cat. 2302820) and centrifuged at 4,500 rpm for 10 min. The aqueous phase was transferred into a new 2-ml tube, 200 μl 5 M ammonium acetate and 1.5 ml ice-cold ethanol were added followed by centrifugation at 10,000g for 10 min to pellet DNA. The pellet was washed twice with 1 ml 80% ethanol. After removing as much ethanol as possible, the DNA pellet was left to dry at 42 °C for 1 min, and then resuspended in 50 μl TE buffer.

Illumina sequencing
Genomic DNA prepared as in 'Genomic DNA isolation from Xenia sp.' was fragmented into about 400 bp, and libraries were made with ThruPLEX DNA-Seq kit (TaKaRa) according to the manufacturer's manual. These libraries were sequenced using the NEXseq500 platform with NextSeq 500/550 High Output Reagent Cartridge v2 (Illumina).

Nanopore sequencing
Genomic DNA was used to build Nanopore sequencing libraries with Ligation Sequencing Kit (SQK-LSK108, Oxford Nanopore Technologies), following the manufacturer's manual. For the first three runs, genomic DNA was not fragmented, to generate long reads. To obtain more reads, for the fourth run of Nanopore sequencing, genomic DNA was sheared to 8-10 kilobases by g-TUBE (Covaris, 520079

Hi-C
To perform Hi-C on Xenia sp. tissue, we modified a previously published protocol for nuclear in situ ligation 13 , as described in detail.
Nuclear permeabilization and chromatin digestion (step 2). The pellet from step 1 was resuspended in 10 ml ice-cold Hi-C lysis buffer (10 mM Tris, pH 8.0, 10 mM NaCl, 0.2% NP-40, 1× protease inhibitors cocktail (Roche, 04693132001)) and rotated for 30 min at 4 °C followed by centrifugation at 1,000g for 5 min at 4 °C. The pellet was resuspended with 1 ml ice-cold 1.2× NEB3.1 (120 μl NEB3.1 to 880 μl ddH 2 O) buffer and transferred to a 1.5-ml microcentrifuge tube followed by centrifugation at 1,000g for 5 min at 4 °C. The pellet was washed again with 1 ml ice-cold 1.2× NEB3.1 followed by centrifugation. After removing the supernatant, 400 μl 1.2× NEB3.1 buffer and 12 μl of 10% SDS were added to the pellet. P200 pipette tip was used to thoroughly resuspend and dissociate the pellet. The mixture was then incubated at 65 °C for 10 min at 950 rpm in a Thermomixer (Eppendorf). After cooling the mix on ice for 5 min, 40 μl 20% Triton X-100 was added to the mixture to neutralize the SDS. After carefully mixing by pipetting with a P200 pipette tip and inverting the tube several times, the mixture was then incubated at 37 °C for 60 min with rotation (950 rpm) in a Thermomixer. To digest the crosslinked genomic DNA, 30 μl of 50 U/μl BglII (NEB R0144M) was added to the mixture and incubated overnight at 37 °C with rotation at 950 rpm in a Thermomixer.

Fill in 5′ overhang generated by BglII digestion with biotin (step 3).
A nucleotide mix containing dATP, dGTP and dTTP was made by adding 1 μl each of 100 mM dATP, dGTP and dTTP into 27 μl ddH 2 O. To the 480.0 μl BglII-digested nuclear preparation from the above step 2, 4.5 μl of the nucleotide mix, 15 μl 1 mM biotin-16-dCTP (Axxora, JBS-NU-809-BIO16) and 10 μl 5 U/μl Klenow (NEB, M0210L) were added followed by incubation at 37 °C for 90 min with intermittent gentle shaking at 700 rpm for 10 s after every 20 s using Thermomixer. The tube was also taken out and inverted every 15-20 min. After this incubation, the mixture was kept on ice.

Reverse crosslink and DNA isolation (step 5).
To the reaction mixture from step 4, 25 μl of 20 mg/ml proteinase K (Invitrogen, 25530-049) was added and the mixture was divided equally into 8×1.5-ml microcentrifuge tubes (about 950 μl per tube). The tubes were then incubated overnight at 65 °C with rotation at 950 rpm in a Thermomixer. The next day, 3 μl 20 mg/ml proteinase K was added to each tube followed by incubation at 65 °C for 2 h with mixing in Thermomixer. The mixtures were combined into one 50-ml conical tube. After cooling down to room temperature, 10 ml phenol (pH 8.0) (Sigma) was added and mixed by vortex for 2 min. The mixture was then centrifuged for 10 min at 3,000g (Sorvall Lynx 6000 centrifuge). The supernatant containing the DNA was mixed with 10 ml phenol:chloroform (1:1) (pre-warmed to room temperature) and vortexed for 2 min. The whole mixture was then transferred to a 50-ml MaXtract High Density tube (Qiagen, 129073) and centrifuged at 1,500g for 5 min (Sorvall Lynx 6000 centrifuge). The top phase containing the Hi-C DNA was transferred to a 50-ml conical tube and the volume (usually about 10 ml) was adjusted to 10 ml with 1× TE as needed. To pellet the DNA, 1 ml 3 M Na-acetate, 5 μl 15 mg/ml GlycoBlue (Invitrogen AM9515) and 10 ml isopropanol were added to the mixture and incubated at −80 °C for >1 h. The DNA was then pelleted by centrifugation at 17,000g for 45 min at 4 °C (Sorvall Lynx 6000 centrifuge). The Hi-C DNA pellet was resuspended in 450 μl 1× TE and transferred to a 1.5-ml microcentrifuge tube followed by adding 500 μl phenol:chloroform (1:1). After mixing by vortex, the mix was centrifuged at 18,000g for 5 min at room temperature. The top aqueous layer was collected into another tube followed by adding 40 μl 3M Na-acetate, 1 μl 15 mg/ml GlycoBlue (Invitrogen AM9515, 300 μl) and 1 ml ice-cold 100% ethanol. After incubating at −80 °C for >30 min, the DNA was centrifuged at 21,000g for 30 min at 4 °C. The DNA pellet was washed with freshly prepared 70% ethanol and air-dried, followed by dissolving in 45 μl EB (10mM Tris, pH 8.0). The contaminated RNA in the DNA preparation was digested by adding 0.5 μl 10 mg/ml RNaseA and incubated at 37 °C for 30 min.

Remove biotin from the free DNA (unligated DNA) ends (step 6).
To remove the biotin at the free DNA ends, 1.0 μl 10 mg/ml BSA (NEB, 100×), 10 μl 10× NEB 2.1 buffer, 1 μl 10 mM dATP, 1 μl 10 mM dGTP and 5 μl T4 DNA polymerase (NEB M0203S), and 42 μl water were added to 40 μl (about 3 μg) Hi-C DNA preparation from step 5. The mixture was divided into two equal aliquots in 2 PCR tubes and incubated at 20 °C for 4 h. Then, 2 μl of 0.5 M EDTA was added to each of the two tubes to stop the reaction. The Hi-C DNA was then purified using the Clean and Concentrator Kit (ZYMO, D4013) followed by elution with 50 μl EB.
Sequencing adaptor ligation (step 9). The 50-μl beads from step 8 was mixed with 3.75 μl sequencing adaptor (TruSeq RNA Sample Prep Kit v.2), 10 μl 1× T4 DNA ligase buffer, 3 μl T4 DNA Ligase (30 U/μl) (Thermo Scientific, EL0013) and incubated at room temperature for 2 h. The beads were collected by a magnetic stand followed by washing twice with 400 μl 1× binding buffer + 0.05% Tween, 200 μl 1× binding buffer, and then 100 μl EB. The beads were resuspended in 40 μl EB. To release the DNA from the beads, the mixture was incubated at 98 °C for 10 min and then centrifuged at 500 rpm to pellet the streptavidin beads.

Sequencing library preparation (step 10). TruSeq RNA Library Prep
Kit was used to make DNA sequencing library (eight PCR cycles were used) and the DNA was sequenced by NextSeq 500.

scRNA-seq
For each of the six scRNA-seq library preparation,  39 and Hydra 40 scRNA-seq datasets (Supplementary Table 1). This suggests that, even using v.3 chemistry, the presence of seawater and/or Xenia-sp.-specific features may contribute to the reduced scRNA-seq quality.
We noticed the mapping rate in v.3 chemistry is lower than in v.2 chemistry. We sequenced more reads for the v.3 libraries, because v.3 captured more total cells and more RNA molecules per cell. Although we sequenced more for the v.3 libraries, we obtained lower sequence saturation (on average, 79.6% in v.3 libraries and 92.6% in v.2 libraries). Because the v.2 and v.3 reagent contents are proprietary information, it is difficult for us to assess why the two methods gave different results. Regarding our library preparation, the v.3 method entailed 22% of the total volume coming from the cell suspension in the Ca 2+ -free seawater, while in the v.2 method, 17.4% of the total volume came from the Ca 2+ -free seawater cell suspension. We therefore know that one difference between the two methods is that the salt concertation in v.3 library preparation is higher than that in the v.2 library preparation. The higher salt concentration in v.3 could lead to a higher RNA extraction efficiency in the v.3 library preparation, which could contribute to the difference between our v,2-and v,3-based scRNA-seq. Althought the 10x platform worked well for the Xenia sp. we studied here, it is important to keep in mind that modifications may be needed for successful scRNA-seq for other marine cnidarians.

Quantification of endosymbiotic Xenia cells by microscopy and FACS
To quantify the endosymbiotic cell percentage in Xenia, we first applied a microscopy-based strategy. By imaging cryo-preserved tissue sections stained with 1 μg/ml DAPI that labelled all nuclei, we determined the total number of Xenia cells per section by counting the number of Xenia cell nuclei: these nuclei are easily differentiated from the alga nuclei when overlapped with the autofluorescence signal in far red channel from algae. The number of Xenia cells containing alga is estimated by counting the number of algae surrounded by Xenia tissue. The estimated percentage by this method is on average 2-6%, depending on whether the sections were taken from stalks or tentacles (Extended Data  Fig. 4d). The limitation of this method is that some algae that appear to be inside the tissue may be between Xenia cells and not inside cells. Therefore, this estimate could represent an upper limit of the percentage of alga-containing Xenia cells.
In the second method, we used FACS to separate free algae and algae contained inside the Xenia cells. Xenia polyps were dissociated into single-cell suspension with the same preparation method as described in 'scRNA-seq'. The cells were fixed with 1% (final concentration) formaldehyde on ice for 1 h, followed by 0.2% Triton X-100 permeablization and 1 μg/ml DAPI staining. We first separated free algae and alga-containing Xenia cells according to the algae autofluorescence in the Cy5.5 channel. Free algae and algae inside Xenia cells should have different forward scatter (FSC) and side scatter (SSC) signals because the alga inside Xenia cells is enclosed by the Xenia cellular membrane structure. Thus, we used FSC and SSC to further gate the total population of algae into two subpopulations. Microscopy analyses showed that this gating separated free algae and alga-containing Xenia cells.
To determine the total Xenia cell number, Xenia cells together with algal cells were gated according to DAPI-positive signal followed by gating with the Cy5.5 signal. The total Xenia cells were calculated as alga-free Xenia cells plus the alga-containing Xenia cells. On the basis of these FACS analyses, we were able to estimate the percentage of alga-containing Xenia cells in Xenia polyps to be about 2% of total Xenia cells. The illustration of this FACS sorting can be found in Extended Data Fig. 7b-g. Because the procedure of single-cell dissociation may cause an alga-containing Xenia cell to lose its alga, the approximately 2% of alga-containing Xenia cells obtained by the FACS method probably represents an underestimation. Thus, we estimate the fraction of alga-containing Xenia cells to be about 2-6%.
Bulk RNA-seq Total RNA was isolated from 3 polyps, 32 tentacles or 6 stalks by RNeasy Plus Mini Kit (Qiagen). To obtain additional transcriptomes from different cell types, we dissociated coral tissue into individual cells according to a previously published method 41 and subjected the dissociated cells to OptiPrep-based cell separation 42 . Cells with different densities were separated into four layers, and RNA was isolated from each layer with RNeasy Plus Mini Kit (Qiagen). For transcriptome of FACS-isolated alga-containing and alga-free cells, three polyps were dissociated with the same protocol as used in the scRNA-seq and the dissociated cells were subjected to FACS. Cy5.5-positive and -negative cells were collected as alga-containing and alga-free cells, respectively, and used for total RNA extraction as above. cDNA libraries were built according to TruSeq Stranded mRNA Library Prep Kit (Illumina) and subjected to Illumina NextSeq 500 for sequencing. For gene annotation, paired-end sequencing of 75 bp for each end was used. For FACS-isolated bulk-cell transcriptomes, single-end sequencing of 75 bp was used.
Xenia regeneration, BrdU labelling and EdU pulse-chase Individual Xenia sp. polyps were placed into a well of 24-well cell-culture plate (Corning) containing 1 ml artificial seawater from our aquatic tank. The polyps were allowed to settle in the well for 5-7 days before cutting away the tentacles. After cutting, there were a lot algae released into the seawater, which together with the free algae living inside the cavity of the coral could serve as alga reservoirs for the uptake of algae during regeneration.
For the BrdU labelling experiments, 0.5 mg/ml BrdU was added into the well 2 d before sample collection. The BrdU-labelled stalks were fixed by 4% PFA overnight, followed by washing with PBST (PBS+0.1% Tween 20) twice for 10 min each. The stalk was then balanced with 30% sucrose overnight followed by embedding in OCT, frozen in dry ice bathed in ethanol and subjected to cryo-sectioning. The slides were washed with PBS 3 times for 5 min each time followed by treating with 2 M HCl containing 0.5% Triton X-100 for 30 min at room temperature. The slides were then incubated with PBST (0.2% Triton X-100 in PBS) 5 min for 3 times each followed by blocking with 10% goat serum and then incubating with mouse anti-BrdU antibody (ZYMED, 18-0103, 1:200 dilution in 10% goat serum) at 4 °C overnight. Slides were washed with PBST 3 times for 10 min each followed by incubation with the secondary antibody (Invitrogen) for 1 h at room temperature and washing with PBST 3 times for 10 min each. The nuclei were counterstained with Hoechst 33342 and the signal was visualized using a confocal microscope (Leica). Clear BrdU signal in the nucleus labelled by Hoechst was counted as a BrdU + cell. If the Xenia BrdU + nucleus was juxtaposed to an alga, it was counted as an alga-containing BrdU + Xenia cell.
For EdU pulse-chasing experiments, the regenerating Xenia stalks were incubated with 1 mM EdU during regeneration day 3 and day 4. After washing out EdU, the coral was incubated with artificial seawater and samples were collected on regenerating days 7, 9, 11, 13, 15, 17 and 19. The samples were dissociated into single-cell suspensions followed by fixing with 1% formaldehyde at 4 °C overnight as described in 'scRNA-seq'. The fixed cells were pelleted at 800g for 5 min and further fixed with 4% PFA for two days to block the autofluorescence in the 488-nm channel. Then, the EdU click chemistry was carried out using the Click-iT EdU Cell Proliferation Kit (Invitrogen, C10337) according to manufacturer's protocol. The cells were further stained with DAPI, and then analysed by FACS as described in Extended Data Fig. 7 and 'Quantification of endosymbiotic Xenia cells by microscopy and FACS'.

Whole-mount RNA ISH
To perform RNA ISH on Xenia, we modified the whole-mount RNA ISH protocol for zebrafish 43 .
For making gene-specific sense or anti-sense probes, we designed primers (Supplementary Table 7) to genes of interest for PCR to amplify gene fragments from Xenia sp. cDNA. The T3 promoter sequence was added to the 5′ of the reverse primers so that the PCR products could be directly used for synthesizing anti-sense RNA probes by T3 RNA polymerase (Promega, P2083) using DIG RNA Labelling Mix (Roche, 11277073910). DIG-labelled RNA probes were purified by RNA Clean and Concentrator-5 (ZYMO), heated to 80 °C for 10 min, immediately transferred on ice for 1 min, and then diluted in Prehyb + buffer (50% formamide, 5× saline-sodium citrate buffer (SSC, 0.75M NaCl, 0.075M sodium citrate), 50 μg/ml heparin, 2.5% Tween 20, 50 μg/ ml single-stranded DNA (Sigma, D1626)) to a final concentration of 0.5 μg/ml, and stored at −20 °C until use.
Xenia polyps were relaxed in Ca 2+ -free seawater for 30 min and then fixed in 4% PFA in Ca 2+ -free seawater overnight at 4 °C. Fixed polyps were washed with PBST (0.1% Tween 20 in PBS) twice for 10 min each, and then incubated in 100% methanol at −20 °C overnight.The next day, the tissues were washed sequentially in 75%, 50% and 25% methanol for 5 min each and then washed in PBST for 10 min. They were then treated with 50 μg/ml proteinase K in PBST for 20 min followed by a brief wash in PBST. The tissues were post-fixed in 4% PFA at room temperature for 20 min and then washed with PBST 2 times for 10 min each. Prehybridization was performed in Prehyb + at 68 °C for 2 h, followed by incubation with probes in Prehyb + overnight at 68 °C. To probe gastrodermis markers, 2% SDS (final concentration) was added to help the probes to penetrate the tissue. After probes were removed, samples were washed sequentially in 2× SSC (0.3 M NaCl and 0.03 M sodium citrate) containing 50% formamide for 20 min twice, 2× SSC containing 25% formamide for 20 min, 2× SSC for 20 min twice, and 0.2× SSC for 30 min 3 times each, all at 68 °C. Then, samples were washed in PBST at room temperature for 10 min and incubated in DIG blocking buffer (1% ISH blocking reagent (Roche, 11096176001) in maleic acid buffer (0.1 M maleic acid, 0.15 M NaCl, pH 7.5) for 1 h at room temperature, followed by incubation in anti-DIG antibody (anti-digoxigenin-AP (Roche, 11093274910)) at 1:5,000 dilution in DIG blocking buffer overnight at 4 °C. The next day, the samples were washed in PBST for 10 min 3 times each at room temperature, then in 9.5T buffer (100 mM Tris-HCl pH 9.5, 50 mM MgCl 2 , 100 mM NaCl, 0.1% Tween 20) for 10 min 3 times each at room temperature. Hybridization signals were revealed by incubation in BCIP/NBT buffer (1 SIGMAFAST BCIP/NBT tablet (Sigma, B5655) in 10 ml H 2 O)) at 4 °C until brown-purplish colours were sufficiently dark. For this study, the colour development took 48 h. The samples were then washed in PBST twice for 10 min each. The samples were post-fixed in 4% PFA overnight at 4 °C, followed by washing in PBST twice for 10 min each, and then washed in methanol for 3 h at room temperature. The tissues were kept in PBS and imaged using SMZ1500 microscope (Nikon) under Ring Light System (Fibre-Lite). For cross-sections of stalks, the whole-mount sample was processed for cryo-section as described in 'Xenia regeneration, BrdU labelling and EdU pulse-chase'.

RNAscope ISH assay for LePin and Granulin 1 expression
To visualize RNA expression in endosymbiotic cells, we used the ultrasensitive RNAscope ISH approach (Advanced Cell Diagnostics (ACD)). LePin-or Granulin-1-specific oligonucleotide probes were ordered from ACD (see Supplementary Table 7 for further information). The fluorescent RNAscope assay was carried out by RNAscope Multiplex Fluorescent Reagent Kit v.2 (ACD) according to the manufacturer's protocol. The chromogenic assay was carried out by RNAscope 2.5 HD Duplex Detection Kit (ACD), according to manufacturer's protocol. Both assays used the cryo-section of the fixed Xenia polyp prepared according to the manufacturer's protocol.

Genome assembly
Sequencing data from Nanopore were used to initiate the genome assemble by Canu (v.1.7) 44 . The assembled genome was further polished with Illumina short reads by Nanopolish (v.0.9.2, https:// github.com/jts/nanopolish) with 5 cycles, which resulted in 1,482 high-quality contigs for the diploid genome. The diploid genome assembly was separated into haploid by HaploMerger2 45 . The haploid genome assembly was further subject to Hi-C assisted scaffolds by 3D de novo assembly pipeline, Juicer (v.1.5) 14 . By aligning all the Illumina genomic sequencing data with the assembled genome, we found 0.45% single nucleotide polymorphism (SNP) within the whole assembled genome of the Xenia sp.

Gene annotation
The funannotate genome annotation pipeline (v.1.3.3, https://github. com/nextgenusfs/funannotate) was used to annotate the Xenia sp. genome. In brief, transcriptome data were assembled by Trinity (v.2.6.6) 46 and used to generate the gene models based on the presence of mRNA by PASApipeline (v.2.3.2) 47 . These gene models were used as training sets to perform de novo gene prediction by AUGUSTUS (v.3.2.3) 48 and GeneMark-ES Suite (v.4.32) 49 . All gene models predicted by PASApipeline, AUGUSTUS and GeneMark were combined and subjected to EVidenceModeller to generate combined gene models 50 . The predicted genes were filtered out if more than 90% of the sequence overlapped with repeat elements as identified by RepeatMasker and RepeatModeler (http://www.repeatmasker.org). PASA was further used to add 3′ and 5′ untranslated region sequences to the remaining predicted genes. Pfam  52 were used to annotate the function of these gene models. Among all the predicted genes, 23,939 (82.5%) gene models were supported by transcriptome data because they have detectable reads (reads number >0). Among these models, 20,397 have read numbers >5.

Phylogeny tree analysis
We used OrthoFinder (v.2.2.7) to find orthologues from different species on the basis of protein sequences from 13 species listed Fig. 1d, and inferred the species tree 53,54 . In brief, 'orthofinder -S diamond -t 22 -M msa -f fasta_files' was used to generate the result. Diamond (v.0.9.21) was used for sequence search and OrthoFinder grouped 308,348 genes (83.8% of total) into 19,244 orthogroups. One thousand six hundred and one orthogroups, according to previously reported method 55 , with a minimum 10 species having single-copy genes, were used to infer the species tree. These orthogroups were subjected to multiple sequence alignment by MAFFT (v.7.407) and columns with more than eight gaps were trimmed. The trimmed alignment with 73.6% data occupancy (see Source Data for Fig. 1d) was used to infer the maximum likelihood unrooted species tree by FastTree (v.2.1.10) with the default configuration in OrthoFinder. This species tree was further rooted by the STRIDE algorithm, which has been demonstrated to correctly root the species tree spanning a wide range of time scales and taxonomic groups 56 .

Single-cell clustering and marker gene identification
The raw single-cell sequencing data were de-multiplexed and converted to FASTAQ format by Illumina bcl2fastaq (v.2.20.0) software. Cell Ranger (v.3.1.0, https://support.10xgenomics.com/ single-cell-gene-expression/software/overview/welcome) was used to de-multiplex samples, process barcodes and count gene expression. The sequence was aligned to the annotated Xenia sp. genome and only the confidently mapped and non-PCR duplicated reads were used to generate gene expression matrix for each library with 'cellranger count' command. The expression matrix of Cell-Ranger-identified cells from each library was read into R and further analysed with Seurat (v.3.0.2) 57 . Cells with UMI numbers less than 400 or mitochondria gene expression >0.2% of total reads were excluded for downstream analysis. To further remove outliers, we calculated the UMI number distribution detected per cell and removed cells in the top 1% quantile.
To remove batch effect and integrate data from different libraries, we applied the Seurat v.3 method for data integration 57 . For each dataset, we identified the top 1,000 genes with the highest dispersion. We used the top 1,000 genes in the non-regeneration sample as anchor features to identify anchors between different non-regeneration datasets. The first 20 dimensions were used to generate the integrated data. Dimensional reduction was carried out on the integrated data, and used for further clustering analysis. Clustering and marker gene identification in non-regeneration condition was further performed with Seurat v.3. The cell clusters in regeneration samples were identified with the label transfer method in Seurat v.3. All violin plots were generated using Seurat VlnPlot function.

Identification of Xenia sp. cells performing endosymbiosis with Symbiodiniaceae
The bulk transcriptome data of FACS-isolated alga-containing or alga-free Xenia cells were aligned to Xenia sp. genome by STAR (v.2.5.3a) 58 . Individual gene expression (reads per kilobase of transcript, per million mapped reads) for each sample were calculated by RSEM (v.1.3.0) 59 . The gene-expression levels of each bulk RNA-seq of FACS-isolated cells were compared with the gene-expression levels calculated using average UMI number for each gene in each cell cluster identified by scRNA-seq. The Pearson correlation coefficient was calculated for each comparison.

Pseudotime analysis
To infer the trajectory of endosymbiotic Xenia cells, we integrated scRNA-seq data of regenerating and non-regenerating samples using Seurat v.3. All cells belonging to the endosymbiotic cell cluster (cluster 16, total of 382 cells) were subjected to Monocle (v.2.10.1) 29 analyses. To find the variable genes among these cells for downstream analysis, we grouped these cells into three subclusters with Monocle cluster-Cells function (with default setting for most parameters, except for num_clusters = 4, which generated 3 clusters). Each of these three subclusters contains 247, 53 or 82 cells. The top 1,000 differentially expressed genes between these three subclusters were used as ordering genes to construct the trajectory by DDRTree algorithm. The differentially expressed genes along pseudotime were detected using the differentialGeneTest function in Monocle. The cell numbers in each of the five predicted endosymbiotic cell states are state 1 = 36, state 2 = 109, state 3 = 155, state 4 = 45 and state 5 = 37.

RNA velocity
RNA velocity estimation was carried out using the velocyto.R program (http://velocyto.org, v.0.6), according to the instructions 30 . In brief, velocyto used raw data of the regeneration sample to count the spliced (mRNA) and unspliced intron reads for each gene to generate a .loom file. This .loom file was loaded into R (v.3.6.1) using the read.loom. matrices function and used to generate the RNA velocity map. The RNA velocity map was projected into the t-SNE space that was identified by Seurat.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Data availability
We have uploaded all raw genomic, bulk RNA-seq and scRNA-seq data to NCBI (BioProject PRJNA548325). The genome files are available at http://cmo.carnegiescience.edu/data; we have also made the genome data interactive using UCSC genome browser, http://genome.ucsc. edu/cgi-bin/hgTracks?hubUrl=http://cmo.carnegiescience.edu/gb/ hub.txt&genome=xenSp1. We allow anyone interested to explore the predicted proteomes of Xenia and 14 other cnidarian using our blast server: http://c-moor.carnegiescience.edu:4567. All scRNA-seq analyses and results are available at GitHub: https://github.com/ciwemb/ endosymbiosis. Select intermediate RDS objects are available at: http:// cmo.carnegiescience.edu/data. We have worked to prototype a web portal to organize all the above links. This work-in-progress has a goal of making research findings, experimental protocols and computational data available to the scientific community. As the portal involves information beyond this study, we are still working with colleagues to best design it so that it will be easy to use and informative. The portal can be accessed at: http://cmo.carnegiescience.edu. Source Data are provided with this paper.

Code availability
Corresponding author(s): Minjie Hu, Chen-Ming Fan, Yixian Zheng Last updated by author(s): Apr 12, 2020 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see Authors & Referees and the Editorial Policy Checklist.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability Raw sequence data for this study is available in NCBI BioProject under accession PRJNA548325. Assembled genome and gene annotation are available at http:// cmo.carnegiescience.edu/data. The scRNA analysis code is available at https://github.com/ciwemb/endosymbiosis nature research | reporting summary

October 2018
Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection.

Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences
For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf

Life sciences study design
All studies must disclose on these points even when the disclosure is negative.

Sample size
No statistical methods were used to predetermine sample size. We followed standards in the biology field.
Data exclusions For single cell RNA-seq analysis, based on pre-established criteria for single-cells, in order to remove empty droplet, or droplet with potential dead cells or potential doublets, cells with UMI numbers less than 400 or mitochondria gene expression >0.2% were filtered out. To further remove outliers, we calculated the UMI number distribution detected per cell and removed cells in the top 1% quartile.

Replication
Each experiment was replicated with multiple independent animals. To draw a conclusion, at least two independent experiments were carried. All replicates were successful.
Randomization Xenia colonies or polyps were randomly chosen from the aquarium tank Blinding Quantification of LePin signal was blinded by de-identifying samples.

Reporting for specific materials, systems and methods
We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response.

Validation
The BrdU antibody was validated by a lot of studies listed in the manufactory's website: https://www.thermofisher.com/ antibody/product/BrdU-Antibody-clone-ZBU30-Monoclonal/03-3900. It has been applied in IF, IHC, FACS in Chemical, Chicken, Mouse, Rabbit and Rat. We validated it by the lack of staining when BrdU was not added into the sample.

Animals and other organisms
Policy information about studies involving animals; ARRIVE guidelines recommended for reporting animal research Laboratory animals Xenia sp. was cultured in laboratory aquarium tank. We can not yet tell their age and sex.

Wild animals
The Xenia sp. used in this study was originally from the wild, but we obtained it from an aquarium shop in Baltimore.

Field-collected samples
The study didn't involve samples collected from field

Ethics oversight
The study of Xenia or some other cnidaria does not yet have ethical oversight.
Note that full information on the approval of the study protocol must also be provided in the manuscript.