Transcriptomic Profiling of Zebrafish Hair Cells Using RiboTag

Introduction

Hearing loss is a genetically heterogeneous disorder, with mutations in over 150 genes estimated to underlie genetic non-syndromic hearing deficits (Van Camp and Smith, 2017). Of these, a large proportion affect genes that are preferentially expressed in the sensory cells of the inner ear, namely the mechanosensory hair cells (HCs) (Elkon et al., 2015). Zebrafish have served as an excellent model system for functional analysis of genes in HC function, as they possess HCs both in the inner ear as well as within an external lateral line system, are easy to manipulate genetically, and generate large numbers of progeny within a short gestational period (Nicolson, 2005, 2017; Erickson and Nicolson, 2015). However, the study of the molecular changes induced by manipulation of HC-expressed genes in zebrafish has been limited due to a paucity of models that allow cell type-specific molecular analysis of changes in gene expression. Specifically, across vertebrate species, the auditory, vestibular, and lateral line sensory organs are comprised of a variety of cell types, of which HCs make up only a small percentage (Hertzano and Elkon, 2012; Jiang et al., 2014; Matern et al., 2017). Therefore, due to their relative scarcity, cell type-specific approaches such as manual cell sorting, fluorescence activated cell sorting (FACS), or single cell RNA-Seq (scRNA-Seq) are necessary to analyze gene expression in HCs. These methods have been used in both mice and zebrafish to analyze HC gene expression changes that occur in mutant animals, during development and regeneration, or after exposure to noise or ototoxic drugs (McDermott et al., 2007; Hertzano and Elkon, 2012; Jiang et al., 2014; Steiner et al., 2014; Burns et al., 2015; Elkon et al., 2015; Scheffer et al., 2015). However, both scRNA-Seq and cell sorting-based techniques require dissociation of tissues to obtain a single cell suspension. Tissue dissociation can induce significant cellular stress due to loss of lateral inhibition and cell-cell contact, and combined with the prolonged time associated with tissue processing, may lead to confounding changes in gene expression (Sanz et al., 2009; Gay et al., 2013, 2014).

To avoid dissociation-induced molecular changes, recent studies have developed techniques in both mice and zebrafish models to extract RNA from specific cell types within intact organs (Heiman et al., 2009; Sanz et al., 2009; Gay et al., 2013; Tryon et al., 2013; Erickson and Nicolson, 2015; Roh et al., 2017). These approaches rely on pulldown of RNA from a cell type of interest via tagged ribosomes, or directly tagged RNA. As an example, the RiboTag mouse model expresses a component of the 60S subunit of the ribosome with a C-terminal hemagglutinin tag (RPL22-HA) (Sanz et al., 2009). Using this model, Cre-induced expression of RPL22-HA can be used to capture actively translated RNA from a cell type of interest via immunoprecipitation. The BACarray and NuTRAP mice models work in a similar way to RiboTag, however ribosomes are tagged with a green fluorescent protein (GFP). Additionally, the NuTRAP mice co-express nuclear tagging proteins that allow for concomitant epigenetic profiling of a Cre expressing cell type. In 2013, Tryon et al. also introduced several RiboTag/TRAP models to study cell type-specific gene expression in zebrafish (Tryon et al., 2013). Similar to the BACarray and NuTRAP mouse models, these zebrafish models rely on tagging the 60S ribosomal subunit through non-inducible tissue-specific expression of an rpl10a-GFP fusion gene. Collectively, in mouse and zebrafish, RiboTag, BACarray, and NuTRAP can be used for cell type-specific isolation of RNA via pulldown of labeled ribosomes. In these models, the immunoprecipitated RNA is enriched for actively translated transcripts, and is referred to as the “translatome” rather than the whole cellular transcriptome.

An alternative method to ribosomal pulldown is to isolate cell type-specific RNA from a complex tissue environment by thiouracil (TU) tagging. This method relies on tissue specific expression of uracil phosphoribosyltransferase (UPRT), an enzyme capable of integrating 4-thiouracil into newly synthesized RNA, which can then be affinity purified. Both mice and zebrafish TU-tagging models have been developed, and rather than obtaining only actively translated RNA as in the RiboTag/TRAP models, this method allows for capture of all newly synthesized RNA within a cell type of interest (i.e., the transcriptome) after the application of 4-thiouracil (Gay et al., 2013, 2014; Erickson and Nicolson, 2015). As with the ribosomal pulldown techniques, affinity purified RNA is enriched for the cell type of interest, and also contains some RNA from other tissues. In zebrafish, TU-tagging has been used to identify HC expressed transcripts. However, this model was only able to identify a small number of genes, and only those genes with very high expression levels in HCs were found to be significantly enriched in the immunoprecipitated RNA compared to input (Erickson and Nicolson, 2015).

