Loss of neurogenesis in Hydra leads to compensatory regulation of neurogenic and neurotransmission genes in epithelial cells

Hydra continuously differentiates a sophisticated nervous system made of mechanosensory cells (nematocytes) and sensory–motor and ganglionic neurons from interstitial stem cells. However, this dynamic adult neurogenesis is dispensable for morphogenesis. Indeed animals depleted of their interstitial stem cells and interstitial progenitors lose their active behaviours but maintain their developmental fitness, and regenerate and bud when force-fed. To characterize the impact of the loss of neurogenesis in Hydra, we first performed transcriptomic profiling at five positions along the body axis. We found neurogenic genes predominantly expressed along the central body column, which contains stem cells and progenitors, and neurotransmission genes predominantly expressed at the extremities, where the nervous system is dense. Next, we performed transcriptomics on animals depleted of their interstitial cells by hydroxyurea, colchicine or heat-shock treatment. By crossing these results with cell-type-specific transcriptomics, we identified epithelial genes up-regulated upon loss of neurogenesis: transcription factors (Dlx, Dlx1, DMBX1/Manacle, Ets1, Gli3, KLF11, LMX1A, ZNF436, Shox1), epitheliopeptides (Arminins, PW peptide), neurosignalling components (CAMK1D, DDCl2, Inx1), ligand-ion channel receptors (CHRNA1, NaC7), G-Protein Coupled Receptors and FMRFRL. Hence epitheliomuscular cells seemingly enhance their sensing ability when neurogenesis is compromised. This unsuspected plasticity might reflect the extended multifunctionality of epithelial-like cells in early eumetazoan evolution.


Supplt-1: Detailed methods describing the RNA-seq assembly of the different transcriptomes and the quantification of the abundances of transcripts.
An extended description of the procedures relating to RNA-seq experiments are provided in sections S1.1 to S1.5. Further details concerning the command lines and program versions are provided in section S1.6. S1.1 De-novo assembly of the Hv_Jussy transcriptome used for spatial gene profiling and drug / heat-shock treatments For spatial gene profiling, 25 animals from the Hv Jussy strain were dissected for each replicate as depicted in Fig. 1B, each body slice being ~250 µm thick, and reads were mapped to a Hv_Jussy de-novo transcriptome. All tissue samples were immediately placed in RNALater (Qiagen) and total RNAs were extracted the same day (RNAeasy mini kit, Qiagen). All conditions were collected in biological triplicates over different weeks. An average of ~30 (± 2.5 SD) million reads per sample was sequenced. Before de-novo assembly, sequencing adapters and trans-splice leaders (Stover and Steele 2001;Derelle, et al. 2010) were removed using cutadapt (Martin 2011), reads were corrected using SEECER (Le, et al. 2013) and cd-hit-454 (Li and Godzik 2006). Finally, digital normalization was performed using two rounds of the Trinity normalization tool. The resulting dataset was assembled using Trinity (Grabherr, et al. 2011;Haas, et al. 2013) with default options and Velvet/Oases (Zerbino and Birney 2008;Schulz, et al. 2012) using a kmer value of 35 and final coverage of 5. As post-processing, the velvet (430'345 sequences) and trinity (773'456) assemblies were pooled and, for each sequence, the longest open reading frame (ORF) of more than 100 amino acids encoded by the forward strand was extracted using EMBOSS tools (Rice, et al. 2000). Near duplicates sequences were removed using a 98% similarity threshold on both amino acids and nucleotides Cd-hit and Cd-hit-est (Li and Godzik 2006) yielding 131'430 sequences. Finally, only sequences aligned by blast (Camacho, et al. 2009) that match the Hm-105 strain genome (Hydra vulgaris group) (Chapman, et al. 2010) with at least 90% similarity and E-values lower than 10 -10 were selected for further steps (107'422 sequences).
Discriminating between sequences originating from the same gene (e.g. alleles or natural variants) and close paralogs is a difficult task in absence of a comprehensively assembled genome. In order to facilitates this aspect in downstream analyses, remaining sequences were realigned to the Hm-105 strain genome (Chapman, et al. 2010) using GMAP (Wu and Watanabe 2005) and assigned to genomic loci, considering here that one single base overlap was sufficient for an assignation into the same locus. A final round of clustering was performed per locus with a 95% sequence identity threshold on deduced protein sequences. The assembly resulted in 79'773 sequences (51'633 and 28'140 from the Trinity and Velvet/Oases assembly, respectively), arising from 25'878 loci.

