Multidimensional quantitative analysis of mRNA expression within intact vertebrate embryos

Traditional in situ hybridization approaches based on catalytic reporter deposition (CARD) yield high-contrast images of mRNA expression domains within whole-mount vertebrate embryos (Tautz and Pfeifle, 1989; Harland, 1991; Lehmann and Tautz, 1994; Kerstens et al., 1995; Nieto et al., 1996; Pernthaler et al., 2002; Denkers et al., 2004; Kosman et al., 2004; Thisse et al., 2004; Clay and Ramakrishnan, 2005; Barroso-Chinea et al., 2007; Acloque et al., 2008; Piette et al., 2008; Thisse and Thisse, 2008; Weiszmann et al., 2009; Ruf-Zamojski et al., 2015). However, the intensity of the staining is qualitative rather than quantitative; furthermore, spatial resolution is often compromised by diffusion of reporter molecules prior to deposition (Tautz and Pfeifle, 1989; Thisse et al., 2004; Acloque et al., 2008; Piette et al., 2008; Thisse and Thisse, 2008; Weiszmann et al., 2009), and multiplexing is cumbersome, requiring serial staining of each target mRNA (Lehmann and Tautz, 1994; Nieto et al., 1996; Denkers et al., 2004; Kosman et al., 2004; Thisse et al., 2004; Clay and Ramakrishnan, 2005; Barroso-Chinea et al., 2007; Acloque et al., 2008; Piette et al., 2008). These strengths and weaknesses all derive from the enzyme-mediated deposition process responsible for signal amplification. Direct-labeled probes offer complementary trade-offs, avoiding signal amplification to enable quantitative, high-resolution, multiplexed studies in thin samples (Kislauskis et al., 1993; Femino et al., 1998; Levsky et al., 2002; Kosman et al., 2004; Capodieci et al., 2005; Chan et al., 2005; Raj et al., 2008), but often generating insufficient signal to achieve the needed contrast in thick samples such as whole-mount vertebrate embryos.

To quantify relative mRNA expression levels for defined anatomical regions within vertebrate embryos, it is necessary to destroy the sample morphology. Current approaches employ some combination of microdissection (Nawshad et al., 2004; Redmond et al., 2014; Treutlein et al., 2014), cell dissociation (Manoli and Driever, 2012; Jean et al., 2015; Petropoulos et al., 2016), homogenization (Axelsson et al., 2007; de Jong et al., 2010; Pena et al., 2014), fluorescence-activated cell sorting (Manoli and Driever, 2012; Treutlein et al., 2014; Allison et al., 2016), magnetic-activated cell sorting (Treutlein et al., 2014; Allison et al., 2016; Taylor et al., 2016) or lysis (Nawshad et al., 2004; de Jong et al., 2010; Laranjeiro and Whitmore, 2014; Redmond et al., 2014; Treutlein et al., 2014; Jean et al., 2015; Allison et al., 2016; Petropoulos et al., 2016), followed by RNA quantitation using quantitative real-time polymerase chain reaction (qPCR) (Nawshad et al., 2004; Axelsson et al., 2007; Laranjeiro and Whitmore, 2014; Pena et al., 2014; Jean et al., 2015), RNA sequencing (Treutlein et al., 2014; Allison et al., 2016; Petropoulos et al., 2016), in situ hybridization flow cytometry (Allison et al., 2016; Taylor et al., 2016), microarray hybridization (de Jong et al., 2010; Redmond et al., 2014; Jean et al., 2015) or hybridization barcoding (Laranjeiro and Whitmore, 2014; Pena et al., 2014). Owing to this fundamental trade-off between anatomical context and quantitation, there is an unmet need for multiplexed quantitative analysis of mRNA expression with high-resolution within intact specimens.

