Zebrafish were propagated and maintained as described (Westerfield, 2007), on recirculating, monitored, and filtered system water, on a 14:10 light/dark cycle, at 28.5°C. Procedures involving animals were approved by the Animal Care and Use Committees of the University of Idaho and of the University of Texas Health Science Center at Houston. Wild-type (WT) zebrafish were of a strain originally provided by Scientific Hatcheries (now Aquatica Tropicals, Plant City, FL) and AB (RRID:ZIRC_ZL1) from the Zebrafish International Resource Center at the University of Oregon. The lws:PAC(H) transgenic line harbors a PAC clone that encompasses the lws locus, modified such that a GFP-polyA sequence, inserted after the lws1 promoter, reports expression of lws1, and an RFP (dsRedExpress)-polyA sequence, inserted after the lws2 promoter, reports expression of lws2 (Tsujimura et al., 2010). This line was the kind gift of Shoji Kawamura and the RIKEN international resource facility. The thrb2:tdTomato transgenic line expresses the tdTomato reporter under control of the thyroid hormone receptor beta 2 promoter, resulting in tdTomato in all adult LWS cones (Suzuki et al., 2013). This line was the kind gift of Rachel Wong. In this study larval (4 dpf) and adult (0.5–1.5 years; both sexes) zebrafish were used.

Dissociation and FACS were carried out as previously reported (Sun et al., 2018a,b), for bulk RNA-Seq and for qPCR. In brief, adult lws:PAC(H) and thrb2:tdTomato zebrafish were collected near the time of light onset but maintained in the dark (dark-adapted), euthanized with MS-222, and retinal tissues dissected away from other ocular tissues including the RPE and collected into cold RNAse-free phosphate-buffered saline (PBS). Retinas were dissociated for 10 min at 37°C in a filtered buffer containing papain, trypsin, neutral protease, catalase, and superoxide dismutase. The reaction was quenched with heat-inactivated fetal bovine serum, and samples resuspended in DNAseI for 10 min at room temperature. Samples were pelleted and resuspended in RNAse-free PBS prior to FACS. Lws:PAC(H) samples were sorted with a SONY Cell Sorter SH800 based upon GFP and RFP fluorescence, and collected into TRIzol LS or lysis buffer from the Machery–Nagel RNA extraction kit (Sun et al., 2018b). Thrb2:tdTomato samples were sorted using the same instrument and conditions, but based upon tdTomato fluorescence intensity and the scatter characteristics (Sun et al., 2018b).

The quality of isolated RNA was evaluated using the Bioanalyzer 2100 RNA 6000 Nano assay (Agilent Technologies). For samples with a RIN score greater than 8.0, sequencing libraries were constructed according to the manufacturer’s protocol using the TruSeq^® RNA Library Preparation Kit (Illumina) with 5 μg total RNA as input whereby ribosomal RNA was removed by poly-A selection using Illumina. Illumina sequencing adapters were ligated to each sample. Ligated fragments were then amplified for 12 cycles using primers incorporating unique dual index tags. Libraries were sequenced at the National Eye Institute (NEI) on the Illumina HiSeq 2500 platform, and raw sequencing reads were demultiplexed by NEI. Raw bulk RNA-Seq reads were trimmed for Illumina adapters and quality using Trimmomatic v0.36 (Bolger et al., 2014). Trimmed reads were aligned to the zebrafish reference genome GRCz11 using STAR v2.5.2a (Dobin et al., 2013) and Salmon v1.0.0 (Patro et al., 2017). Approximately 1–3 million reads were mapped per sample library (3 libraries per condition) to zebrafish reference genome GRCz11, with the exception of one lws1:GFP + library, which still provided minimal depth for downstream analyses, and so was included in subsequent analyses to increase statistical power. Additional quality control metrics were evaluated using FastQC v0.11.5¹ and HTStream (Petersen et al., 2015). Data were imported into R (R Core Team, 2017)² using tximport (Soneson et al., 2015). Differential expression (DE) analysis was then carried out using DESeq2 (Love et al., 2014). Log₂FC, as well as a moderated Log₂FC (to normalize for transcripts displaying very low levels of expression) were determined for lws1:GFP vs. lws2:RFP DE analyses. As additional strategies for identifying transcripts DE in LWS1 vs. LWS2 cones, DE analyses were also carried out for the entire LWS cone population (thrb2:tdTomato +) vs. lws1:GFP and vs. lws2:GFP. Gene ontology analyses were performed using gProfiler.