S1.2 RNA-seq analyses upon loss of neurogenesis using the thermosensitive Hv_Sf1 strain (HU/HS/Col treatments)
For a given condition the central body columns of 35-40 Hv_Sf1 polyps were dissected and pooled together, from untreated (3, 4, 6 and 10 starvation days), HU-treated (0, 1, 3, 7 days post-HU), HS-treated (7 days post-HS), Col-treated (10 days post-Col) animals (Fig. 4A). Each condition was sampled in 3 or 4 biological replicates, representing 37 samples in total. Libraries were prepared with the Low Sample TruSeq total RNA preparation protocol from Illumina (San Diego, CA, US). Pools of 4 or 5 multiplexed libraries were loaded per lane of a HiSeq2500 sequencer (Illumina) and single-end sequenced up to 100 nt. After trimming adapters and trans-spliced leader using cutadapt, reads from control and treated Hv_Sf1 were mapped to the Hv_Jussy transcriptome as the observed sequence similarities are nearly 100% between these two strains. S1.3 Stem-cell specific RNA-seq using transgenic Hv_AEP strains  For cell-type specific transcriptomics, we used reads obtained from FACS-sorted cells from body columns of transgenic animals obtained from Hv_AEP strains expressing actin::eGFP (strain called here Endo-GFP) ), or actin::eGFP (Ecto-GFP) (Wittlieb, et al. 2006) in endodermal and ectodermal epithelial cells respectively, or Cnnos1::GFP (Cnnos1-GFP) in i-cells (Hemmrich, et al. 2012). These strains were kindly provided to us by Thomas Bosch (Kiel, Germany). Four biological replicates were prepared per condition. The body columns from 300-400 Hv_AEP transgenic polyps were dissociated with pronase (6 mg/ml) in Gierer dissociation medium (Gierer, et al. 1972). GFP positive cells from the Cnnos1-GFP strain were sorted with a FACS Area (Beckton-Dickinson), GFP positive cells from the Ecto-GFP and Endo-GFP strains (Hemmrich, et al. 2012;Buzgariu, et al. 2014) with a MoFlow Astrios (Beckman Coulter). The sorted cells (3.10 5 -6.10 5 ) were centrifuged, resuspended and kept in RNACell protect (Qiagen) until RNA extraction with RNeasy Plus kit (Qiagen). In addition to these FACS-sorted samples, two samples were prepared from unsorted body columns. A de-novo transcriptome was assembled from the 12 FACS-derived samples using Trinity after adapter and trans-spliced leaders removal (cutadapt) and in-silico reads normalization. It yielded 61'501 transcripts, arising from 44'306 putative loci (according to trinity naming scheme). Reads trimmed with cutadapt from the different Hv_AEP transgenic strains were mapped to the Hv_AEP transcriptome for the quantification steps (see section S1.6).

S1.4 Quantification of transcript levels and other data analysis
Mapping steps were performed separately for each library using Bowtie2 (Langmead and Salzberg 2012) with strand specificity and otherwise default options. Count tables were produced by counting the total number of mapped reads aligning to each reference sequence. Inter-sample library normalizations and statistical analyses were performed using DESeq2 version 1.6.3 (Love, et al. 2014) with default options. Most graphs were produced using the ggplot2 (Wickham 2009) and ggtern packages (www.ggtern.com). When biological replicates required to be averaged (as on ternary plots), geometric means of normalized read counts were used.

S1.6 Detail of the command lines and the versions of the programs used
De-novo transcriptome assembly and quantification of transcripts levels (Hv_Jussy strain, spatial gene profiling) Steps 1 to 3 were ran in parallel for every libraries. Libraries were then pooled for the remaining steps.

Altogether re-normalization (Trinity normalisation tool version r20131110)
Input file: all_libraries_normalized.fa Output file: all_libraries.fa_normalized_K25_C50_pctSD200.fa ./util/normalize_by_kmer_coverage.pl --seqType fa --JM 500G --single all_libraries.fa --max_cov 50 --SS_lib_type F --output altogether_normalization The output has been renamed norm_reads.fa for clarity.  All open reading frames of more than 300 nucleotides are obtained (possibly multiple ORFs per initial sequences). Getorf also annotate the header of each sequence with the coordinates of the detected ORF. If the ORF is encoded in the 3'  5' direction, a (REVERSE SENSE) flag is added to the header.

Blast 2.2.29+ alignment to Hm-105 genome using a 90% threshold to remove contaminants
The Hydra genome has been produced in two flavours. We use here the version referred in Chapman et al. 2010 as the "RP" assembly.