We have shown previously that in situ hybridization chain reaction (HCR; Fig. 1A) (Dirks and Pierce, 2004; Choi et al., 2010) enables straightforward multiplexing, high contrast and subcellular resolution when mapping target mRNAs within complex specimens (Choi et al., 2014, 2016). In situ HCR uses DNA probes complementary to mRNA targets to trigger the self-assembly of fluorophore-labeled DNA HCR hairpins into tethered fluorescent amplification polymers. Using a library of orthogonal HCR amplifiers, signal amplification is performed for all targets simultaneously. Here, we demonstrate the crucial property that the amplified HCR signal is proportional to the number of target mRNAs per subcellular imaging voxel (Fig. 1B), enabling accurate and precise relative quantitation within intact vertebrate embryos.

We demonstrate in situ mRNA quantitation by examining somitogenesis in the zebrafish embryo. In all vertebrates, somites provide one of the first outward appearances of the metameric body plan, periodically pinching off from the presomitic mesoderm (PSM) in bilaterally symmetrical pairs as precursors to the axial muscles and vertebral column (Oates et al., 2012). To date, detailed studies of the gene dynamics underlying somitogenesis have relied heavily on examination of one-channel (Henry et al., 2002; Oates and Ho, 2002; Mara et al., 2007; Gomez et al., 2008; Ferjentsik et al., 2009; Choorapoikayil et al., 2012; Schroter et al., 2012) and two-channel (Holley et al., 2000; Oates and Ho, 2002; Jülich et al., 2005) CARD images that display expression of one or two mRNAs per embryo. Here, using in situ HCR, we demonstrate quantitative subcellular analyses of four mRNAs in the same embryo. Read-out analyses reveal quantitative changes in gene co-expression ratios as somites mature and move away from the PSM, and read-in analyses reveal the anatomical locations where distinct gene co-expression relationships occur.

Subcellular mRNA quantitation within thick samples such as whole-mount vertebrate embryos has not been rigorously verified for any method. We chose to meet this challenge by exploiting the ease of multiplexing using in situ HCR: we redundantly detect the same target mRNA using two distinct probe sets, each of which carries initiators for orthogonal HCR amplifiers labeled with spectrally distinct fluorophores (Fig. 2A). This experimental design provides an avenue for validating that HCR signal scales linearly with the number of target mRNAs per voxel, without requiring knowledge of the absolute number of targets in any voxel.

Using this approach, we perform two-channel redundant detection of four different target mRNAs within whole-mount zebrafish or mouse embryos (Fig. 2B), observing highly correlated subcellular voxel intensities (Fig. 2C; Pearson correlation coefficient 0.91≤r≤0.97 for 2×2×2 µm voxels). Control experiments suggest the absence of significant penetration or crowding effects, supporting the conclusion that the HCR signal is quantitative (see section S2.1 in the supplementary material). For each target mRNA in Fig. 2C, variation along the diagonal indicates biological variation in expression levels between voxels; the accuracy is high (a clear linear relationship with intercept near zero) and the precision is very good (scatter of ∼5% to 20% of the dynamic range depending on the target mRNA). The scatter in each channel arises from non-uniformity in the background per voxel, as well as from non-uniformity in the signal generated per target molecule [resulting from variation in probe hybridization yields (Raj et al., 2008; Shah et al., 2016) and HCR amplification polymer lengths (Dirks and Pierce, 2004; Choi et al., 2010, 2014)]. Scatter is reduced by averaging neighboring pixel intensities, enabling dramatic improvements in precision while still retaining subcellular resolution (Fig. 2D). This averaging effect is also evident in the high precision achieved using HCR to perform relative and absolute quantitation of miRNAs in northern blots (Schwarzkopf and Pierce, 2016), where the size of the voxel is effectively the size of the band being quantified. We expand upon the origins and properties of this important averaging effect in section S2.2 of the supplementary material.

We next tested our ability to use in situ HCR to discriminate a known difference in mRNA levels between embryos. The FlipTrap line Gt(desma-citrine) (Trinh et al., 2011) provides an ideal test setting as the expression level for citrine is expected to be approximately twofold higher in homozygous versus heterozygous embryos (expressed from two or one transgenic alleles), whereas desma expression is expected to be similar in both genotypes (expressed from both alleles regardless). Extending our two-channel redundant detection approach, we detect citrine with one probe set and desma with a second probe set (Fig. 3A). Comparing mean citrine expression levels within regions of high expression (Fig. 3B) yields a homo/hetero expression ratio of 2.0±0.5 (mean±s.d., n=3 embryos; Fig. 3C). This roughly twofold difference in expression level is also evident in the distinct slopes observed in scatter plots of subcellular voxel intensities (Fig. 3D). The ability to discriminate twofold changes in mRNA expression in an anatomical context is important, for example, in evaluating perturbations or candidate drugs intended to alter gene expression.

