Whole-Organism Developmental Expression Profiling Identifies RAB-28 as a Novel Ciliary GTPase Associated with the BBSome and Intraflagellar Transport

Primary cilia are specialised sensory and developmental signalling devices extending from the surface of most eukaryotic cells. Defects in these organelles cause inherited human disorders (ciliopathies) such as retinitis pigmentosa and Bardet-Biedl syndrome (BBS), frequently affecting many physiological and developmental processes across multiple organs. Cilium formation, maintenance and function depend on intracellular transport systems such as intraflagellar transport (IFT), which is driven by kinesin-2 and IFT-dynein motors and regulated by the Bardet-Biedl syndrome (BBS) cargo-adaptor protein complex, or BBSome. To identify new cilium-associated genes, we employed the nematode C. elegans, where ciliogenesis occurs within a short timespan during late embryogenesis when most sensory neurons differentiate. Using whole-organism RNA-Seq libraries, we discovered a signature expression profile highly enriched for transcripts of known ciliary proteins, including FAM-161 (FAM161A orthologue), CCDC-104 (CCDC104), and RPI-1 (RP1/RP1L1), which we confirm are cilium-localised in worms. From a list of 185 candidate ciliary genes, we uncover orthologues of human MAP9, YAP, CCDC149, and RAB28 as conserved cilium-associated components. Further analyses of C. elegans RAB-28, recently associated with autosomal-recessive cone-rod dystrophy, reveal that this small GTPase is exclusively expressed in ciliated neurons where it dynamically associates with IFT trains. Whereas inactive GDP-bound RAB-28 displays no IFT movement and diffuse localisation, GTP-bound (activated) RAB-28 concentrates at the periciliary membrane in a BBSome-dependent manner and undergoes bidirectional IFT. Functional analyses reveal that whilst cilium structure, sensory function and IFT are seemingly normal in a rab-28 null allele, overexpression of predicted GDP or GTP locked variants of RAB-28 perturbs cilium and sensory pore morphogenesis and function. Collectively, our findings present a new approach for identifying ciliary proteins, and unveil RAB28, a GTPase most closely related to the BBS protein RABL4/IFT27, as an IFT-associated cargo with BBSome-dependent cell autonomous and non-autonomous functions at the ciliary base.


Introduction
The cilium is a conserved organelle, inferred to have existed in the last eukaryotic common ancestor (LECA) and now present in most extant protists, as well as all multicellular animals. Motile cilia generate cell movement or fluid flow, whereas non-motile (primary) cilia have evolved as specialised 'antennae' that capture extracellular sensory cues and orchestrate extrinsic signal transduction pathways linked to development (e.g., Sonic hedgehog) [1,2]. Cilium dysfunction in humans is associated with a growing number of so-called ciliopathies that affect virtually all physiological and developmental functions [3]. For example, Bardet-Biedl syndrome (BBS) includes retinal degeneration, cystic kidneys, obesity and skeletal anomalies (polydactyly) as primary ailments [4].
Cilia are subdivided into distinct subcompartments, each with unique structural and functional features, as well as molecular compositions [5]. The canonical cilium of 9 doublet microtubules (MTs) extends from a mother centriole-derived basal body, which connects via distal appendages (transition fibers) to the plasma membrane. The proximal-most 0.2-1.0 μm of the axoneme, called the transition zone, functions in early ciliogenesis, and together with basal body structures provides a permeability barrier that separates the ciliary cytosol and membrane from the cell body [6][7][8]. Additional subregions include the inversin and distal tip compartments, as well as the ciliary pocket, which is a depression of the periciliary membrane where the basal body is rooted [5]. Many ciliopathy proteins and associated complexes localise to particular ciliary subcompartments, where they conduct subdomain-specific functions [5,9].
Cilia rely on various intracellular transport systems to sort and deliver the protein cargo required for cilium formation, maintenance and function [10]. The best understood is intraflagellar transport (IFT), which consists of large macro-molecular assemblies that move bidirectionally between the ciliary base and tip, driven by kinesin-2 anterograde (base to tip) and IFT-dynein retrograde (tip to base) motors [11][12][13]. Associated with the motors-and essential for IFT-are the IFT-A and IFT-B complexes, which likely serve as cargo adaptors [11]. The IFT-associated BBS complex (BBSome) also tethers ciliary cargo and regulates the coupling of IFT-A and IFT-B complexes [14]. Also important are membrane trafficking pathways that regulate vesicle formation and transport between post-Golgi sorting stations and the periciliary membrane, as well as endocytic retrieval and recycling events at the periciliary membrane and within the ciliary pocket [15][16][17][18][19]. Various IFT and ciliary membrane trafficking regulators have been identified, including small GTPases of the RAB, ARF and ARL families, that function during early cilium formation as well as transport events post-ciliogenesis [10,20,21].
C. elegans represents a powerful genetic model for investigating cilium formation and function [22]. Hermaphrodite worms possess 60 ciliated cells (of 960 total), all of which are sensory neurons. The non-motile sensory cilia extend from the dendritic tips and many are contained within bilateral chemo-and thermo-sensory cuticular organs, supported by glial cell (sheath and socket) processes that establish environmentally exposed channels [42][43][44]. C. elegans cilium morphologies range from the canonical rod-like to forked, multi-branched and membrane-expanded structures [22]. Worm cilia also possess ultrastructural features conserved in vertebrate/mammalian cilia; for example, amphid (head) and phasmid (tail) channel cilia possess long A-tubule extensions that establish a proximal axonemal region or 'middle segment' of 9 outer doublet MTs and a 'distal segment' of 9 outer singlet MTs [22]. Because many ciliary genes and pathways are conserved in nematodes, and complete loss of cilia is non-lethal [26], C. elegans has been a leading metazoan model for discovering new ciliary genes and uncovering new insight into ciliary transport, function and disease mechanisms.
In this study, we identified a unique expression profile for ciliary genes using a series of RNA-Seq libraries generated specifically to improve annotation of the transcriptome [45][46][47]. We confirmed that our clustering analysis identifies known ciliary proteins, including several not previously studied in C. elegans, and uncovers novel conserved ciliary proteins. One of these proteins is RAB-28, which is expressed exclusively in ciliated cells, where it associates with the periciliary membrane and behaves as an IFT cargo via BBSome-and nucleotide binding-dependent mechanisms. Overexpression of predicted active or inactive forms of RAB-28 leads to variant-specific ciliary and cell non-autonomous sensory pore morphogenesis defects. Together, our work provides a novel approach to finding new ciliary proteins, and uncovers a functional association between the BBSome, IFT and the orthologue of the cone-rod dystrophy protein, RAB28.