Single-cell library preparation and data analysis are fully described in Santhanam et al. (2022). In brief, two female and one male wild-type AB strain zebrafish, age 7 months, were euthanized with MS222 on ice and isolated retina-RPE preparations were pooled in a 3:1 mix of Leibovitz’s L-15 medium (Gibco) and Earle’s Balanced Salt Solution (EBSS; Gibco). Cells were dissociated with papain solution (50 U/mL; Worthington Biochemical) for 60 min at 28°C with gentle trituration. Digestion was halted by addition of 2x volume of 0.1% BSA in L15: EBSS medium, and cells were counted and tested for viability. Cell samples were submitted for 3′ scRNA library preparation and sequencing through the Baylor College of Medicine Single Cell Genomics Core facility. The single-cell library was prepared using a Chromium Next GEM Single Cell 3′ Reagent Kit v2 (10x Genomics, Pleasanton, CA, USA). The single-cell library was sequenced with Illumina HiSeq2500.

The single-cell library sequences were initially analyzed using the 10x Genomics CellRanger V2.1.1.0 pipeline. Sequences were aligned to the zebrafish reference genome GRCz11 using CellRanger count, and quality-checked using FASTQC V0.11.9. Over 93% of reads mapped successfully to the genome. The initial alignment and analysis were performed through the Texas Advanced Computing Center (TACC) Lonestar5 computing service. Subsequent analyses of aligned data used the Seurat V2.1.1 (Butler et al., 2018) package in R V3.6.1. The data were initially filtered to remove low-abundance genes (expressed in fewer than 10 cells), doublets and cells with >5% mitochondrial genes. The dataset contained between 200 and 4,000 genes per cell, and 13,551 cells were sequenced/analyzed. PCElbowPlots were performed and 20 principal components were used for downstream analysis of each dataset. PC1 to PC20 were used to construct nearest neighbor graphs in the PCA space followed by Louvain clustering and non-linear dimensional reduction by TSNE to visualize and explore the clusters. Expression levels are expressed in a base 2 log scale.

Stock solutions of tri-iodothyronine (T3) were prepared in DMSO (Sigma), and maintained at −20°C in the dark. Embryos were obtained from WT crosses, with the time of light onset considered the time of fertilization. 0.003% phenylthiourea (PTU) was added to system water at 10–12 h post-fertilization (hpf) to inhibit melanin synthesis [13]. Prior to T3 treatment, embryos were dechorionated using fine forceps, and the 1000X T3 stock solution was added to system water for a final concentration of 100 nM (DMSO final concentration was 0.1%). Controls were treated with 0.1% DMSO. Treatments took place from 48 to 96 hpf, and solutions were refreshed every 24 h (Mackin et al., 2019).

Total RNA from larval (4 dpf) zebrafish tissues was extracted using the Machery-Nagel kit, and then the Superscript III/IV (Invitrogen) was used to synthesize cDNA template with random primers. Gene-specific primer pairs for qPCR are provided in Supplementary Table 1. Amplification was performed on a StepOne Real-Time PCR system using SYBR Green or Power Track SYBR Green master mix (Applied Biosystems). Quantification of transcript abundance was relative to the reference transcript (β-actin), using the ddCT method. Graphing and statistics were performed in Excel. Sample groups were evaluated for normal distributions using the Shapiro–Wilk test. For comparisons showing normal distributions, p-values were calculated using Student’s t-test, and for comparisons not showing normal distributions, p-values were calculated using Mann–Whitney tests.

Hybridization chain reaction procedures were carried out according to the manufacturer’s instructions (Molecular Instruments) (Choi et al., 2014). In brief, zebrafish tissues were fixed overnight in phosphate-buffered 4% paraformaldehyde at 4°C. Tissues were then washed in PBS, dehydrated in MeOH, and stored in MeOH at −20°C at least overnight. Tissues were rehydrated in a graded MeOH/PBS/0.1% Tween 20 series, permeabilized with proteinase K (larvae only), and post-fixed with 4% paraformaldehyde prior to hybridization. Hybridization was done overnight at 37°C. Tissues were washed with the manufacturer’s wash buffer, and then 5XSSCT (standard sodium citrate with 0.1% Tween-20), and the amplification/chain reaction steps were performed following the manufacturer’s protocol. Probe sets were designed and generated by Molecular Instruments and can be ordered directly from their website.