Briefly:
a. When a sequence matches multiple locations in the genome, as simple score consisting of the number of matches minus the number of mismatches is used to rank alignments (field1 -field2 in PSL format) and only the best alignment for a given transcript is kept. The name of the genomic contig is retrieved as well as genomic coordinates where the match occurs and the transcript name. The file is then sorted. b. Within matches to a genomic contig, overlapping transcripts are identified using the lowest_coordinate and highest_coordinate fields and assigned identical locus numbers.
Finally, the original full-length RNA sequences (either from velvet or trinity assemblies produced as described in subsections 6a and 6b) were retrieved and the transcriptome obtained was named hydra_transcriptome_jussy.fa

Quantification
After indexing of the transcriptome produced in the previous step, bowtie2 v2.2.1 was used using the command line, for each library in parallel using the output from cutadapt (see above).
bowtie2 --nofw -x hydra_transcriptome_jussy -U lib001_cutadapt.fastq -S 001.sam The total number of reported reads mapping to all contigs were counted and reported in a count table.

De-novo transcriptome assembly and quantification of transcript levels (Hv_AEP strain, stem-cell specific quantifications)
READS PRE-PROCESSING
The file tspl.conf used here is strictly identical to the one used for Hv_Jussy containing known trans-spliced leader and illumina 3' TruSeq adapter (see above)

Quantification
After indexing of the transcriptome produced in the previous step, bowtie2 v2.2.4 was used using the command line, for each library in parallel using the output from cutadapt (see above).
bowtie2 --nofw -x aep_transcriptome --U ./lib001_AEP_cutadapt.fastq --S lib001.sam The total number of reported reads mapping to all contigs were counted and reported in a count table.

Quantification of transcript levels (Hv_Sf1 strain, drugs and heat-shock treatments)
READS PRE-PROCESSING With the file tspl.conf being strictly identical to the one used for Hv_Jussy containing known trans-spliced leader and illumina 3' TruSeq adapter (see above)

Quantification
After indexing of the transcriptome produced in the previous step, bowtie2 v2.2.4 was used using the command line, for each library in parallel using the output from cutadapt (see above). The reads from Hv_Sf1 animals were mapped to the reference Hv_Jussy transcriptome as described above.
bowtie2 --nofw -x hydra_transcriptome_jussy -U lib001_drugs_sf1_cutadapt.fastq -S lib001.sam The total number of reported reads mapping to all contigs were counted and reported in a count table.

Hydra strains (Hv_Basel, Hv_Jussy, Hv_AEP)
Supplt-2: Expression patterns of six NG Dlx,DMBX1,Ets1, and nine NT (Arminin 1b-l, Arminin 01798, FMRFRL, GPRA0.6, Inx1, Inx8, MEP1B, PAH, PW peptide) genes analysed by ISH in three H. vulgaris strains, starved for either three days (Hv_Sf1) or 10 days (Hv_Basel, Hv_AEP). See the corresponding RNA-seq profiles in Fig. 5F, 5G and in Supplemental Table-1. All tested genes show a similar spatial expression pattern and an identical cell-type signature in the three H.vulgaris strains: the neuropeptide Hym-355 in the apical and basal nerve nets, the Dickkopf related gene Dlp-1 gradually in gland cells along the body column, the metallloendopeptidase MEBP1 in the gland cells homogenously along body column, the homeobox gene prdl-b in dividing nematoblasts, the innexin Inx8 in the i-cells. Concerning the ectodermal epithelial cells, Inx1, PW peptide and Dlx are homogenously expressed along the axis ,while GPRA0.6 transcripts are detected at the base of tentacles. The endodermal epithelial cells express the genes encoding the receptor FMRFRl, the antimicrobial peptides Arminin 1b-l and Arminin 01798. In contrast, the genes encoding the transcription factors Ets1, DMBX1 and the Phenyl Alanine Hydroxylase -PAH are expressed in the epithelial cells from both the ectodermis and the gastrodemis. Note the bipolar expression pattern of Ets1, identified only at the base of tentacles and in the foot region of all three strains. Picturing were performed with an SZX10 Olympus steromicroscope (low magnification) and a Zeiss Axioplan2 microscope (high magnification).

Supplt-3: Impact of HU treatment on the maintenance of the nervous system in homeostatic tissues and on the regeneration of apical and basal nervous systems in Hydra
Supplt-3: In each animal, the nervous system can be observed at one pole in homeostatic (nonregenerating) condition (arrows), and at the other pole in regenerating context (arrowheads). Note that HU given before mid-gastric bisection (HU 4d) dramatically affects neurogenesis in regenerating tissues either basal (D) or apical (E) without affecting the nervous system in homeostatic tissues (apical in D and B, basal in E and C). By contrast HU given immediately after bisection (HU 8d) only slightly affects neurogenesis in regenerating tissues, either basal (F) or apical (G). Finally HU given continuously (HU 12d) affects both neurogenesis and patterning processes in regenerating tissues (H, I). Also note the extended basal nervous system in I. dpa: days post-amputation.