Results
Whole-organism developmental expression profiling reveals a ciliumspecific gene expression pattern To provide a complementary approach to ciliary gene discovery, we took advantage of the temporally-invariant birth of all C. elegans cells and tissues, including ciliated neurons, during development [48] (Fig 1A). Nearly all of the 60 ciliated neuronal cell types in C. elegans hermaphrodites are born within a discrete embryonic time period 300-450 minutes post-fertilisation, with cilium formation occurring very shortly thereafter (Fig 1A). We hypothesised that ciliogenesis genes are highly expressed during this time period and therefore distinguishable from genes required for general neuronal formation and development, which are expressed during a broader time window (Fig 1A). Using an available whole-organism developmental series of RNA-seq libraries from C. elegans [45][46][47], we confirmed this hypothesis: many wellcharacterised cilia genes are highly expressed in the early embryo, display peak expression in the late embryo and first larval stage, and show greatly reduced expression during subsequent developmental stages (Fig 1B).
Next we sought to identify novel ciliary/ciliogenic genes displaying a similar cilia-related expression pattern. Using a set of 41 well-characterised ciliary component genes as "baits", representing the cilia-related gene expression profile during development (sheet 1 in S1 Table), we queried our ciliary transcriptome to identify other gene "preys" with similar expression profiles across the RNA-Seq libraries (Fig 1B and sheet 2 in S1 Table). Hierarchical clustering of genes based on temporal expression reveals a tight cluster of 34 of the bait genes with 151 prey genes (cluster 1), many of which are uncharacterised (S1 Fig and sheet 3 in S1 Table). We also filtered the gene list to only include those with human orthologues [49]. To validate our dataset, we determined if filtered cluster 1 is enriched for genes with >12 hits in the ciliary database Cildb [23,24] (sheets 3 and 4 in S1 Table). We found that conserved Cildb-represented genes are significantly enriched (84.5 fold) in filtered cluster 1 compared to the entire genome (p<0.0001, chi-squared test; sheet 4 in S1 Table). We further compared the cluster 1 gene set to two previous expression studies [25,50] and see significant overlap between lists (p<0.0001; sheet 4 in S1 Table).
Because our genes were identified by shared temporal gene expression, we reasoned that they may share common promoter regulatory elements beyond the previously identified Xbox motif [27,28] and so we performed a motif discovery analysis on the conserved genes in cluster 1 (sheet 3 in S1 Table). Not surprisingly, we find several motifs that are significantly enriched, including matches to the GAGA factor binding site (S2 Fig). As might be expected, the 14-bp X-box was not identified by the motif elucidation program despite its enrichment in our dataset (sheet 4 in S1 Table), likely due to the number of degenerate bases in its consensus sequence and because the program identifies motifs only up to 8-bp long [26][27][28]51]. None of the elements we identified have been previously associated with cilia gene expression.