Whole, HCR-processed, adult (0.5–1.5 years) retinas were mounted in glycerol and imaged with a 20X dry lens, 40X water-immersion lens, or 40x oil-immersion lens. 3 μm-step sizes were used for 20X images. Z-series for 40x images were taken with ≤1 μm step size. Whole larval eyes were removed from HCR-processed embryos, and the sclera removed by microdissection. Eyes were mounted in glycerol and imaged with a 20X dry lens using a Nikon-Andor spinning disk confocal microscope and Zyla sCMOS camera. A z-series encompassing the entire globe of the embryo eye was obtained with 3-μm step sizes, using Nikon Elements software. Z-stacks were flattened by max projection, and brightness/contrast adjusted in FIJI (ImageJ). Cross sectional images for adult whole mounted retinas were obtained using the orthogonal view function in FIJI.

Presence/absence of gngt2a expression in dorsal regions of larval eyes was determined using whole eye stacks. Brightness/contrast were adjusted to determine presence/absence of signal when signal was dim. Presence within dorsal retina was defined as signal localized in the center of the region dorsal to the center of the lens. A proportion test was used to determine statistical significance.

Hybridization chain reaction fluorescence intensity quantification was performed using FIJI, using the approach of Thiel et al. (2022). For each image the color channels were split, and background measurements were performed on the channel reporting expression of the gene of interest at 3 positions in the z dimension. Whole eyes were traced in the segmentation editor to create an “object” for measurement. The 3D intensity measure plugin was used to obtain “intensity sum,” also known as integrated density, as well as total volume. The single channel reporting expression of the gene of interest was used as the “signal” for 3D analysis. The mean fluorescence intensity of the background was multiplied by the volume of the eye to determine total background. Corrected total fluorescence was calculated by subtracting total background from measured intensity sum (Thiel et al., 2022).

Graphing and statistics were performed in Excel. Sample groups were evaluated for normal distributions using the Shapiro–Wilk test. For comparisons showing normal distributions, p-values were calculated using two-tailed Student’s t-test, and for comparisons not showing normal distributions, p-values were calculated using Mann–Whitney tests. ^*** denotes p < 0.001, ^** denotes p < 0.01, * denotes p < 0.05.