To more effectively isolate RNA from zebrafish HCs, we have developed a transgenic zebrafish RiboTag model to evaluate the translatome of zebrafish inner ear and lateral line HCs; the Tg(myo6b:GFP-2A-rpl10a-3xHA) zebrafish [from here on referred to as Tg(myo6b:RiboTag)]. This model carries a construct to drive expression of both an HA-tagged Rpl10a protein and GFP in HCs under control of the HC-specific myo6b promoter (Seiler et al., 2004; Obholzer et al., 2008). To our knowledge, this is the first zebrafish model to allow for HC-specific gene expression analysis via two methods: (1) tissue dissociation and cell sorting based on GFP expression, and (2) immunoprecipitation of HA-tagged ribosomes to enrich for HC expressed transcripts. We use both RT-qPCR and RNA-Seq to analyze gene expression in our Tg(myo6b:RiboTag) model. We show that immunoprecipitated RNA from our Tg(myo6b:RiboTag) model is significantly enriched for known HC expressed transcripts, indicating that this model is effective in enriching for the HC translatome. Additionally, a comparison of our translatome dataset with a previously published zebrafish HC transcriptome dataset (generated using sorted HCs) shows that similar gene expression results can be obtained using the Tg(myo6b:RiboTag) model without cell sorting. Finally, we use the Tg(myo6b:RiboTag) model to identify novel HC expressed transcripts, and demonstrate that RiboTag immunoprecipitation helps to avoid gene expression changes that are induced by dissociation. Overall, we believe that this model will be highly useful for studying the normal development and function of zebrafish HCs, as well as changes to HC gene expression in response to different conditions.

Materials and Methods

Zebrafish Husbandry

Zebrafish were grown at 28°C in E3 embryo media (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂, and 0.33 mM MgSO₄) using standard methods. Work performed at the National Institute of Health was approved by the NIH Animal Use Committee under animal study protocol #1362-13. At the University of Maryland School of Medicine (UMSOM), all procedures involving animals were carried out in accordance with the NIH Guide for the Care and Use of Laboratory Animals and have been approved by the Institutional Animal Care and Use Committee at the University of Maryland, Baltimore (protocol numbers 0514001 and 1116003).

Vector Construction and Generation of the Tg(myo6b:GFP-2A-rpl10a-3xHA) Zebrafish