We next apply in situ HCR to multidimensional analyses of zebrafish somitogenesis to explore the power of mRNA quantitation in an anatomical context. Fig. 4A displays a four-channel image for four target mRNAs expressed as the somites emerge and are displaced from the PSM: two cyclic segmentation genes (her7 and her1) and two muscle genes (myod1 and tpm3). Fig. 4B presents expression profiles in a strip spanning five regions of interest (somites S7, S8, S9 and S10, and the PSM), revealing the relative expression levels from the older somite (S7) to the tissue soon to become a somite: ∼100%, 75%, 50%, 25%, 0% for myod1 (magenta curve) and ∼0%, 0%, 0%, 50%, 100% for her7 (red curve). These profiles provide the first quantitative glimpse of the strongly anti-correlated bulk expression trends for myod1 and her7, but do not take advantage of the subcellular resolution of the data.

To perform subcellular quantitative analyses, we display scatter plots of voxel intensities reminiscent of multiplexed fluorescence-activated cell sorting (FACS) data. Fig. 4C demonstrates read-out from a region of interest in the four-channel image (in this case, a rectangle containing somites S9 and S10, and a longitudinal stripe of adaxial cells) to scatter plots of normalized voxel intensities for pairs of target mRNAs. Strikingly, these expression scatter plots reveal well-defined expression clusters with differing slopes and amplitudes corresponding to voxels with related expression characteristics. For example, the her7-myod1 quadrant reveals one cluster with low myod1 and variable her7 levels, a second cluster with low her7 and variable myod1 levels, and a third cluster with correlated variable expression of both targets. When examining quantitative expression scatter plots generated by read-out, a researcher is free to select any subset of voxels as an expression cluster of interest for subsequent read-in.

Unlike FACS analysis, where the display of expression clusters would be the final product, in situ HCR maintains the anatomical context, permitting us to map expression clusters of interest back into the embryo to interactively investigate their physical positions. Fig. 4D demonstrates read-in from four selected expression clusters in the her7-myod1 quadrant (shaded lilac, yellow, orange and cyan) back into the embryo, segmenting the image into four expression domains. The lilac expression cluster maps to the posterior region of S10, and the yellow expression cluster maps to S9 and the anterior region of S10, placing the boundary between these expression domains in the middle of the youngest somite. The orange and cyan expression clusters map to the adjacent regions of the longitudinal stripe of adaxial cells. Cluster shading can be propagated from the her7-myod1 quadrant to the other three quadrants, which is akin to projecting the voxel intensities onto four axes, but without the difficulty of four-dimensional visualization. For example, the well-separated lilac and yellow clusters identified in the her7-myod1 quadrant partition the complex cluster near the origin in the myod1-tpm3 quadrant into two clusters with differing slopes (Fig. 4D). Similar read-in analyses can be performed for expression clusters identified in any quadrant or combination of quadrants, permitting a detailed examination of the spatial organization of distinct genetic circuit states (Figs S34-S40).