Supplt-4: Modulations of spontaneous contractile activity in HU-treated polyps.
Supplt-4A: Control polyps and 4d HU-treated polyps (as indicated in Fig. 2A) were bisected at midgastric position and let to regenerate for nine days in Hydra Medium. The contractile activity of intact, head-and foot-regenerated animals was then recorded as follows: 10 animals were transferred in 1 ml HM, left for 5 minutes and recorded every 0.5 second for 15 minutes with an SZX10 Olympus stereomicroscope equipped with a DP73 Olympus camera. The number of spontaneous contraction bursts, which comprise several coordinated contractions that reduce the animals to a small ball, was counted and the contractile activity was expressed as number of contraction burst per hydra per 15 minutes. Untreated animals that regenerated their head or foot exhibit a contractile activity similar to that observed in intact polyps. As expected the HUtreated animals that regenerate heads or feet lacking a nervous system contract significantly less frequently than control animals having regenerated the same body part. However their contractile activity still persists, significantly more reduced when the apical nervous system did not regenerate properly when compared to the impact of a limited basal nervous system regeneration. A Mann-Whitney U statistical test was used to establish the p-values. See the corresponding movies as indicated below. previously performed ISH on whole polyps and characterized genes as specifically expressed either in the nematocyte lineage (48, coral color), or in the gland cells (18 transcripts, blue color) or in nerve cells (8 transcripts, green color). We selected these genes and tested their RNA-seq expression levels in each celltype specific transcriptomes (Ecto-GFP, Endo-GFP, Cnnos1-GFP) as well as. in a RNA-seq transcriptome established on whole unsorted gastric columns. The results were normalized for each gene as such as a relative expression level of 100% is obtained in the unsorted gastric column samples. The analysis shows that (i) transcripts attributed to nematocyte genes are mainly found in the Ecto-GFP fraction, (ii) transcripts attributed to gland cell genes are found in the Endo-GFP fraction, and (iii) transcripts attributed to nerve cell genes are found in both the Ecto-GFP and the Endo-GFP fractions. All these results are consistent with the location of these different cell types in the epidermis and the gastrodermis. Except for the nematocyte transcripts, these contaminations do not exceed 20% of the expression level measured in the whole body column. Finally the Cnnos1-GFP fraction appears as not contaminated.

Supplt-9: Expression analysis of the gland cell genes found up-regulated after the loss of neurogenesis
Supplt-9: (A) RNA-seq profiles of the gland cell genes identified as upregulated upon the loss of neurogenesis in drug-(HU, colchicine) or HS-treated animals (see Fig.5B-5E). The graphs show for each gene the highly variable spatial distribution of transcript levels (blue dots), the almost exclusive detection of transcripts in the Endo-GFP fraction (green dots) and the increased expression levels 11 days after HU or HS treatments but the dramatically low level of transcripts after colchicine (red dots). Y axis indicates the number of k-reads. (B) WM-ISH expression pattern of four gland cell genes (Dkk1/2/4-A, Dkk1/2/4-C or MEP1B,Kazal1) in Hv_Basel polyps either untreated and starved for 10 days (control) or exposed to HU (as in Fig. 4A) and fixed 7 days later (HU). As previously reported by Augustin et al. (2006), Dkk1/2/4-C exhibits a graded expression in the gastrodermis, whereas Dkk1/2/4-A, MEP1B and Kazal1 display a homogenous expression all along the body axis except the head and the foot regions. Note that upon the loss of neurogenesis, these genes are still expressed by the gland cells, which are not eliminated by the HU treatment. Quantitative assessments cannot be made from WMISH patterns as the permeability of the tissue layers decreases after HU exposure.

Supplt-10: Quantitative analysis and cell-type expression of the candidates genes up-regulated upon loss of neurogenesis.
Supplt-10: (A) Q-PCR validation of the up-regulation of two NG and seven NT epithelial genes detected after the loss of neurogenesis. Real-time Q-PCR was performed on mRNA extracted seven days after HU withdrawal, corresponding to 10 days of starvation for control animals. The relative expression levels of the candidate genes were normalized against three housekeeping genes (EF-1γ; TBP1, YWHAZ). On this graph we represented the relative fold change in expression upon HU treatment assessed by both Q-PCR (pinkred dots) and RNA-seq (blue dots). The obtained data were multiplied by a factor as such as the average values of control conditions are equal to "1" for each gene. The horizontal dashed line represents the twofold threshold. Note the good correlation between the RNA-seq and the qPCR quantifications.