Plasmid construction was based on the Tol2/Gateway zebrafish kit (Kwan et al., 2007). The pME-GFP-2A-rpl10a-3xHA was a generous gift from Dr. Brant Weinstein at the NIH. The p5E-pmyo6b entry clone used to drive expression in hair cells has been described previously (Kindt et al., 2012). These two clones were used along with the tol2 kit gateway clones p3E-polyA (#302) and pDest (#394) to create the myo6b:GFP-2A-rpl10a-3xHA expression construct. To generate a stable transgenic fish line, plasmid DNA at 50 ng/μL and tol2 transposase mRNA at 20 ng/μL were injected into zebrafish embryos as previously described to create Tg(myo6b:GFP-2A-rpl10a-3xHA)^idc10(Kwan et al., 2007). Each line was grown to the F1 generation and outcrossed to confirm single copy integration.

Immunohistochemistry and Confocal Imaging

Immunohistochemistry was performed on whole-mount larvae. Wildtype or transgenic larvae were fixed with 4% paraformaldehyde in phosphate buffered saline (PBS) for 4 h at 4°C. After 5 × 5 min washes in PBS, followed by a 5 min wash in H₂O, larvae were permeabilized with ice cold acetone (at −20°C) for 5 min. Larvae were then washed in H₂O for 5 min, followed by a 5 × 5 min washes in PBS, and then blocked overnight with PBS containing 2% goat serum and 1% bovine serum albumin (BSA). A primary rat anti-HA antibody (Roche) was diluted at 1:750 in PBS containing 1% BSA, and larvae were incubated in the solution for 4 h at room temperature. After 5 × 5 min washes in PBS to remove the primary antibody, an Alexa 568 secondary antibody diluted at 1:1,000 (Life Technologies) was added in PBS containing 1% BSA and incubated overnight at 4°C. After 5 × 5 min washes in PBS to remove the secondary antibody, larvae were rinsed in H₂O and mounted in Prolong Gold (Life Technologies). Fixed samples were imaged on an inverted Zeiss LSM 780 laser-scanning confocal microscope with a 63 × 1.4 NA oil objective lens. Excitation wavelengths of 488 and 546 nm were used to excite GFP and Alexa 568, respectively.

Tg(myo6b:RiboTag) Translatome Immunoprecipitation

The Tg(myo6b:RiboTag) immunoprecipitation protocol was modified from the RiboTag immunoprecipitation protocol described in Sanz et al. (2009). Briefly, GFP⁺ Tg(myo6b:RiboTag) zebrafish larvae were euthanized at 5 days post fertilization (dpf) using MS-222/Tricaine and rinsed with system water. Groups of fifty larvae were then either flash frozen and stored at −80°C or used immediately for immunoprecipitation. Larvae were resuspended and homogenized in 1 mL of supplemented homogenization buffer (50 mM Tris-HCl pH.7, 100 mM KCl, 12 mM MgCl₂, 1% Nonidet P-40, 1 mM 1,4-Dithiothreitol, 1X protease inhibitor cocktail, 200 U/mL RNAsin, 100 μg/mL cycloheximide, 1 mg/mL heparin) by douncing on ice. Homogenates were spun at 10,000 g for 10 min at 4°C to remove particulates, and a small sample of clear supernatant was reserved for total RNA isolation (input control, IN). Remaining supernatant was incubated with 5 μg HA antibody (BioLegend) at 4°C under gentle rotation for 4–6 h. After incubation, supernatants were incubated with 300 μL of rinsed Invitrogen Dynabeads Protein G magnetic beads (Thermo Fisher) overnight, rotating. Following incubation, bound beads were rinsed three times with 800 μL high salt buffer (50 mM Tris-HCl pH.7, 300 mM KCl, 12 mM MgCl₂, 1% Nonidet P-40, 1 mM 1,4-Dithiothreitol, 100 μg/mL cycloheximide) at 4°C for 10 min, rotating. After washing, 350 μL of buffer RLT from the RNeasy Plus Micro kit (Qiagen) was added to the bound beads or reserved input sample and vortexed for 30 s to dissociate bound RNA. RNA was then extracted according to manufacturer's instructions (using 16 μL of nuclease free water for elution) and stored at −80°C. RNA quality and concentration was assessed using the Agilent RNA Pico kit (Agilent Technologies) at the UMSOM Genomics Core Facility.

Fluorescence Activated Cell Sorting

Zebrafish dissociation and FACS were performed as previously described (Elkon et al., 2015). Briefly, approximately 150 Tg(myo6b:RiboTag) larvae per replicate were euthanized at 5 dpf and dissociated in cold trypsin-EDTA solution (0.5 g/L trypsin, 0.2 g/L EDTA, Sigma) by trituration with a p1000 pipette tip on ice for 20 min. Dissociation was then halted by adding HBSS supplemented with 10% fetal bovine serum (FBS) and 100 μg/mL DNaseI. Cells were filtered through a 70 μm cell strainer (Fischer Scientific) and pelleted by centrifugation at 2,000 rpm for 10 min at 4°C. Cells were washed once in HBSS, resuspended in HBSS supplemented with 10% FBS, and filtered into glass tubes through a 35 μm cell strainer (Falcon). Flow cytometry analyses were performed at the University of Maryland Marlene and Stewart Greenebaum Comprehensive Cancer Center Flow Cytometry Shared Service. Samples of GFP positive and negative cells (HCs and the rest of the fish) were collected using a BD FACSAria II (BD Biosciences), and a small aliquot of each sorted population was re-analyzed to determine cell purity. RNA was extracted from sorted cells using Trizol LS Reagent (Thermo Fischer Scientific) and the Direct-zol™ RNA MiniPrep Plus (Zymo Research).

RT-qPCR

RT-qPCR was performed as described previously with minor modifications (Matern et al., 2017). RNA from 5 dpf Tg(myo6b:RiboTag) zebrafish input and IP samples, or sorted cells, was reverse-transcribed using the Maxima First Strand cDNA Synthesis Kit (Thermo Fisher Scientific), and qPCR was performed using the Maxima SYBR Green/ROX qPCR Master Mix (Thermo Fisher Scientific). To account for low levels of RNA obtained from sorted cells, a preamplification step using PerfeCTa PreAmp Supermix (Quantbio) was added for the RNA-Seq validation and immediate early gene expression experiments. Expression values were normalized to actb1 expression (see Supplementary Table 4 for primer sequences).

RNA Sequencing and Informatics

RNA from 5 dpf Tg(myo6b:RiboTag) zebrafish IN and IP samples was submitted in biological quadruplicates for RNA-Seq at the UMSOM Institute for Genome Sciences. Only samples with RNA integrity numbers (RIN) > 8 were used for sequencing. Libraries were prepared from 25 ng of RNA using the TruSeq RNA Sample Prep kit (Illumina) per manufacturer's instructions, with the exception of an additional PCR cycle. Samples were sequenced on an Illumina HiSeq 4000 with a 75 bp paired-end read configuration. Between 100 and 150 million reads were obtained for each sample, and reads were aligned to the zebrafish genome (Danio rerio.GRCz10) using TopHat version 2.0.8 (maximum number of mismatches = 2; segment length = 30; maximum multi-hits per read = 25; maximum intron length = 50,000) (Kim et al., 2013). The number of reads that aligned to the predicted coding regions were determined using HTSeq (Anders et al., 2015), and only genes with CPM (reads count per transcripts per million mapped reads) values > 0.01 in all IN and IP replicates were called as expressed (17,164 genes). One IP sample had a high intergenic content suggestive of DNA contamination and was excluded from the analysis. See Supplementary Table 1 for alignment statistics. Significant differential expression was assessed using DEseq (Anders and Huber, 2010). RNA-Seq data were submitted to the Gene Expression Omnibus database (GEO accession GSE102861), as well as the gEAR Portal (UMgEAR.org). Gene ontology of highly enriched and depleted gene sets was performed using the Gene Ontology (GO) database (http://www.geneontology.org) (Harris et al., 2004). For HC enriched and depleted gene sets, top branches of GO terms are shown. Anatomical structure enrichment was performed using the Zebrafish Expression Ontology of Gene Sets (ZEOGS) tool (Prykhozhij et al., 2013), with a corrected p-value cutoff of 0.10. Sorted HC dataset transcript IDs from Steiner et al. were converted to gene IDs using bioDBnet (https://biodbnet-abcc.ncifcrf.gov) (Mudunuri et al., 2009).

Results

Generation of Tg(myo6b:GFP-2A-rpl10a-3xHA) Zebrafish to Create the Tg(myo6b:RiboTag) Model

To create a RiboTag zebrafish model that would allow for enrichment of the inner ear and lateral line HC translatome, we utilized the rpl10a gene, which has previously been used to effectively immunoprecipitate ribosomes in zebrafish via a GFP tag (Tryon et al., 2013). However, instead of an Rpl10a-GFP fusion protein, here we use an Rpl10a-3xHA fusion protein for the ribosomal pulldown. To drive Rpl10a-3xHA expression in HCs we utilized the HC-promoter myo6b, which has been shown to specifically drive expression of downstream genes in zebrafish inner ear and lateral line HCs (Obholzer et al., 2008). The construct is engineered to express Rpl10a-3xHA together with GFP bi-cistronically using the viral P2A peptide, resulting in the myo6b:GFP-2A-rpl10a-3xHA construct (Figure 1A). This construct was injected into zebrafish embryos at the one-cell stage to create a stable Tg(myo6b:GFP-2A-rpl10a-3xHA) transgenic line, referred to here as our Tg(myo6b:RiboTag) model. To confirm Rpl10a-3xHA was localizing properly in HCs, we immunostained our Tg(myo6b:RiboTag) zebrafish with an anti-HA antibody at 5 or 6 dpf, when HCs in the inner ear and along the lateral line have developed. As expected, we observed both GFP and HA signal localizing to the cytosol of inner ear and lateral line HCs (Figures 1B,C). Consistent with previously observed Rpl10a-GFP localization in xef1α>TRAP zebrafish (Tryon et al., 2013), the HA staining in the HCs of the Tg(myo6b:RiboTag) zebrafish also shows nucleolar localization (Figure 1C, white arrows), indicating that our model is localizing Rpl10a-3xHA properly.

Figure 1. The Tg(myo6b:RiboTag) zebrafish. (A) Schematic representation of the GFP-2A-rpl10a-3xHA construct driven by the HC-specific myo6b promoter. (B) Representative image of a live Tg(myo6b:RiboTag) zebrafish at 5 dpf showing GFP expression in the inner ear (IE) and lateral line neuromasts (N). (C) Immunohistochemistry using an HA antibody showing that Rpl10a-HA expression overlaps with GFP expression specifically within the HCs of the inner ear and neuromasts. White arrows denote nucleolar Rpl10a-HA staining.

Ribosome Immunoprecipitation Isolates the HC Translatome From Tg(myo6b:RiboTag) Zebrafish

In order to test whether the Tg(myo6b:RiboTag) zebrafish could be used to analyze gene expression of the inner ear and lateral line HCs, we next adapted the RiboTag immunoprecipitation protocol described in Sanz et al. (2009) to capture the HC translatome at 5 dpf (Figure 2A) using fresh or frozen embryos. We chose 5 dpf as the time point for our experiment, as this time point allows us to compare our results to other previously published HC gene expression datasets also generated using 5 dpf larvae. This technique yields ~50 ng of immunoprecipitated RNA (IP, average RNA concentration = 3.5 ± 2.3 ng/μL in 16 μL elution volume) and 660 ng of input control RNA (IN, average RNA concentration = 41.3 ± 24.6 ng/μL in 16 μL elution volume) per fifty homogenized Tg(myo6b:RiboTag) larvae (n = 19 IPs). Because of the amount of homogenate used for immunoprecipitation, these quantities of RNA correspond to ~330 ng of input RNA and ~1 ng immunoprecipitated RNA per larvae. Interestingly, while analysis of the IN samples showed characteristic ratio of 1.6–2 between the 28s to 18s rRNA, representing intact RNA components of the 60S and 40S ribosomal subunits respectively, the IP samples showed ratios >2.5, indicating reduced levels of 18s rRNA (Figure 2B). This observation is consistent with the zebrafish ribosomal immunoprecipitation results in Tryon et al. (2013), and is thought to be a result of tagging the 60S subunit of the ribosome (of which Rpl10a is a component). Overall, these results indicate that high quality RNA in amounts suitable for downstream analyses such as RT-qPCR and RNA sequencing (RNA-Seq) can be obtained using the Tg(myo6b:RiboTag) fish along with our immunoprecipitation protocol.

Figure 2. HC ribosome immunoprecipitation pulls down HC-specific transcripts. (A) Explanation of the HA-tagged ribosome immunoprecipitation protocol (see Methods for more information). (B) Agilent Bioanalyzer PicoChip outputs showing that input (IN) and immunoprecipitated (IP) RNA quality is comparable between sample types despite reduced levels of 18S rRNA in the IP samples. This observation is consistent with previous zebrafish ribosome immunoprecipitation protocols (Tryon et al., 2013). (C) RT-qPCR results showing that transcripts for known HC-expressed genes such as atoh1a and myo6b are significantly enriched by HA-tagged ribosome immunoprecipitation, whereas transcripts for genes not specifically expressed in HCs are either not significantly enrichened (sox2, myod1) or depleted (rho, vil1). Error bars represent fold change ± standard deviation, and statistical significance was assessed by two-tailed Welch's t-test (n = 10). *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001.

After establishing our protocol, we next aimed to validate whether the immunoprecipitation step was efficiently pulling down RNA from HCs. For this analysis, we used RT-qPCR to determine the relative transcript abundance of known HC expressed genes between IN and IP samples. As expected, immunoprecipitation of HA-tagged ribosomes resulted in significant enrichment of transcripts for the HC expressed genes atoh1a, which encodes a transcription factor necessary for HC development, and myo6b, the promoter of which is used to drive rpl10a-3xHA expression in this model (Figure 2C) (Seiler et al., 2004; Millimaki et al., 2007). We did not observe significant differences in transcript abundance for the supporting cell expressed gene sox2 or the muscle expressed gene myod1. However, we observed significant depletion of the eye specific gene rho and gut specific gene and vil1 in the RNA obtained from the IP compared to the IN. These results indicate that the protocol adapted for the Tg(myo6b:RiboTag) ribosome immunoprecipitation is able to specifically enrich for transcripts of HC expressed genes, while also depleting, to a varying extent, transcripts expressed in other cell types.

RNA Sequencing of IP and in Samples From Tg(myo6b:RiboTag) Zebrafish

After confirming the efficiency of immunoprecipitation protocol, we next sought to perform an unbiased and systematic analysis of the utility of the Tg(myo6b:RiboTag) model as a tool to study the zebrafish HC translatome. We therefore performed RNA-Seq on IP and IN samples in atleast three biological replicates from 5 dpf Tg(myo6b:RiboTag) zebrafish. In total, 17,164 genes were detected as expressed in our IN and IP samples based on our criteria (IN and IP counts per million [CPM] > 0.01 in all samples), allowing for a direct comparison of expression levels (see Supplementary Data Sheet 1). Of these, transcripts for 2,379 genes were significantly enriched in IP samples compared to IN (fold change > 2, IP CPM > 1, false discovery rate [FDR] < 0.05), and transcripts of 2,258 genes were significantly depleted in IP samples compared to IN (fold change < 0.5, IN CPM > 1, FDR < 0.05) (Figure 3A). As an internal control, we first analyzed transcript enrichment and depletion of the same HC and non-HC expressed genes used in our RT-qPCR analysis (Figure 2C) and saw parallel results (Figure 3B). Our RNA-Seq analysis also found that transcripts for the known HC genes atoh1a and myo6b are significantly enriched in the IP samples compared to IN.

Figure 3. HC ribosome immunoprecipitation can reliably detect HC enriched and depleted transcripts by RNA-Seq. (A) Bar graph showing the distribution of significantly enriched and depleted transcripts in the IN vs. IP samples binned by fold change range. For genes with depleted transcripts in the IP samples, the number of genes per bin is as follows: 2x − 5x = 1726, 5x − 10x = 422, 5x − 10x = 109, and 50x + = 1. For IP enriched gene transcripts, the number of genes per bin is as follows: 2x − 5x = 1685, 5x − 10x = 450, 5x − 10x = 223, and 50x + = 21. (B) RNA-Seq fold change enrichment and depletion of HC expressed and non-expressed genes replicates the results obtained by RT-qPCR (Figure 2C). Statistical significance was assessed using DEseq (see Methods). *FDR < 0.05, **FDR < 0.01, ***FDR < 0.001.

As shown in Figure 3A, the majority of transcripts found to be significantly enriched or depleted in the IP compared to the IN samples had a two to five-fold-change in transcript abundance. This large number of genes with small differences in transcript abundance between IP and IN could be attributed to the likelihood that many HC-expressed genes are also expressed in other cell types. Because our protocol isolates RNA from whole, intact larvae, concomitant enrichment and depletion of RNA expressed in both HCs and other cell types could lower the overall fold change values for HC enrichment and depletion. Therefore, in order to restrict our analyses to genes with a preferential expression in HCs compared to other tissues, we focused on genes with high fold-change transcript enrichment in the IP samples. A gene ontology (GO) analysis of the genes with five-fold or greater transcript enrichment in IP compared to IN samples (n = 694) identified “detection of mechanical stimulus” as a top enriched GO term (Table 1, GO associated with HC function in bold). This same analysis was also performed on the genes with five-fold or greater depletion in the IP compared to IN samples (n = 532), resulting in top GO terms including “detection of light stimulus” and “visual perception” (Table 2, in bold). Together with our RT-qPCR, our RNA-Seq, and GO analysis of five-fold enriched transcripts support that our Tg(myo6b:RiboTag) model and immunoprecipitation protocol can enrich for the zebrafish HC translatome.

Table 1. Gene ontology analysis of hair cell enriched genes.

Table 2. Gene ontology analysis of hair cell depleted genes.

Immunoprecipitation Enriches for Known HC Expressed Transcripts

We next wanted to assess the efficiency of our model for enriching for other genes known to be expressed in zebrafish HCs. Therefore, we analyzed our transcripts with five-fold or greater enrichment (n = 694) and depletion (n = 532), from here on name “enriched” and “depleted” gene sets, using the Zebrafish Expression Ontology of Gene Sets (ZEOGS) tool (Prykhozhij et al., 2013). This online tool can be used to detect overrepresented anatomical structures based on known gene expression patterns from the Zebrafish Information Network (ZFIN) (Howe et al., 2013). When analyzing the enriched gene set, the top result for overrepresented anatomy is “neuromast,” and also includes the inner ear HC types “hair cell anterior macula” and “hair cell posterior macula” (Table 3). We also performed a ZEOGS analysis on the highly depleted genes, resulting in the top term of “retinal photoreceptor layer” and no results related to the inner ear or neuromast hair cells (Supplementary Table 2). This indicates that the Tg(myo6b:RiboTag) immunoprecipitation is capturing transcripts of genes experimentally validated to be expressed in neuromast and inner ear HCs.

Table 3. ZEOGS analysis of HC enriched transcripts.

Zebrafish have proved to be a highly useful animal model for studying the normal development and function of HCs, as well as for dissecting the mechanistic consequences of gene mutations that cause human hearing loss (Lieschke and Currie, 2007). For this reason, we next probed our RNA-Seq dataset to detect zebrafish homologs of known human syndromic and non-syndromic deafness causing genes reported on the Hereditary Hearing Loss Homepage (http://hereditaryhearingloss.org) (Van Camp and Smith, 2017). In total, we found that 36 deafness genes were significantly enriched in our zebrafish HC IP samples compared to IN (fold change > 2, p-value < 0.05), of which 27 remained significant after taking into account multiple testing (Table 4, bold values denote FDR < 0.05). Importantly, only 13 of these genes have been previously reported in the literature to be expressed in the developing zebrafish inner ear or HCs, which indicates that the Tg(myo6b:RiboTag) model is effective in identifying genes not previously known to be expressed in zebrafish HCs.

Table 4. Fold change enrichment of human deafness gene homologs in Tg(myo6b:RiboTag) IP samples.

Tg(myo6b:RiboTag) Gene Expression Profile Is Comparable to Previously Published FACS-Based HC Transcriptomes

In order to ensure that similar HC gene enrichment could be obtained from our Tg(myo6b:RiboTag) immunoprecipitation method, we next wanted to compare our translatome dataset to a previously published cell sorting experiment. To do this, we utilized a dataset generated by comparing gene expression in sorted zebrafish HCs to sorted skin cells using microarray technology (Steiner et al., 2014). Setting a five-fold transcript enrichment criterion in HCs compared to skin cells within the microarray dataset, we identified 1,041 unique gene IDs with significant enrichment in sorted HCs (p-value < 0.05 from Steiner et al., 2014 results). We then analyzed the fold change values of this set of genes within our Tg(myo6b:RiboTag) dataset. Of the 1,041 genes identified as significantly enriched in the Steiner et al. dataset, 762 were identified as expressed in our dataset and could be used in this analysis. This is secondary to our strict CPM cutoff to determine gene expression (CPM > 0.01 for all IP and IN samples), as all 1,041 genes are found within our dataset before CPM filtering. Similar to the Steiner et al. dataset, the transcripts encoding these 762 genes also had a higher enrichment in HCs according to our dataset (IP vs. IN, mean Log₂FC = 0.51), compared to all other transcripts detected as expressed (mean Log₂FC = 0.085) (Figure 4A). This result suggests that similar HC expression results can be obtained using either cell sorting or the Tg(myo6b:RiboTag) zebrafish.

Figure 4. Tg(myo6b:RiboTag) immunoprecipitation avoids gene expression changes from cell sorting. (A) Sets of genes found to be enriched in sorted HCs from Steiner et al. (2014) have significantly higher fold change enrichment in the IP samples compared to remaining expressed genes (background). Significance was determined by Wilcoxon's test. (B) Cell sorting experiment utilizing the Tg(myo6b:RiboTag) model. Populations of both GFP positive (GFP+) and GFP negative (Neg, GFP–) were collected via FACS based on GFP expression. Post sort analysis shows high purity of negative and positive populations (98.2 and 95.3%, respectively). (C–E) RT-qPCR analysis of zebrafish homologs of immediate early and heat shock protein encoding genes in sorted HCs (GFP+, C), sorted non-HCs (GFP–, D), and immunoprecipitated RNA (IP, C) vs. RNA extracted from whole, non-dissociated larvae (IN). Error bars represent fold change ± standard deviation. Statistical significance was assessed by two-tailed Welch's t-test (n = 3). *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001.

A possible concern with dissociation and FACS-based approaches for analysis of gene expression relates to induction of changes in gene expression secondary to cellular trauma, loss of tissue context (e.g., cell-cell contacts and lateral inhibition), and length of time from dissociation to RNA extraction. Indeed, a recent manuscript published by van den Brink et al. described robust induction of immediate early gene expression (i.e., Fos, Jun, and Egr1), as well as heat shock protein genes, in mouse muscle stem cells that were dissociated for single cell RNA-Seq (van den Brink et al., 2017). We therefore wanted to compare the expression of the zebrafish homologs of selected immediate early and heat shock protein encoding genes between sorted HCs and immunoprecipitated HC transcripts in the Tg(myo6b:RiboTag) IP. For this analysis, we utilized the concomitant myo6b promoter driven expression of GFP in the HCs of the Tg(myo6b:RiboTag). Five-day-old larvae were dissociated and cells were separated by FACS to GFP-positive (HCs) and GFP-negative (cells from the rest of the fish) populations (Figure 4B). Using RT-qPCR, we found that the immediate early genes egr1, fosab, fosb, and junb, as well as the zebrafish heat shock protein gene hsp70.1 were enriched in sorted HCs compared to their level of expression in intact fish (Figure 4C). Similarly, all five genes were highly expressed in the sorted GFP-negative cells in comparison to their expression in intact zebrafish (Figure 4D), suggesting that their expression is likely induced by the dissociation and cell sorting process. Finally, a similar RT-qPCR analysis performed on the Tg(myo6b:RiboTag) IP vs. IN samples showed that only two of these transcripts, egr1 and fosb are enriched in the IP (Figure 4E). Taken together, these data suggest that (1) immediate early genes are induced in the sorted samples, and (2) that this technical artifact is avoided by measuring gene expression using the Tg(myo6b:RiboTag) zebrafish.

Immunoprecipitation of the HC Translatome Reveals Novel HC Expressed Genes

By comparing our data to genes previously identified as HC-specific in zebrafish, we have demonstrated that the Tg(myo6b:RiboTag) immunoprecipitation method allows for detection of HC expressed genes in zebrafish. However, many of the transcripts that were detected as highly enriched in the IP were not previously reported in the literature as expressed in zebrafish HCs. For this reason, we selected 10 genes that have not been previously characterized as expressed in zebrafish HCs from among our top 100 significantly enriched transcripts (see Supplementary Table 3) for validation of expression by RT-qPCR. These genes were acin1a, cnga1, cnga3a, cnga3b, dynlrb2, grin2db, loxhd1a, loxhd1b, onecut1, and zbtb20. First, we utilized independent Tg(myo6b:RiboTag) IP and IN samples, and found that eight out of the 10 genes validated as significantly enriched in IP compared to IN RNA samples (Figure 5A). To then further confirm the expression of these genes in HCs, we utilized RNA extracted from sorted HCs (post-sort purity = 95.3%, see Figure 4C). All 8 genes that validated as enriched in our IP vs. IN samples were also detected as expressed in sorted zebrafish HCs (cycle threshold [CT] values ≤ 30, Figure 5B). These data further show the utility of the Tg(myo6b:RiboTag) zebrafish in detecting HC expressed genes, as well as for the discovery of new genes with potential novel functions in HCs.

Figure 5. HC RiboTag immunoprecipitation reveals expression of novel HC expressed genes. (A) Validation of 10 genes found to be highly enriched in HCs from the Tg(myo6b:RiboTag) RNA-Seq using RT-qPCR in independent HC immunoprecipitation (IP) and input (IN) samples. Error bars represent fold change ± standard deviation, and statistical significance was assessed by two-tailed Welch's t-test (n = 3). *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001. (B) Validation of the same 10 highly enriched genes selected from the Tg(myo6b:RiboTag) RNA-Seq in sorted HCs by RT-qPCR. Black bar indicates the cycle threshold (CT) cutoff used to denote reliable expression (CT ≤ 30). Error bars represent CT standard deviation between biological replicates (n = 3).

Discussion

Cell type-specific expression analysis has gained popularity in the past decade, and is of particular importance to understanding the biology of rare cell types. The HCs of the zebrafish lateral line and inner ear are one such example of a rare cell type present within a complex tissue environment. The Tg(myo6b:RiboTag) zebrafish described here was generated to express both GFP and HA-tagged ribosomes in inner ear and lateral line HCs (Figures 1A-C). Thus, this model enables HC-specific gene expression analysis through two independent approaches: (1) tissue dissociation and cell sorting based on GFP expression, and (2) immunoprecipitation of HA-tagged ribosomes to enrich for HC-expressed transcripts. While HC-specific transcriptome analysis using flow cytometry is not novel, this is the first application of a RiboTag approach for immunoprecipitation of the zebrafish HC translatome, as well as the first model to allow HC-specific expression analysis using two separate approaches in zebrafish. Through the analyses presented in this manuscript, we have demonstrated the utility of the Tg(myo6b:RiboTag) model in enriching for HC-expressed transcripts (through the immunoprecipitation approach), while concurrently avoiding changes to gene expression that occur secondary to tissue dissociation and the cell sorting process. The immunoprecipitation protocol can be performed using either fresh or frozen tissues, making it a convenient technique for analysis of gene expression at multiple time points or treatment conditions, and removing the need for multiple cell sorting sessions. Our subsequent RNA-Seq and validation experiments revealed not only the expected enrichment of known zebrafish HC expressed transcripts, but also enrichment of transcripts with potential new functions in inner ear and/or lateral line HCs.

In light of these results, some important considerations should also be taken into account when using this model. First, unlike cell sorting where, in ideal conditions, the RNA will contain genetic material only from the cell types of interest, all ribosome immunoprecipitation models, including the Tg(myo6b:RiboTag) zebrafish, are not completely effective at isolating the translatome of a cell type of interest. Immunoprecipitated samples contain, to a varying extent, RNA from other cell types. Indeed, our average IP RNA yield (56 ng per 50 larvae at 5 dpf) was much larger than would be expected if this technique was specifically acquiring only actively translated RNA from HCs.

Second, the inclusion of transcripts from other cell types limits the interpretation of the data. The HC translatome results described in this manuscript were generated by comparing immunoprecipitated HC RNA to the whole fish transcriptome, detecting 2,379 genes as more than two-fold enriched in HCs. This method of comparison is not ideal for detecting all HC expressed transcripts, as it is estimated that an average 11,000–13,000 genes are expressed in all cell types (Ramsköld et al., 2009). Some transcripts which are expressed in HCs may not be detected as enriched if they are expressed in many other cell types included in the sample, as the effect of the depletion may (although not always), outweigh the effect of the enrichment. We therefore recommend that the Tg(myo6b:RiboTag) model be used to analyze gene expression in HCs between different conditions (e.g., drugs or noise) through direct comparison of the IP RNA. The IN RNA can be used to validate whether the changes in gene expression originate from HCs, and we recommend a two-fold enrichment cutoff in IP vs. IN to denote HC expression.

Third use of this model inherently relies on the expression of HA-tagged ribosomes in HCs under control of the myo6b promoter. This promoter is highly specific to HCs to all time points, which is critical to the successful use of this method and the purity of the immunoprecipitated RNA. However, analysis of the HC translatome in different conditions will therefore rely on (1) the presence of HCs and (2) the expression of myo6b. As examples, HC gene expression in very early developmental time points prior to the onset of myo6b, mutants in which HCs never develop or do not express myo6b, or treatments that rapidly decrease HC number or abolish myo6b expression, may be difficult or impossible to analyze with the Tg(myo6b:RiboTag) zebrafish.

Finally, the immunoprecipitation protocol is limited by the total amount of tissue that can be homogenized in one sample (we recommend not exceeding 10% the total volume of homogenization buffer). For analyses of the HC translatome at later developmental stages, different dissection techniques can be applied to reduce total tissue volume and “pre-enrich” for HCs (i.e., inner ear dissection, skin peeling, etc.). These dissections may also be applied to isolate and compare gene expression differences in specific HC subtypes (i.e., inner ear vs. neuromast). Overall, each of these points should be considered when designing experiments using this model.

Our results show that the Tg(myo6b:RiboTag) zebrafish immunoprecipitation protocol detailed in the methods yields quantities of RNA appropriate for downstream analyses of HC gene expression such as RT-qPCR and RNA-Seq. Additionally, this protocol is able to enrich for the zebrafish HC translatome, resulting in marked increases in transcript abundance of known HC expressed genes in the IP samples compared to RNA extracted from whole larvae (IN). Indeed, our RNA-Seq analysis of immunoprecipitated RNA shows enrichment for genes related to inner ear and lateral line HC function, while also showing depletion of genes involved in the function of other organs. Analysis of significantly enriched genes by RT-qPCR revealed high expression of genes not previously validated as expressed in zebrafish HCs, such as loxhd1a, loxhd1b, and dynlrb2, among others. Loxhd1a and loxhd1b are both homologs of the human deafness gene LOXHD1, the causative gene of DFNB77 (Grillet et al., 2009). Additionally, dynlrb2, a gene that encodes for a dynein light chain protein, has been previously identified as highly expressed in mouse HCs at both the transcript and protein levels, although its function in HCs has yet to be elucidated (Jiang et al., 2001; Scheffer et al., 2015; Shen et al., 2015; Hickox et al., 2017). Zebrafish may therefore be an appropriate model to study the roles that these genes play in normal HC function. Overall, taking into account the considerations outlined above, we believe that the Tg(myo6b:RiboTag) zebrafish represent an easy to use, versatile and sensitive model for studying inner ear and lateral line HC gene expression.

Author Contributions

AB and KK generated the model, MM, YO, KK, and RH designed and interpreted the experiments, MM, AB, YO, and NP performed the experiments, YS analyzed the gene expression data, and MM, KK, and RH wrote the manuscript.

Funding

This work was supported by NIDCD/NIH R01DC013817 and R01DC03544 (RH), NIDCD/NIH Training Grant DC00046 (MM), NIDCD/NIH F31DC016218 (MM), and NIDCD/NIH Intramural Research Program Grant 1ZIADC000085-01 (KK).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Dr. David Raible from the University of Washington, as well as Dr. Doris Wu from the NIDCD, for their helpful comments on the manuscript. They would also like to thank Dr. Amiel A. Dror for contributing to the design of Figure 2A, Dr. Brant Weinstein and Dr. Mayumi Miller at NICHD for providing the RiboTag clone, as well as Dr. Xiaoxuan Fan, director of the University of Maryland School of Medicine Flow Cytometry Facility (SOM/UMGCC), for his help and expertise with the FACS experiments.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2018.00047/full#supplementary-material

References