Post-ciliogenesis chemoreceptor co-expression profiling
We observed that chemoreceptors (G protein-coupled receptors or GPCRs) [52], many of which are known to be cilium-localised, are largely absent from our ciliary dataset (sheet 3 in S1 Table). We hypothesised that such genes might be transcribed post-ciliogenesis. Using a statistical predictive modelling strategy similar to that used above, we identified 80 genes sharing the expression profile of srg-36, a cilium-localised dauer pheromone receptor required for entry into the dauer larval diapause life stage (Fig 1C and sheet 5 in S1 Table) [53]. Of the 80 genes, 27 are serpentine transmembrane receptors, representing a significant 7-fold enrichment over all such genes in the genome (p<0.0001, chi-squared test). Analysis of the upstream regulatory regions revealed that 25 of the 28 chemoreceptor genes (including srg-36) possess an E-box regulatory element (S3 Fig) [54]. This suggests that their similar temporal expression profiles stem from regulation by the same, or similar, transcription factors such as the E-boxbinding Basic Helix-Loop-Helix bHLH transcription factor [54]. Interestingly, 4 of the predicted serpentine receptors, including srg-36, are known to be expressed in ciliated sensory  [53,55,56]. We note that the rise and fall of chemosensory/GPCR gene expression trails that of the cilia-related genes (Fig 1C), indicating that these two expression profiles of related developmental processes are temporally discrete.
All reporters employing the endogenous gene promoter show enriched or almost exclusive expression within ciliated cells (S4A Fig), thus validating the predictions from our coexpression profiling. Specifically, maph-9, rpi-1, ccdc-104, rab-28, and fam-161 are exclusively expressed in ciliated neurons, with ccdc-149 expressed in most ciliated neurons as well as pharyngeal neurons, touch receptor neurons, and motor neurons as previously reported (Figs 2 and S4A) [59]. Furthermore, the translational reporters reveal that all 7 proteins localise at the base of, or within, cilia (Fig 2 and Fig 3A). FAM-161::GFP signals are found at the ciliary base and within the proximal part of the ciliary axoneme, including the transition zone (TZ) (Fig  2). This localisation is similar to that of the mammalian protein implicated in retinitis pigmentosa, and is consistent with the suggested role of FAM161b in cargo delivery to photoreceptor outer segments (cilia) [60,61]. Both GFP reporters for RAB-28 (driven by bbs-8 gene or endogenous promoter), the orthologue of the small ciliary GTPase RAB28 linked to human autosomal-recessive cone-rod dystrophy [62][63][64], localise along the entire cilium. (Fig 2 and Fig 3A). Ciliary axoneme localisations are also observed for the GFP-tagged RPI-1, MAPH-9, YAP-1 and CCDC-104 translational reporters; in contrast, CCDC-149::GFP is absent from the ciliary axoneme, although a pool of signal is evident at the ciliary base (Fig 2). The ciliary localisation of CCDC-104::GFP is consistent with the mammalian ciliary base and axonemal localisations reported for the ARL3 interacting protein, CCDC104/BARTL1 [65]. We also found that RPI-1::GFP appears to localise to ciliary and cytoplasmic microtubules, similar to its mammalian development, hatching and post-embryonic development are highlighted above the x-axis with the worm stages illustrated below. Most ciliated neurons are born in late embryogenesis with ciliogenesis occurring shortly after these cells are born. Worm illustrations modified from WormAtlas (www.wormatlas.org). Q/G/K/P/H/V/T/M refer to different embryonic cell lineages (http://www.wormatlas.org/celllineages.html). (B) Normalised expression levels for all 'bait' cilium-associated genes (red) and the top 276 target genes (black, p-value < 5e-7) over the developmental series of analysed RNA-Seq libraries. Expression peaks for all genes during late embryo or larval stage 1, which matches closely the birth of the ciliated neurons shown in (A). This cilia-specific gene expression pattern shows much lower levels of expression from larval stage 2 to young adult. (C) Box and whisker plots of the normalised expression across the whole-organism developmental RNA-Seq libraries for the chemosensory and ciliary gene sets. Expression of the ciliary set immediately precedes that of the chemosensory set. Peak expression for ciliary genes is at the late embryo-L1 larval stage, whereas expression peaks at L1-L2 larval stages for chemosensory genes. EE; early embryo. LE; late embryo. L1; larval stage 1. L2; larval stage 2. L3; larval stage 3. L4; larval stage 4. YA; young adult. orthologue (RP1) implicated in retinitis pigmentosa (Fig 2) [66]. Our finding of ciliary localisations for YAP-1, MAPH-9 and CCDC-149 are the first in any cellular system or organism, and is consistent with a known ciliogenesis role for mammalian YAP, and cilia-related phenotypes in MAP9 (MAPH9 orthologue)-disrupted zebrafish and dachshunds [67][68][69][70]. YAP is regulated by the ciliary protein NPHP4, and MAP9 is phosphorylated by the cilia-disassembly/centrosomal kinase PLK1 [71,72].
We also examined genes excluded from filtered cluster 1 due to a reported lack of predicted human orthologues (according to the Ortholist database; [49]). Highly ranked in unfiltered Candidate cilium-associated proteins localise to ciliary structures Representative images of amphid (head) and phasmid (tail) cilia from N2 wild type worms co-expressing GFP-tagged 'translational' reporters for candidate ciliary genes and tdTomato-tagged XBX-1 (IFT protein that localises to the ciliary base and axoneme). All reporters driven by the endogenous gene promoter, except for YAP-1 and RAB-28, where a bbs-8 gene promoter sequence was used (see Fig 3A for a GFP::RAB-28 reporter driven by the endogenous promoter). GFP tags are on the Nor C-terminus, as indicated above the panels. PCMC (periciliary membrane compartment); distal dendrite (den) swelling enriched for endocytosisassociated proteins and vesicles that regulate ciliary membrane homeostasis [19]. The transition zone (TZ) extends from basal body (BB) and functions as a ciliary gate that regulates protein entry and exit to and from cilia [6]. Middle (MS) and distal (DS) segments; characterised by a circular array of 9 doublet (A and B tubules) and 9 singlet (A tubules) microtubules, respectively (not shown). Note that for GFP::CCDC-149, signals are observed at the amphid/phasmid basal bodies (arrows) and as punctae (asterisks) along the dendrite of the OLQ ciliated neuron running parallel to the amphid neurons. Scale bars; 5 μm.
doi:10.1371/journal.pgen.1006469.g002 cluster 1 is an uncharacterised gene with an X-box promoter element [26], tza-3, which we demonstrate is expressed exclusively in ciliated neurons, and encodes a ciliary TZ protein (S4B Fig). Whilst TZA-3 is not conserved, the TZ localisation of this predicted transmembrane protein depends on the highly conserved core TZ scaffolding protein MKS-5 and "MKS module" components MKSR-1 and MKSR-2, but not the "NPHP module" protein NPHP-4 (S4B Fig). These findings are consistent with TZA-3 associating with the MKS module, which plays a role in ciliary gating [6][7][8]. Thus, our expression profiling approach identifies non-conserved ciliary proteins in addition to conserved ciliary proteins.
Of the novel C. elegans ciliary proteins we identified, we sought to characterise RAB-28 further for the following reasons: (1) human RAB28 is associated with a possible ciliopathy (conerod dystrophy [62,63]), (2) many small GTPases play essential roles in ciliogenesis and ciliary membrane trafficking pathways [10], and (3) aside from reports of RAB28 localising at the basal body [62] or the cilium [64], there is a complete absence of molecular studies on this evolutionarily conserved protein that support a functional and mechanistic link to cilia.

RAB-28 undergoes intraflagellar transport (IFT)
We further investigated the ciliary localisation of GFP::RAB-28 using time-lapse imaging and found that this GTPase undergoes continuous, bidirectional, IFT-like movement along ciliary axonemes (S1 Movie). In contrast, the additional non-ciliary GFP::RAB-28 signals present throughout the neurons are diffuse, with no processive movement detected. In C. elegans amphid and phasmid cilia, anterograde IFT is driven by two kinesin-2 motors (kinesin-II and OSM-3) to yield distinct anterograde rates along the middle and distal segments [14,73]. Kymograph analysis of phasmid (tail) cilia confirm that ciliary GFP::RAB-28 moves at IFTassociated velocities, displaying characteristic average anterograde rates of~0.7 μm/sec (middle segment; proximal part of axoneme) and~1.2 μm/sec (distal segment; distal part of axoneme), and retrograde rates of~1.4 μm/sec (along entire cilium length) (Fig 3A and 3B and S1 Movie) [73]. To examine if this ciliary trafficking is truly associated with IFT, we examined GFP::RAB-28 in worms with disrupted CHE-11 (IFT140), which is a component of IFT-A essential for IFT [28,74,75]. Although GFP::RAB-28 is observed within the truncated cilia of che-11 mutants, we could not detect processive movement of the GFP signals, indicating that the trafficking behaviour of RAB-28 within cilia is bona fide IFT (Fig 3A and S1 Movie). In contrast to observations for other IFT-associated proteins in IFT deficient-worms [28,74,75], GFP::RAB-28 does not abnormally accumulate at the ciliary base or tip, or along the axoneme of these animals ( Fig 3A). This lack of accumulation, despite the active transport defect, suggests that our GFP::RAB-28 reporter also freely diffuses within (and between) the ciliary and dendritic compartments, in addition to undergoing IFT. We confirmed this hypothesis using a fluorescence recovery after photobleaching assay in wild type and che-11 mutant worms, which show that GFP::RAB-28 is highly mobile, displaying very rapid exchange kinetics between ciliary and dendritic pools ( Fig 3C). Thus, in addition to active transport via IFT, GFP::RAB-28 undergoes IFT-independent free diffusion, which explains why this reporter does not abnormally accumulate within che-11 mutant cilia.

RAB-28 is an IFT cargo
To determine if RAB-28 is required for IFT, we examined ciliary structure, function and protein transport in worms containing a deletion (gk1040) in rab-28. The gk1040 mutation removes the GTP-binding switch II domain and the farnesylated C-terminal CAAX motif (S5A Fig), both of which are critical for RAB protein function [76]; thus, gk1040 is likely a severe loss-of-function or null allele. We outcrossed the gk1040 background at least 3 times with wild type worms to remove background mutations unlinked to gk1040. We find that ciliated amphid and phasmid sensory neurons in rab-28 mutant worms display a normal dye filling response, which suggests that ciliary structures are intact (short cilia usually abrogate dye uptake due to lack of environmental exposure) [ Thus, unlike the disruption of IFT proteins, loss of RAB-28 does not affect ciliary structure, transport, or function, at least for those proteins and cilia that were analysed. We conclude, therefore, that RAB-28 behaves more like an IFT cargo, or peripherally-associated component of IFT complexes, rather than a central component of the machinery required for bidirectional movement (IFT).
RAB-28 association with the periciliary membrane and IFT is dependent on GTP-binding and the BBSome subunit orthologue, BBS-8 Next, we assessed the requirements of GDP and GTP nucleotide binding for RAB-28 localisation and transport. We made GFP-tagged constructs containing T49N or Q95L mutations in RAB-28, which are predicted to trap the GTPase in GDP (inactive) or GTP (active)-locked states, respectively [81]. These constructs were injected into wild type worms at the same concentration (5 ng/μl) as the GFP::RAB-28(WT) construct described above (Fig 3), and all appear to be expressed at similar levels. Like the RAB-28(WT) reporter, RAB-28(T49N) (hereafter termed RAB-28(GDP)) displays diffuse localisation throughout the entire neuron, including the cilium and its associated periciliary membrane compartment (PCMC) (Fig 4A). However, unlike wild type RAB-28, RAB-28(GDP) fails to undergo detectable IFT, indicating a highly reduced ability to associate with IFT trains (Fig 4A and S1 Movie). In contrast, the localisation of the RAB-28(Q95L) (hereafter termed RAB-28(GTP)) is highly enriched at the PCMC compared with the rest of the cell, forming a 'tunnel-like' localisation indicative of association with the periciliary membrane ( Fig 4A). RAB-28(GTP) also undergoes clearly detectable IFT, albeit at an apparently reduced frequency compared to the wild type reporter (Fig 4A and S1 Movie). Notably, the localisation and IFT behaviour of all three GFP::RAB-28 reporters (WT, GDPlocked, GTP-locked) is not altered when expressed in the rab-28(gk1040) null background (S6  the higher expressing marker, and indicating that the periciliary membrane is the primary site of activated RAB-28 function in nematode sensory neurons (Fig 4A). Taken together, these data show that GTP binding targets RAB-28 to the periciliary membrane and facilitates its association with IFT trains.
These observations for C. elegans RAB-28 are reminiscent of those for mammalian RAB8, whose ciliary targeting also depends on GDP-GTP exchange, via a mechanism involving RABIN8 (RAB8 guanine nucleotide exchange factor) and BBSome functions at the base of cilia [82][83][84]. We therefore examined the localisation and transport behaviour of all three GFPtagged RAB-28 reporters (WT, GDP-locked and GTP-locked) in a likely null allele (nx77) of a C. elegans BBSome gene orthologue, bbs-8. The bbs-8(nx77) mutation does not alter the ciliary targeting and distributions of RAB-28(WT) and RAB-28(GDP) compared to a wild type background ( Fig 4B). In contrast, BBS-8 loss has a striking effect on both of our RAB-28(GTP) reporters, preventing their accumulation at periciliary membranes ( Fig 4B and S7 Fig). Furthermore, the IFT movement of RAB-28(WT) and RAB-28(GTP) is abolished in bbs-8 mutant cilia (Fig 4B and S1 Movie). From these data, we conclude that the IFT-associated BBSome is required for targeting activated RAB-28 (GTP-locked) to periciliary membranes and loading onto IFT trains.

Overexpression of constitutively active or inactive RAB-28 induces distinct cell autonomous and cell-non autonomous sensory pore defects
The differences in localisation and IFT behaviour between the constitutively inactive (GDPlocked) and active (GTP-locked) forms of RAB-28 prompted us to investigate cilium structure and function in worms overexpressing these small GTPase variants. Assessment of dye-filling reveals reduced levels of dye uptake in the amphid neurons of wild type worms expressing RAB-28(GTP), compared to normal dye-filling for wild type worms expressing RAB-28(WT) or RAB-28(GDP) (Fig 5A). Similar dye-filling results were obtained when these variants were expressed in rab-28(gk1040) worms, although worms expressing RAB-28(GDP) may display a very weakly penetrant phenotypic defect (S8A Fig). Despite the distinct dye-filling phenotypes, worms expressing RAB-28(GTP) or RAB-28(GDP) both show a reduced roaming behaviour (Fig 5B), suggestive of cilia-related sensory abnormalities [78].
To investigate the structural basis of these defects, we examined amphid pore ultrastructure in wild type (N2) worms expressing GFP-tagged RAB-28(GTP) or RAB-28(GDP). Somewhat surprisingly, despite the dye uptake phenotype, RAB-28(GTP)-expressing amphid channel cilia are of normal length and morphology, with intact middle and distal segments, as well as TZ compartments (Fig 5C). However, an increased number of ciliary axonemes possess B-tubule seam breaks ('unzipped' microtubules), indicative of incomplete protofilament formation (S8B Fig). A more striking phenotype is that the proximal region of the sheath cell-defined amphid channel is greatly enlarged, and the ciliary axonemes are less tightly bundled (Fig 5C). Also, the channels of RAB-28(GTP) expressing worms display reduced electron densities, potentially indicating abnormally low concentrations of matrix material, which is thought to be secreted by the sheath cell (Fig 5C) [43]. In contrast, although the amphid channel in RAB-28(GDP)expressing worms is not expanded, these worms exhibit a number of abnormally large dense matrix-filled vesicles (MFVs) in the sheath, in that portion of the cell surrounding the channel (Fig 5C). Such electron dense vesicles are not found in the channel region of the amphid sheath cell of wild type (N2), rab-8(gk1040) or RAB-28(GTP) overexpressing worms (Fig 5C), indicating this phenotype is specific to RAB-28(GDP) overexpression. A second phenotype unique to RAB-28(GDP)-expressing worms is that one channel axoneme is missing from the amphid pore, and instead lies adjacent to the channel, embedded in the sheath (Fig 5C). Thus, we conclude from these overexpression data that RAB-28 serves ciliated cell autonomous and sheath cell non-autonomous roles in sensory pore formation and function.

Discussion
In this study, we identified a specific C. elegans ciliary gene expression profile using a set of RNA-Seq libraries, originally produced to annotate the transcriptome, that together provide a temporal landscape of whole-organism gene expression [45][46][47]. A major advantage of C. elegans as a model for ciliary biology is that ciliogenesis occurs during a temporally discrete window of embryonic development (Fig 1) [48], making it well suited for this type of temporal gene expression analysis. The expression profile for known cilia genes peaks at the onset of ciliogenesis and decreases post-ciliogenesis, likely to levels adequate to maintain cilium structure. Using this pattern, we uncover new ciliary associations for several conserved, but poorly characterised proteins, including the orthologues of YAP, MAP9 and CCDC149. The ciliary localisation of MAPH-9 (microtubule-associated protein 9) is consistent with a known functional association between MAP9 and the centrosomal regulatory kinase PLK1 [72], the centrosomal localisation and suggested ciliary functions for MAP9 in zebrafish [69], and a role for MAP9 as a modifier of retinal degeneration in Dachshunds [70]. YAP also has previous links to cilia, including a role in ciliogenesis [67,68] as well as being regulated by the ciliary gate protein, NPHP4 [71]. Although we do not understand the exact ciliary functions of MAP9 and YAP, their suggested roles in other contexts include regulation of ciliary microtubules (MAP9) and cilium-based signalling (YAP) [67][68][69][70][71][72]. Since 64 of the 185 genes in our predicted list of ciliary genes (cluster 1) possess known ciliary associations, we are confident that additional cilia-related components can be uncovered from this dataset.
One of the highest ranked hits in our candidate list of ciliary genes is RAB-28, a small GTPase whose human counterpart is linked to autosomal-recessive cone-rod dystrophy and vision impairment [62,63]. Although recent studies show that mammalian RAB28 is localised to the mammalian photoreceptor basal body and cilium [62,64], its role in photoreceptors and in other tissues and organisms remains unknown. Here, we show that C. elegans RAB-28 is specifically expressed in ciliated sensory neurons, and undergoes IFT. RAB-28 IFT association depends on its GTPase activity; whereas the GDP-locked inactive form localises diffusely with little or no observable IFT, the GTP-locked active form concentrates almost exclusively at the periciliary membrane and undergoes IFT. Moreover, the periciliary membrane association and IFT behaviour of active RAB-28 depends on the BBSome, which functions as an IFT cargo-adaptor [11,14,85]. We conclude that C. elegans BBSome complexes recruit activated RAB-28 to the periciliary membrane, and regulate docking to IFT trains (Fig 6). Our findings also demonstrate that RAB-28 is probably not a 'core' component of the BBSome or IFT particles, since ciliary structures, as well as BBSome and IFT-subcomplex A/B protein localisations, are unaffected in rab-28 mutant worms. Instead, RAB-28 may associate with these complexes as part of a 'peripheral' submodule, performing auxiliary functions unrelated to ciliogenesis, such as those described for metazoan IFT27 (Rab-like 4; RABL4), ARL6 (BBS3) and RAB8, whose ciliary targeting or removal also depend on GTP binding and the BBSome (ARL6, RAB8) [86][87][88][89]. Consistent with a non-essential global requirement for ciliogenesis, IFT and BBSome integrity, patients carrying predicted strong loss-of-function or possibly null alleles in RAB28 do not present with wider ciliopathy symptoms beyond cone-rod dystrophy [62,63].

RAB-28 cell non-autonomous function
Another revealing observation from our work is that overexpression of constitutively active (GTP-locked) or inactive (GDP-locked) RAB-28 cause variant-specific defects in the non-ciliated amphid sheath cell that supports the ciliated endings of amphid neurons. Specifically, RAB-28(GTP) overexpression causes an enlargement of the sheath cell-defined portion of the amphid sensory pore, whereas RAB-28(GDP) overexpression causes the accumulation of abnormally  dense matrix-filled vesicles (MFVs) in the sheath cell. The RAB-28 variant constructs, as well as the RAB-28(WT) construct, are expressed at similar levels, using the same promoter sequence; thus, the observed variant-specific phenotypes are not likely due to differences in expression level, but instead are linked to the nucleotide bound state of RAB-28. Most interestingly, the abnormal sheath cell phenotypes appear to reflect a cell-non autonomous effect because RAB-28 is expressed in ciliated sensory neurons, and not in the amphid sheath cell (S9 Fig). In C. elegans, the amphid channel is fashioned from glial cell (sheath and socket) processes which extend to the nose tip in close proximity to the dendritic processes of the ciliated sensory neurons. The proximal part of the channel forms from a hole in the sheath process which the channel cilia fully penetrate, whereas the distal part of the channel derives from a doughnut-shaped ending of the socket process that fuses with the cuticle and the underlying sheath (Fig 5C). The formation and function of neurons and glia in the amphid pore appear tightly linked, with channel morphogenesis thought to depend on signals from both cell types [43,44,90]. It is also suggested that amphid channel size is regulated via a balance of 'exocytic' membrane delivery via MFVs (increases channel size) and 'endocytic' membrane retrieval (reduces channel size) at the sheath cell plasma membrane [91]. By integrating this model with our findings, we hypothesise that constitutively active RAB-28 may enhance MFV delivery to the sheath cell membrane, leading to an enlarged channel, whereas constitutively inactive RAB-28 blocks MFV delivery and results in abnormally dense MFVs in the sheath. One prediction of this model is that expression of constitutively inactive RAB-28 would reduce amphid channel size, and consistent with this idea, one channel axoneme fails to enter the amphid pore of RAB-28(GDP) expressing worms.
How RAB-28 expressed in ciliated amphid neurons might control morphogenesis events in the supporting non-ciliated sheath cell is unclear. One scenario is that RAB-28 regulates extracellular release of neuronal factors purported to signal to the sheath cell during channel morphogenesis in the embryo [42][43][44] (Fig 6). In support of this notion, activated RAB-28 is almost exclusively concentrated at the periciliary membrane, which is the site of ectosome release in certain nematode sensory neurons; also, RAB-28 expression is elevated in nematode ectosome-releasing cells [92,93]. Future work will be required to investigate possible roles for RAB-28 in cilia-related ectosome release, and to identify sheath cell morphogenic factors released by the ciliated neurons.
A surprising observation is that sensory pore structure and function appears grossly normal in the rab-28 deletion mutant, which is likely a null allele. This indicates that other genes or pathways can compensate for rab-28 loss-of-function, but that the phenotypic effects of the overexpressed dominant active and inactive RAB-28 variants are more severe and cannot be as easily compensated. One possibility is that the cilium structure, and cilium-dependent behavioural phenotypes in worms overexpressing RAB-28(GTP) and/or RAB-28(GDP) derive from excessive activation of downstream effectors or possibly aberrant sequestration of GTPase activating proteins (GAPs) or GDP-GTP exchange proteins (GEFs) common to RAB-28 and other ciliary GTPases. To fully decipher the mechanistic basis of the RAB-28 gain vs loss of function phenotypes, and the physiological relevance of our findings, future efforts should include an examination of worms expressing RAB-28 variants expressed at endogenous or near-endogenous levels and the identification of currently unknown regulators of RAB-28 (GEF and GAP proteins).

Possible functional relationship between RAB28 and IFT27?
Given the close sequence similarity between RAB28 and IFT27, and some of their common functional properties discussed below, we suggest that these two RAB GTPases may have similar or partially overlapping ciliary functions. RAB28 is conserved across metazoans and many ciliated protists (S10 Fig), and inferred to be one of 23 founding RAB family members present in the last eukaryotic common ancestor (LECA) [94,95]. Phylogenetically, RAB28 is most closely related to IFT27 (RABL4), which likely also existed in LECA and is dispersed throughout the eukaryotic lineage [94,95] (S10 Fig). Like the majority of ciliary and IFT components, both RAB28 and IFT27 are restricted to organisms that have cilia (S10 Fig). Interestingly, although RAB28 and IFT27 are present in most metazoans and many protists (e.g., Chlamydomonas and Trypanosoma have both), neither protein appears to be present in Drosophila and IFT27 is missing from C. elegans (S10 Fig). RAB28 and IFT27 appear to share overlapping functional properties. As we show for nematode RAB-28, Trypanosome and mammalian IFT27 associate with IFT in a GTP-dependent manner [87,88,96]. Also, cilium formation is unaffected in mammalian cells with null mutations in IFT27 (similar to C. elegans rab-28 mutants), although disruption of IFT27 in protists results in short flagella [74,96]. In addition, mutations in IFT27 cause BBS, displaying a broad spectrum of ciliopathy phenotypes, including the cardinal retinal degeneration phenotype observed in RAB28 patients [97,98]. Furthermore, both small GTPases are functionally associated with the BBSome, albeit in different ways. The localisation and IFT motility of nematode RAB-28 is BBSome-dependent, but not vice versa (i.e., BBS-5 localisation is normal in rab-28 mutants), whereas mammalian IFT27 regulates BBSome assembly and removal from the cilium. When taken together, the above observations could suggest that in vertebrates RAB28 and IFT27 function at distinct steps of a common pathway, potentially 'upstream' (IFT27) and 'downstream' (RAB-28) of the BBSome [87,88]. Since IFT27 and the BBSome remove proteins from cilia [85,87,88,99], RAB28 may perform a similar function by facilitating the BBSomemediated removal of a more-restricted set of cargo.
In further support of a common, yet tissue-specific, transport role for related ciliary GTPases such as IFT27 and RAB28, it is noteworthy that disruption of the short ARL6/BBS3 isoform in mice and zebrafish results in broad spectrum BBS-related phenotypes, whereas disruption of the long isoform causes a retina-specific phenotype [100]. Transition zone-localised RPGRIP1 and RPGRIP1L represent additional examples of very closely related ciliopathy proteins linked to a broad tissue pathogenesis (RPGRIP1L), versus a retinal-specific disease (RPGRIP1) [101][102][103][104]. Finally, whereas IFT27 is known to regulate the ciliary localisation of sonic hedgehog proteins [87,88], the identity of RAB28-associated cargoes is currently unknown. Interestingly, RAB28 itself may be a cargo of the prenyl-binding protein PDE6D (cGMP phosphodiesterase delta subunit) which immunoprecipitates RAB28 [105]; in agreement with this possibility, loss of PDE6D in the mouse phenocopies the photoreceptor degeneration of RAB28 patients [106]. In light of our worm data, it is possible that PDE6D delta functions with the BBSome, either together or at distinct parts of a transport pathway, to target RAB28 to ciliary membranes.

Concluding remarks
Our study describes a novel method for identifying genes with ciliary functions that complements other studies (bioinformatics, genomics, proteomics) aimed at deciphering a complete 'ciliome'. From our list of candidates, we uncover several novel, evolutionarily-conserved cilium-associated proteins. Our in-depth analysis of RAB-28 in C. elegans reveals it to be a novel IFT-and BBSome-associated protein with cell autonomous and non-autonomous roles in maintaining the structural and functional integrity of cilia, and of the sensory neuron-supporting cells. Further analyses of mammalian RAB28 are required to reveal its functional similarities and differences with IFT27, and possible functional associations with PDE6D and ARL6, to help explain its seemingly specific function in ciliary photoreceptors, and thus specific association with cone-rod dystrophy.

Materials and Methods
Identifying genes with expression correlated with known ciliary genes We curated a list of 41 C. elegans genes with established roles in ciliary function and used these to search for genes with correlated expression. These "bait" genes included: R148. We used expression data derived from RNA-seq for seven developmental stages (EE = early embryo, LE = late embryo, L1 = mid-L1 larvae, L2 = mid-L2 larvae, L3 = mid-L3 larvae, L4 = mid-L4 larvae and YA = young adult) obtained from [45][46][47]. As the expression values were relative to individual isoforms, we first computed an average RPKM (reads per kilobase per million) value for each gene across all isoforms. Using these 41 genes, we next computed the pairwise Pearson correlation of each gene against all 20,363 C. elegans genes in the data set across all libraries and estimated a P value for each test. As we were interested in capturing genes with expression correlates to any of these known ciliary genes, we retained the smallest P value from these 41 tests for subsequent calculations. To avoid circularity, we retained the second-smallest P value for correlation tests involving one of the bait genes. In addition to direct correlation, we sought to identify genes exhibiting a developmental expression pattern similar to known ciliary genes. We observed that all ciliary genes exhibit at least a 10-fold reduction in expression level (relative to their peak expression level) in the latter three developmental stages. A total of 5960 genes showed this expression trend.
Following the identification of those genes exhibiting strongly correlated expression to any of the known ciliary genes, this set was further clustered to search for sub-groupings of genes with strongly coordinated expression. The distance between individual genes was calculated using the 'dist' function in R and hierarchical clustering was performed using single, average and complete linkage. S1 Fig shows a heat map representation of the expression for each of the candidate ciliary genes with the dendrogram derived from clustering using complete linkage and cluster 1 at the top.

Chemosensory expression pattern and promoter element prediction
Using the method described above, we used srg-36 as a bait to identify co-expressed genes. We found 80 genes (p < 0.001), including 27 other predicted chemoreceptors [52]. The 1000 bp upstream from the start codon of all 28 chemoreceptors (including the bait srg-36) were entered into the program MEME [107,108], which found a consensus that matches the published E-box sequence [54]. Separately, the 1000 bp upstream of all the conserved cluster 1 were entered into the MEME suite program DREME [107,108], which found several potential transcription factor binding sites. We used the MEME suite program TomTom to compare identified sites to known sites [108]. Identified promoter elements were compared to consensus splice and trans-splice sites and any that matched were subsequently removed.

Dye-filling and sensory behavioural assays
For the dye-filling assay [109], worms were incubated in 200 μl of DiI solution (Invitrogen; diluted 1:200 with M9 buffer) for 30 min. After incubation, worms were recovered on seeded nematode growth medium plates for a further 30 min and then mounted on slides. Epifluorescence wide-field imaging under the red filter was used to image DiI uptake into the ciliated amphid and phasmid cells. For the roaming (foraging) assay, single worms were placed for 18 hours onto seeded NGM plates and track coverage assessed using a grid reference [109]. For the osmotic avoidance assay, 5-6 worms were placed within a ring of 8M glycerol (Sigma) supplemented with Bromophenol Blue (Alfa Aesar) on unseeded NGM plates and observed for 10 min. Worms that crossed the barrier were removed from the assay. Chemotaxis attraction assays (towards isoamyl alcohol) were performed as previously reported [109], with a chemotaxis index calculated at 30 min and 60 min.

C. elegans fluorescence microscopy, IFT motility and FRAP assays
Strains in Fig 2 and S4B Fig were anaesthetised with 10 mM levamisole in M9 buffer, mounted on slides with 8% agarose pads, and observed by epifluorescence or spinning-disc (WaveFX; Quorum Technologies) confocal microscopy performed on an Zeiss Axio Observer Z1 with a Hamamatsu 9100 EMCCD camera. Image capture and visualisation were performed on Volocity (PerkinElmer). All other worms were immobilised with 40mM tetramisole (Sigma no. L9756) or microbeads (Polysciences no. 00876-15) and mounted on 4% or 10% agarose pads. Epifluorescence images were taken on an upright Leica DM5000B and confocal images on an inverted Nikon Eclipse Ti microscope with a Yokogawa spinning-disc unit (Andor Revolution). Images were acquired using a charge-coupled device camera (iXon+EM-CCD, Andor Technology) and analysed using Image J software. For IFT assays, time-lapse (multi tiff) movies of IFT along phasmid cilia were taken at 200 ms exposure and 4 fps. Separated anterograde and retrograde kymographs were generated from multi tiff files using Icy image analysis software (http://icy.bioimageanalysis.org/) and rates determined using ImageJ [111]. Fluorescence recovery after photobleaching (FRAP) assays were performed using the above confocal microscope with an attached FRAPPA unit (Andor Technology). Samples were bleached using a single pulse of the 488nm laser at 100% with a dwell time of 100 μs. Images were recorded immediately post-bleach and continuously thereafter. Fluorescence intensities measured with image J software following a recently described protocol [109]. Values were normalised to prebleach values and corrected for signal intensity loss during image acquisition.

Transmission Electron Microscopy
Young adult worms were fixed, sectioned and imaged as previously reported [109]. (B) Representative images of the head and tail regions of N2 wild type and rab-28(gk1040) mutant worms following a DiI incorporation assay into amphid and phasmid neurons. Scale bar; 20 μm. (C) Representative images of amphid and phasmid cilia from N2 wild type and rab-28(gk1040) worms expressing OSM-6::GFP. Scale bars; 2 μm. (D) Transmission electron microscopy images of the amphid pore from serial cross sections of N2 wild type and rab-28(gk1040) worms. Low (large panels) and high (small panels) magnification images are shown. Images representative of at least 4 analysed pores for each strain. Both worms show 10 ciliary axonemes in the amphid pore, with each axoneme consisting of a distal segment (DS), middle segment (MS), transition zone (TZ) and periciliary membrane compartment (PCMC). Cartoon shows the amphid channel in cross section and longitudinal orientations (only 3 of the 10 axonemes shown for simplicity in longitudinal cartoon). Numbers above images indicate the position of the section relative to the most anterior section (at '0'); section positions also indicated in cartoon. Scale bars; 200 nm (large panels); 100 nm (small panels). (E) Assessment of rab-28(gk1040) cilia-related sensory behaviours. Shown is a populationbased isoamyl alcohol (IAA) attraction assay (n = 8 for N2 wild type and osm-5; n = 12 for rab-28), as well as single worm assays that measure roaming (n = 35 for all strains) and osmotic avoidance (n = 86 for N2; n = 37 for osm-5; n = 70 for rab-28) behaviours. % avoidance (fraction of worms that avoid -do not cross-the osmotic barrier) is plotted over the time course of the experiment (600 seconds). Chemoattraction index determined at two time points (30 & 60 minutes). osm-5(p813) worms used as a negative control. Ã p<0.01 (unpaired t-test vs WT control at the appropriate assay time point). Although rab-28 worms present with a slightly reduced osmotic avoidance, this behaviour was not statistically significant compared with WT controls (log-rank and Mantel-Cox survival curve tests). (F) rab-28 (gk1040) mutants are indistinguishable from wild-type worms with regards to the cilia-dependent phenotypes of carbon dioxide avoidance and body size. N2 wild-type and rab-28(gk1040) mutant worms respond statistically indistinguishably to a 4 sec puff of 10% CO 2 delivered at 100s (represented by grey bar) by performing a burst of high-amplitude turns (repeated measures ANOVA, pNS). Error bars represent 95% confidence intervals. N2 and rab-28 (gk1040) mutant worms are statistically indistinguishable with regards to length along the midline of the body (ANOVA, pNS). For both graphs n = 112 worms (divided across 3 plates) for wild-type and n = 182 worms (divided across 3 plates) for the rab-28 (gk1040) mutant. (GTP). Like in a wild type background (Fig 3), all three markers are found in the cilium (cil), with RAB-28(GTP) highly enriched at the periciliary membrane. PCMC; periciliary membrane compartment (p). Den; dendrite Scale bars; 3 μm. (TIF) S7 Fig. A GFP::RAB-28(Q95L) marker expressed at low levels localises exclusively to the periciliary membrane in a BBSome-dependent manner. Representative images of whole phasmid neurons and the ciliary region of amphid neurons from worms with an GFP::RAB-28 (GTP) transgene (oqEx304) expressed at low levels. In the rab-28(gk1040) background, RAB-28(GTP) localises exclusively to the periciliary membrane, with weaker signals in the ciliary axoneme. In bbs-8(nx77) worms, no periciliary membrane enrichment is observed for RAB-28 (GTP). Note that the bbs-8 image has been overexposed to show the diffuse (non-specific) localisation of the faint RAB-28(GTP) signals. PCMC; periciliary membrane compartment. Den; dendrite. Scale bars; 3 um.  [42]. Images show that GFP expression is restricted to ciliated neurons and does not include the sheath cell (outlined by white dotted line). This indicates that RAB-28 is not expressed in the amphid sheath cell. Scale bar; 15μm. (TIF) S10 Fig. Phylogenetic distribution of RAB28 and RABL4 (IFT27) across eukaryotes. RAB28 and IFT27 protein orthologues are only found in eukaryotes that are ciliated during at least one life stage (ciliated cells are shown in green, non-ciliated cells in white). Green and white squares, respectively, denote the presence or absence of RAB28 or IFT27 in a given representative species (as presented in Elias et al. [95]). (TIF) S1 Table. Sheet 1: List of all bait genes used to identify ciliary candidates. Sheet 2: List of all genes with given P and R values. Bait genes are highlighted yellow. Sheet 3: List of genes within filtered cluster 1 from the clustering analysis of cilia-related expression patterns. Sheet 4: Quantification of enrichment of cilia-related genes with the cilia-related expression pattern. Sheet 5: List of target and bait genes with the chemosensory-related expression pattern. (XLSX) S1 Movie. Shown are 6 time-lapse videos of the phasmid cilia of N2 wild type, che-11 (e1810) and bbs-8(nx77) worms expressing the indicated GFP::RAB-28 reporter (all under the control of the endogenous rab-28 gene promoter). For illustration purposes, all videos played at 1.75x normal speed (original videos taken at 4 frames per second and played here at 7 frames per second). BB; basal body, cil; ciliary axonemes. Scale bars; 2 μm. (MP4) S2 Movie. Time-lapse video from phasmid cilia of N2 wild type worms expressing GFP:: RAB-28(GTP) at a low level (transgene oqEx304). RAB-28(GTP) localises almost exclusively near or at the periciliary membrane, with additional signals in the ciliary axoneme that move in both directions in an IFT-like manner. To make clear the IFT-like movement of RAB-28 (Q95L) particles, the image was overexposed and the video played at 5x normal speed. Scale bar; 2 μm. (MP4)