Somites are rhythmically and sequentially generated by the presomitic mesoderm (PSM), resulting in a developmental time course being reflected within a single embryo (somite S7 is developmentally more mature than S8, which is more mature than S9, and so on), with new somites emerging at approximately 30-min intervals in zebrafish (Oates et al., 2012). Having previously used expression clusters within the scatter plots to identify expression regions within the image (Fig. 4D), we now reverse the direction of information flow and use each somite within the image to identify expression clusters within the scatter plots. This approach yields quantitative snapshots of gene co-expression changes as the somites form and mature. Fig. 5A depicts four-channel read-out from five regions (S7, S8, S9, S10 and PSM) to expression clusters shaded by their anatomical location (black, cyan, green, yellow and orange). The myod1-tpm3 quadrant reveals striking changes in slope and amplitude during this maturation process (Fig. 5B): a low slope in the PSM jumps to a higher slope in maturing somites S9, S8 and S7; within these three older somites, amplitude increases monotonically with somite maturity. The youngest somite (S10) exhibits expression clusters with both slopes, consistent with the observations of Fig. 4D. The similar slopes and small intercept for the three older somites indicate that the ratio of myod1 to tpm3 expression remains approximately constant as the somites mature. Interestingly, replicate embryos of nominally the same age capture slightly different developmental stages within the oscillatory somitogenesis circuitry, revealing expression clusters with slopes and amplitudes that reflect related but different circuit states (Fig. S42). Consistent with the expectation that somitogenesis is bilaterally symmetrical (Oates et al., 2012), similar expression clusters are observed for left and right somites within an embryo (Figs S42-S46); it will be interesting to see whether more detailed quantitative analyses reveal asymmetries.

Discerning two subcircuits emerging in the expression clusters of Fig. 5B, we project the voxel intensities onto new axes to examine somitogenesis through the lens of a her1-her7 subcircuit and a myod1-tpm3 subcircuit (Fig. 5C). There is a crossover within the youngest somite (S10) between dominance of the her1-her7 subcircuit in the PSM and the myod1-tpm3 subcircuit in maturing somites S9, S8 and S7, leading to nearly orthogonal expression clusters in Fig. 5C. Subcircuit projections can be performed onto axes representing an arbitrary number of circuit elements, compactly summarizing a large quantity of high-dimensional expression information.

qHCR imaging dramatically expands the capabilities of in situ hybridization as a research tool. mRNA expression levels can be quantitatively compared within an anatomical region (Fig. 2), between regions (Fig. 4B) and between embryos (Fig. 3). Quantitative discovery is enabled in two directions: read-out from multi-channel images to co-expression scatter plots reveals well-defined expression clusters with differing slopes and amplitudes (Fig. 4C and Fig. 5); read-in from expression scatter plots to multi-channel images reveals the anatomical locations in which selected co-expression relationships occur (Fig. 4D). In one direction, anatomical locations within the image can be used to identify expression clusters within scatter plots (Fig. 5), and in the other direction, expression clusters within the scatter plots can be used to identify anatomical locations within the image (Fig. 4D). These capabilities result from two crucial properties of in situ HCR: multiplexing spreads the voxel data out onto multiple expression axes to enable identification of well-segregated expression clusters; quantitation generates expression clusters that display slopes and amplitudes revealing similarities and differences between genetic circuit states. Projection of expression scatter plots onto subcircuit axes representing one or more target mRNAs facilitates quantitative dissection of the underlying regulatory circuitry across multiple anatomical regions. Collectively, these capabilities open a new era for in situ hybridization, enabling quantitative bi-directional interrogation of anatomical locations and expression clusters in the study of vertebrate development and disease.

For each target mRNA, a kit containing a DNA probe set, a DNA HCR amplifier (comprising a pair of fluorophore-labeled DNA hairpins), and hybridization, wash and amplification buffers was purchased from Molecular Instruments (molecularinstruments.org). See Table S1 for a summary of sample, probe set and amplifier details; see section S3 in the supplementary material for probe sequences. Sequences for HCR amplifiers B1, B2, B3 and B4 are given in Choi et al. (2014). The size of each probe set was based on the expression level of the target, the affinity and selectivity of the probes, and the level of autofluorescence in the channel corresponding to the target.