To address the hypothesis that LWS1 and LWS2 cone subtypes of the zebrafish are distinct in transcriptional characteristics other than opsin expression, we aimed to identify any genes that are DE in these two cone populations of adult zebrafish. LWS1 vs. LWS2 cones of the lws:PAC(H) transgenic line were FACS-sorted based upon GFP vs. RFP fluorescence (Sun et al., 2018b), and RNA isolated from the sorted cones was sequenced to discover their respective transcriptomes (Figures 1A, B). The results of the dissociation, sorting, evaluation of purity, and RNA quality were reported earlier (Sun et al., 2018b). Based upon a false discovery rate (FDR) < 0.05, and absolute fold-change (FC) > 2, 95 transcripts were enriched in GFP + (LWS1) cones, representing ∼0.6% of the LWS1 transcriptome (Figure 1C, Table 1, and Dataset 1. Note that Tables 1, 2 and Datasets 1–3 show log₂FC rather than absolute FC). Using the same cutoff criteria, 186 transcripts were enriched in RFP + (LWS2) cones, representing ∼1.2% of the LWS2 transcriptome (Figure 1C, Table 2, and Dataset 1). These analyses suggest that these cone subtypes are indeed highly similar, although with sufficient transcriptional differences other than opsin expression, to support our original hypothesis. A moderated Log₂FC approach, to normalize for transcripts expressed at very low levels, returned fewer, but several of the same transcripts (Dataset 2; 19 genes in common with Dataset 1 as LWS1-enriched; 56 genes in common with Dataset 1 as LWS2-enriched). The entire dataset is publicly accessible (GEO accession #GSE232902).³

Transcripts enriched in LWS1 cones included those with functions in phototransduction [gngt2a, encoding one of the γ subunits of transducin; (Lagman et al., 2015)], circadian rhythms (per1b, cry1bb), cell adhesion (plxna1a, ephrin-A1a), and transcriptional regulation (foxg1a, hmgb1b, rorcb) (Table 1 and Datasets 1, 2). Notable LWS1-enriched transcripts included two with functions related to TH signaling (nrip1a, thrab), which is a powerful regulator of lws1 vs. lws2 expression (Mackin et al., 2019). A gene ontology (GO) analysis returned one overrepresented molecular function category (transcription factor binding), and several cellular process categories related to differentiation and neurogenesis (Figure 1D). Transcripts enriched in LWS2 cones also included those with functions in phototransduction (gngt2b, cngb3), cell adhesion (adgrl3.1, nptna), and transcriptional regulation related to nuclear hormone signaling (nr2f2, nr4a3), although none with functions in circadian rhythms (Table 2 and Datasets 1, 2). GO analysis of LWS2-enriched transcripts will be discussed below. Therefore, in addition to the divergent λ_max of LWS1 vs. LWS2 cones (Chinen et al., 2003), the two populations may have further distinctions in phototransduction kinetics, cell-cell contacts, and transcriptional regulation by nuclear hormone receptors.

Selected transcripts were analyzed from two additional sorting experiments by qPCR. The results from Sort #1 were previously reported, and validated the presence of opn1lw1, but absence of opn1lw2 (along with two rh2-type opsin transcripts) in the GFP + (LWS1) cones, and the presence of opn1lw2, but absence of opn1lw1 (along with two rh2-type opsin transcripts) in LWS2 cones (Sun et al., 2018b). Sort #2 verified depletion of gngt2a, arr3a, and neurod1, and enrichment of opn1lw2 and gngt2b in LWS2 (RFP +) cones (Supplementary Figure 1). Curiously, opn1lw1 was not detected as DE by either the RNA-Seq analysis, or by qPCR of Sort #2 (Supplementary Figure 1). Although we only rarely observe co-expression of GFP and RFP reporters in adult lws:PAC(H) retinas (Tsujimura et al., 2010; Mitchell et al., 2015; Stenkamp et al., 2021), we tested whether co-expression was more common for the native transcripts, using multiplex fluorescence in situ hybridization. These studies confirmed that a small fraction of the adult LWS cones indeed express both lws transcripts. Further, we note that the fish used for the RNA-Seq studies were sacrificed in the morning, a time of reduced levels of cone opsin transcript expression (Li et al., 2008), perhaps affecting the likelihood of detecting lws1 as DE.

Another feature of the DE list that represents a possible limitation of the FACS-RNA-Seq approach is the abundance of DE transcripts encoding components of the proteasome (37 of the 185 transcripts enriched in LWS2 cones; Datasets 1, 2). Correspondingly, GO analysis returned numerous overrepresented categories, including one KEGG category related to proteasomal function (Supplementary Figure 2A) The “RFP” in lws:PAC(H) is dsRedExpress (Tsujimura et al., 2010), which has been noted to mis-fold and/or aggregate (Stenkamp et al., 2021), and engage cell stress pathways (Zhou et al., 2011), and so this is a potential explanation for the enriched presence of proteasome components.

We therefore used an additional FACS-RNA-Seq approach and analysis, to validate our findings and identify more transcripts enriched in LWS1 vs. LWS2 cones. Both types of LWS cones were sorted from thrb2:tdTomato transgenic retinas, using the strategy reported in Sun et al. (2018b); (Supplementary Figures 2B, C). Subsequent RNA-Seq (GEO upload in progress) and DE analyses using the thrb2:tdTomato transcripts vs. the lws1:GFP or lws2:RFP transcripts returned lists of transcripts DE in LWS1 cones vs. all LWS cones, and in LWS2 cones vs. all LWS cones (Dataset 3 and Supplementary Figure 2D). Not surprisingly, the list of transcripts enriched in LWS2 cones vs. all LWS cones was dominated by components of the proteasome; however, both lists contained transcripts identified in the prior analysis as DE, and also returned some novel findings. The additional dataset is also publicly available (GEO accession # GSE232902) (see text footnote 3).

As a further means to assess transcriptional distinctions between LWS1 and LWS2 cones, we used scRNA-Seq of WT adult zebrafish retinas. TSNE plotting identified 16 distinct clusters from dissociated retinal cells (Figure 2A), with cone identity assigned to cluster#5, based upon expression of known cone transcripts: opsin markers, gnat2, and pde6c. Within this cluster, individual cells could be identified based upon expression of opn1lw1 (94 cells) or opn1lw2 (199 cells) (Figure 2B). These cells were localized in the TSNE space near each other, and some lws2 transcript was present in some LWS1 cones, suggesting that cones co-expressing lws1 and lws2 were sampled in this study (Figure 2). All transcripts identified within these cones included 573 in the lws1 + cells (31 unique to lws1 +), and 886 in the lws2 + cones (13 unique to lws2 +) (Dataset 4). The scRNA-Seq dataset is publicly accessible (GEO accession #GSE234661).⁴

A DE analysis of lws1 + cells vs. lws1- cells of the scRNA-Seq results provided another means of identifying transcripts enriched in LWS1 cones vs. other retinal cell types (Dataset 4). Transcripts enriched in LWS1 cones included many identified by the bulk RNA-Seq approach, such as gngt2a, arr3a, and ablim3, as well as those not in the bulk RNA-Seq DE lists, including aanat21 and si:busm1-57f23.11 (Figure 2C and Dataset 4). Transcripts enriched in LWS2 cones vs. other retinal cell types included gngt2b, thrb1, and six7. The most consistently identified DE transcripts in LWS1 vs. LWS2 cones, using both approaches, were the two paralogs encoding γ subunits of cone transducin, with gngt2a enriched in LWS1 cones, and gngt2b enriched in LWS2 cones. Known expression patterns of these paralogs appeared to support some degree of cone subtype specificity, with gngt2a found in ventral/peripheral retina, and gngt2b found in dorsal/central larval retina (Lagman et al., 2015), patterns shown to reflect those of lws1 and lws2, respectively (Takechi and Kawamura, 2005; Tsujimura et al., 2010; Ogawa and Corbo, 2021). The mapping of gngt2a and gngt2b paralogs on the TSNE graph also suggested coordinated co-expression of these transcripts within specific LWS cone subtypes, although the gngt2s were also more broadly associated with other retinal cell type clusters (Figure 2B).

We wished to further evaluate potential DE transcripts through curating a “short list” to prioritize for additional study (Table 3). Transcripts were prioritized based upon (1) appearing in more than one DE list; (2) evidence in the literature (or in the ZFIN database) of expression patterns consistent with DE in LWS1 vs. LWS2 cones; (3) levels of expression in the bulk RNA-Seq dataset suggesting the transcript would be detectable by in situ hybridization; (4) potential cone-specific functional relevance (phototransduction, circadian rhythm, cell adhesion); and (5) potential relevance in regulation of lws1 vs. lws2 expression (nuclear hormone signaling) (Mitchell et al., 2015; Mackin et al., 2019). We noted that si:busm1-57f23.1 within our short list, was also detected as significantly downregulated in a zebrafish thrb mutant (Volkov et al., 2020). We reasoned that other DE transcripts identified in thrb-/- vs. WT in this previous study may also be DE in LWS1 vs. LWS2 cones, and so added these to the short list.

We selected eight of these transcripts (Table 2; LWS1-enriched: gngt2a, nrip1a, vax1, vax2, si:busm1, cry3a. LWS2-enriched: gngt2b, nr2f2) to test two hypotheses: (1) LWS1 and LWS2 cones are functionally distinct (already supported by the outcomes of the RNA-Seq analyses); and (2) multiple genes are coordinately regulated in LWS cones by thyroid hormone. We also wished to use these transcripts to evaluate heterogeneity within the LWS cone subtypes, as such heterogeneity was noted by Aramaki et al. (2022) for M opsin dominant vs. S opsin dominant cone types in mouse. Our initial approach was to focus upon the first hypothesis by evaluating expression patterns of these transcripts in whole adult wild-type retinas using multiplex in situs.

We next used a TH treatment protocol demonstrated to increase lws1 at the expense of lws2, in individual cones (Mackin et al., 2019), to evaluate the response of the eight transcripts to TH, thereby testing our second hypothesis. This protocol involved treatment of zebrafish embryos at 48 hpf with 100 nM T3, or the DMSO vehicle, and collecting whole larvae at 96 hpf for measurement of relative abundance of transcript (qPCR) and changes in expression pattern (multiplex fluorescence in situ) as the experimental endpoints.

In total, the bulk RNA-Seq results showed 95 transcripts enriched in LWS1 cones and 186 transcripts enriched in LWS2 cones, a finding that supports our first hypothesis that these cone subtypes differ beyond the level of opsin expression. To our knowledge, any such distinctions in the human LWS vs. MWS cone populations, other than opsin expression, have not been noted, although this is largely due to the challenges of detecting LWS-expressing vs. MWS-expressing cones within a dataset (Lukowski et al., 2019). The tandemly-replicated LWS/MWS opsin genes of primates display 98% homology at the level of transcript (Onishi et al., 2002) and are difficult to distinguish using standard RNA-Seq approaches. Interestingly, the scRNA-Seq approach of Peng et al. (2019) permitted this analysis for macaque cones, with the conclusion that these LWS and MWS cones are transcriptionally distinct only for the opsin genes, in contrast to our findings for the zebrafish.

Some of the transcripts enriched in zebrafish LWS1 vs. LWS2 cones suggest the possibility of further functional differences between the LWS cone subtypes. For example, the presence of the paralogous gamma transducin enriched in each cone subtype may reflect differences in phototransduction kinetics and/or recovery. However, to our knowledge such distinctions have not been experimentally tested. In addition, because multiple transcripts encoding factors involved in circadian rhythms were found to be enriched in LWS1 cones [cry1bb, per3, (Dataset 1) cry1ba, aanat2, per2 (Dataset 4)], and no circadian related transcripts were enriched in LWS2 cones, LWS1 cones may have specialized functions in circadian rhythm. The majority of LWS1 cones are located within the ventral domain of the retina, exposed to direct sunlight, and therefore are in an ideal position to obtain circadian information (Takechi and Kawamura, 2005; Ogawa and Corbo, 2021). While the transcript tested using HCR in situ did not appear to be ventrally enriched, differences in sample collection timing may have been a factor (near lights-on for RNA-Seq, and midafternoon for adult whole retina collection), and retina-specific spatial expression data for the other circadian genes have not been reported in the literature.

Other differences in gene expression may give insight into the regulatory landscape of the retina. Multiple nuclear receptors and proteins that interact with nuclear receptors were predicted to be DE between LWS1 and LWS2 cones, including nrip1a and nr2f2. Given that ligands of nuclear hormone receptors (retinoic acid and TH) can regulate expression of lws1 vs. lws2 (Mitchell et al., 2015; Mackin et al., 2019), nrip1a and nr2f2 may be considered candidates for participation in this regulation. The expression of vax2 in LWS1 cones also proves interesting. Vax2 is instrumental in regulating RA metabolism by altering the expression of RA-catabolizing and RA-synthesizing enzymes in developing mouse retina (Alfano et al., 2011). RA is known to be important in regulating cone opsin expression in multiple species (Prabhudesai et al., 2005; Roberts et al., 2005) including zebrafish, and specifically for regulating lws1 vs. lws2 (Mitchell et al., 2015), and RA receptors are known to heterodimerize with TH receptors (Berrodin et al., 1992). RA signaling takes place in ventral retina of juvenile zebrafish, spatially coinciding with a “transition zone” where LWS cones switch from lws2 to lws1 expression as the retina grows (Mitchell et al., 2015). Further, experimentally athyroid juvenile zebrafish display only a small ventral patch of lws1 expression, which also spatially coincides with the RA signaling domain (Mackin et al., 2019). As vax2 is known to be present in the zebrafish retina long before LWS cone development (French et al., 2009), this gene may play an upstream role in lws regulation by spatially tuning RA levels which, in turn, impart dorsoventral location information to developing LWS cones along with TH. While the presence of vax2 in larval zebrafish retina is not surprising and likely is involved in the regulation of many genes, vax2 expression in adult zebrafish was unexpected as vax2 is not expressed in adult mouse retina (Barbieri et al., 1999). We speculate that vax2 may be important in the maintenance of correct topography of cone subtypes in the adult zebrafish, perhaps even after regeneration (Stenkamp et al., 2021).

We initially hypothesized that TH could be a master regulator of transcriptional differences between the LWS cone subtypes, LWS1 and LWS2. For this hypothesis to be true, DE genes would be coordinately regulated with lws1 and lws2 such that genes enriched in LWS1 cones would be upregulated by TH and genes enriched in LWS2 cones would be downregulated by TH. Our observations did not match this hypothesis. While gngt2a and gngt2b show some features of coordinated regulation with lws1 and lws2, respectively, other transcripts such as si:busm1 show the opposite–si:busm1 was predicted to be enriched in LWS1 cones but appears downregulated by TH. Nrip1a is enriched in LWS1 cones and is upregulated by TH, similar to lws1, but its expression domain extends beyond the lws1 domain and even beyond the photoreceptor layer. Other transcripts that are specifically enriched in LWS1 cones such as vax1 and vax2 display no transcriptional response to TH treatment. Therefore, TH is likely involved in regulation of transcripts other than lws1 and lws2, but in a more complex manner than originally hypothesized, and perhaps independently of the lws opsin genes.

Transcriptional heterogeneity among photoreceptor subtypes was recently described in detail for the zebrafish as the “partitioning” of expression of paralogs of photoreceptor components that emerged through whole genome duplications (Ogawa and Corbo, 2021). For example, the authors took note of partitioning of a paralog of the gamma subunit of transducin, gngt2a, within LWS1 and RH2-4 cones (Ogawa and Corbo, 2021), which express the most red-shifted members of the tandemly-replicated lws and rh2 cone opsin gene arrays, respectively (Chinen et al., 2003). While our studies are largely consistent with this concept, the patterns of expression of the gngt2 paralogs in cone subtypes appear more nuanced. Some LWS2 cones express gngt2b, but some express gngt2a or possibly both paralogs. Expression of gngt2a in adult retina is not limited to the ventral domain of LWS1 cones. Similarly, some but not all LWS1 cones express vax1. Some but not all LWS2 cones express nr2f2. For each of these LWS cone subtypes, there appears to be a great deal of transcriptional heterogeneity across the expanse of the retina.

Further heterogeneity is revealed within cone responses to treatment with exogenous TH–some but not all of these genes can be regulated by TH. Based on these observations, it appears that TH may be regulating these genes not because it is particularly regulating the entirety of LWS1 and LWS2 cone phenotype but rather because TH is an important regulator of spatial dynamics of gene expression in the zebrafish retina. Recent work in mouse supports this thyroid hormone-mediated spatial gradient regulation. It was found that trβ2, a TH receptor shared by mice and zebrafish, can control the expression and chromatin state of multiple cone genes that are expressed in a dorsoventral gradient in the mouse retina (Aramaki et al., 2022). This influential study also reveals heterogeneity within populations of mouse cones that express both sws1 and mws opsins, in support of the concept that not all cones expressing a particular opsin are identical (Aramaki et al., 2022). The present study extends this concept to the lws1- expressing and lws2-expressing cone populations of the zebrafish. We aim to identify the TH receptor(s) that mediate the effects of TH in the zebrafish, as well as determine whether any are directly interacting with elements on the regulated gene.

While our approach was able to expand our understanding of LWS1 and LWS2 cone biology, we were limited in our ability to detect all of the potential transcriptional variability in these photoreceptors due to the limitations of sensitivity of bulk and scRNA-Seq. The bulk RNA-Seq results were further restricted by the effectiveness of the cell sorting procedure, and by the presence of proteasome-related transcripts in the LWS2 cone samples. Therefore we leveraged the bulk RNA-Seq datasets with information derived from a scRNA-Seq approach. We evaluated eight transcripts for their spatial patterning and TH responses. It is possible that many of the other differentially expressed genes may be regulated by TH. Further, these evaluations were done at the transcript level and do not reveal protein expression and/or localization, which would add another layer of support to the hypothesis that LWS1 and LWS2 cones are functionally distinct beyond opsin expression. In the future, we plan to perform single cell RNA-Seq on control and TH treated retinas to determine widespread effects of TH on the retinal cell transcriptomes. Due to the similarity of the LWS1 and LWS2 cone transcriptomes, LWS cone studies might also benefit from a manual photoreceptor collection technique, as in Angueyra et al. (2023), in order to increase the signal-to-noise ratio in the transcriptome data and reveal more potentially subtle gene expression differences.

Our results build upon our previous studies that showed lws1 and lws2 are regulated by TH (Mackin et al., 2019) by suggesting that multiple genes within the cones expressing these tandemly replicated opsins are non-stochastically regulated, and that these genes may also be regulated by TH. Further, our findings add to the growing literature that shows TH is a major regulator of spatial patterning in the retina and that its role is conserved across multiple species.