Procedures for the care and use of zebrafish embryos were approved by the Caltech IACUC and the USC IACUC. Wild-type Danio rerio embryos were either strain TL [for crowding studies (see section S2.1.2 in the supplementary material)] or AB strain (all other wild-type studies); transgenic embryos were obtained from a previously described screen (Trinh et al., 2011). With the exception of sample mounting, in situ HCR was performed using the protocol for whole-mount zebrafish embryos provided by Choi et al. (2016). To mount embryos fixed at 26 hpf, 1 ml of 1% agarose (Invitrogen, 16500-100) solution (w/v) was placed in an imaging dish with a #1.5 coverglass bottom (WillCo Wells, GWST-5040, 0.01 mm) and then a custom negative plastic mold (Megason, 2009) was placed on top to create grooves in the agarose. Following solidification of the agarose at room temperature (5 min) and removal of the mold, a hair loop was used to position 10-12 embryos in grooves for imaging. For embryos fixed at 10 hpf, embryos were de-yolked (to avoid high autofluorescence from the yolk) and then flattened onto a #1.5 coverslip (Sigma-Aldrich, Z692247) using a hair loop. To immobilize the samples, one or two drops of 1% agarose were added on top of the embryos and allowed to solidify.

Procedures for the care and use of mouse embryos were approved by the USC IACUC. Mouse embryos were Mus musculus strain 129S4.Cg-Tg(Wnt1-cre)2Sor/J. In situ HCR was performed using the protocols for whole-mount mouse embryos (E9.5) provided by Choi et al. (2016).

Images were acquired within whole-mount zebrafish embryos using a Zeiss 710 NLO inverted confocal microscope with an LD LCI Plan-Apochromat 25×/0.8 Imm Korr DIC objective, with the exception that for Figs S3-S7, images were acquired using a Zeiss 700 inverted confocal microscope with a Plan-Apochromat 20×/0.8 M27 objective. Images were acquired within whole-mount mouse embryos using a Zeiss 710 NLO inverted confocal microscope with an LD C-Apochromat 40×/1.1 W Korr M27 objective. See Table S2 for a summary of the confocal microscope settings used for each target mRNA. Images are displayed without subtracting background. Images for Figs 4, 5 and Figs S29-S47 are displayed as mean intensity projections over five focal planes to facilitate comparison of multiple targets expressed at closely related depths within the somites.

Image analysis was performed as detailed in section S1.3 of the supplementary material and in Tables S3 and S4, including calculation of raw voxel intensities, characterization of signal and background, calculation of normalized voxel intensities, signal-to-background analysis, and signal-to-signal analysis.

See Fig. 6 for an illustration of the read-out/read-in workflow.

Acquire and normalize data. Step 0a: Perform a multiplexed in situ HCR experiment. Image the expression patterns for n target mRNAs using an n-channel experiment in a single specimen (Choi et al., 2016) (three channels depicted in Fig. 6). Step 0b: Calculate normalized voxel intensities for each channel. Average pixel intensities to create raw voxel intensities for each channel (see section S1.3.1 in the supplementary material) and then calculate the normalized signal estimate for each voxel in each channel (see section S1.3.3 in the supplementary material). The same normalization should be applied to all images to enable quantitative comparison between embryos.

Read-out from anatomical space to expression space. Step 1a: Select anatomical regions of interest within a multi-channel image (two regions of interest depicted in Fig. 6: tan and gray rectangles). Step 1b: Plot expression scatter plot for each pair of channels: for each voxel in the regions of interest, plot the normalized signal intensity shaded by region color (tan or gray dots in Fig. 6; see also examples in Fig. 4C and Fig. 5B).

Read-in from expression space to anatomical space. Step 2a: Select expression clusters of interest (magenta and orange clusters selected in Ch3-Ch1 scatter plot, dots not selected are shaded blue). Step 2b: Redisplay the multichannel image showing only the voxels in one expression cluster at a time (images bounded by magenta, orange or blue rectangles; see also example in Fig. 4D).

Read-out/read-in analyses and associated subcircuit analyses can be performed using the Matlab script (readoutreadin.m), which is available for download at molecularinstruments.org, including a user guide that provides step-by-step instructions for running the script, screen shots that illustrate these steps, a sample image for testing the script and sample output files.

This work is dedicated to our friend and colleague, Robert M. Dirks, who tragically lost his life in the Metro-North train crash in New York in 2015. We thank L. A. Trinh and B. R. Wolfe for helpful discussions, M. Schwarzkopf for experimental assistance, and L. A. Fletcher, C. Paquette and D. Mayorga for fish care. The findings are those of the authors and should not be interpreted as representing the official views or policies of the U.S. Government.