Recruitment of transcription factor ETS1 to activated accessible regions promotes the transcriptional program of cilia genes

Abstract Defects in cilia genes, which are critical for cilia formation and function, can cause complicated ciliopathy syndromes involving multiple organs and tissues; however, the underlying regulatory mechanisms of the networks of cilia genes in ciliopathies remain enigmatic. Herein, we have uncovered the genome-wide redistribution of accessible chromatin regions and extensive alterations of expression of cilia genes during Ellis–van Creveld syndrome (EVC) ciliopathy pathogenesis. Mechanistically, the distinct EVC ciliopathy-activated accessible regions (CAAs) are shown to positively regulate robust changes in flanking cilia genes, which are a key requirement for cilia transcription in response to developmental signals. Moreover, a single transcription factor, ETS1, can be recruited to CAAs, leading to prominent chromatin accessibility reconstruction in EVC ciliopathy patients. In zebrafish, the collapse of CAAs driven by ets1 suppression subsequently causes defective cilia proteins, resulting in body curvature and pericardial oedema. Our results depict a dynamic landscape of chromatin accessibility in EVC ciliopathy patients, and uncover an insightful role for ETS1 in controlling the global transcriptional program of cilia genes by reprogramming the widespread chromatin state.


INTRODUCTION
Ciliopathies are a group of multiorgan disorders whose aetiologies lie with primary ciliary dysfunction, leading to an expanding spectrum of diverse developmental and degenerati v e diseases ( 1 , 2 ). Although individual ciliopathy diagnoses are rare, their total contribution to human disease burden is surprisingly equal to that of several more common diseases, with an estimated incidence of ∼1:1000 ( 3 ). To date, > 20 phenotypically distinguishable symptoms have been associated with ciliopathies. Although individual diseases are known for the most commonly affected organ, for each of these ciliopathies remar kab le phenotypic variability has been observed even between members of the same family, making clinical diagnosis a challenge. For example, Ellis-van Cre v eld syndrome (EVC) is one of the skeletal ciliopa thies tha t not only exhibits the general clinical symptoms of ciliopathies, such as pol ydactyl y, characteristic skeletal fea tures, congenital cardiopa thy, and de v elopmental and intellectual delays, but may also involve unique symptoms, such as dental anomalies and ectodermal dysplasia (dysplastic nails and hair abnormalities) ( 4 ). Emerging studies have identified numerous cilia genes involving primary cilia impairment that may explain EVC ciliopathy pa thogenesis, indica ting tha t ciliopa thies share common pathway defects in which the cilia act as a nexus ( 5 ).
Cilia, whose dysfunction results in ciliopathies, are microtubule-based structures located on the surface of most vertebrate cells. The assembly and disassembly of cilia rely on rigorous transcriptional programs that integrate cooperati v ely to direct cilia and tissue formation and de v elopment (6)(7)(8). Maintenance of cilia formation and function involv es e xtensi v e di v ersity in gene e xpression and dynamic interaction between the cilia regulation axis and multimeric protein complexes that impact the downstream Hedgehog (HH), Notch and Wingless-Int (WNT) signalling pathways (9)(10)(11)(12). Ther efor e, these cilia genes are crucial to human development, especiall y w hen considering their collecti v e contributions to the de v elopment of ciliopathies. To date, hundreds of cilia genes have been implicated in established ciliopathies; meanw hile, potentiall y hundreds more cilia genes associated with ciliary structures and / or functions could result in as yet unknown or novel ciliopathies ( 13 ). Studies based on candidate gene approaches have provided insightful views for understanding the pathogenesis of ciliopathies ( 14 , 15 ). Ne v ertheless, additional diseases hav e emerged that were not traditionally considered to be ciliopathies; therefor e, pr evious r esear ch focusing on individual genes is insufficient to fully explain the transcriptional alterations that happen in ciliopathy-critical e v ents. Conv ersely, v ery little information exists regarding the gene regulatory networ ks involv ed in ciliopathies. Ther efor e, to understand the molecular basis of ciliopathies, it is critical to identify the regulatory elements governing the changes in cilia gene expression.
Alterations in epigenetic landscapes are known to contribute to gene regulatory networks ( 16 ). As a substantial part of epigenetic regulation, the collaboration of chromatin accessibility sites with regulatory elements profoundly changes the chromatin state, which subsequently influences the downstream tr anscriptional progr am ( 17 ). Indeed, among the genomic factors that may affect biological behaviour at a gi v en locus, chromatin accessibility is the most important, and pioneer factors that target closed chromatin can lead to its opening (18)(19)(20). The use of the innovati v e epigenetic Assa y f or Transposase-Accessible Chromatin with sequencing (ATAC-seq) has advanced our understanding of the machinery of epigenetic modification and gene expr ession r egulation in processes such as tumorigenesis and embryonic de v elopment (21)(22)(23). Although previously identified genetic alterations have been recognized in ciliopathies, the heterogeneity in the specific genetic alterations between patients provides convincing evidence that epigenomic regulation is crucial to understanding ciliopathy pathogenesis. The disruption of the structure of topo-lo gicall y associa ted chroma tin domains (TADs) led to abnormal limb formation through the rewiring of epigenomemediated gene-enhancer interactions in a case of ciliarela ted ciliopa thy ( 24 ). Epigenetic mechanisms, including DNA methylation and histone modifications, were found to be vital for cooperation between cilia genes and the regulatory network in polycystic kidney disease (PKD)-associated ciliopathies ( 25 , 26 ). The widespread phenotypic changes in ciliopa thies indica te tha t e xtensi v e alterations in the gene expression programs that control cilia disassembly and signalling response must be required for the disease to de v elop ( 27 ). Ther efor e, it is imperati v e to understand the specific regulatory changes controlling gene expression that lead to ciliopa thy pa thogenesis.
In this study, we explored the global redistribution of chromatin accessibility leading to broad transcriptome effects in EVC ciliopathy patients. We observed a significant increase in accessibility sites, defined as EVC ciliopathyactivated accessibility regions (CAAs), and modulation of the kinetics of neighbouring cilia genes. Furthermore, the transcription factor (TF) ETS1 was recruited to the CAAs and acted as a fundamental regulator for global alterations of the chromatin state. Suppression of ETS1 reduced the activity of CAAs and disturbed the expression program of cilia genes. Consistent with the phenotype of ciliopathy patients, disruption of ets1 expression in zebrafish caused bod y curva ture and heart d ysfunction. Collecti v ely, our findings re v eal the role of ETS1-dri v en chromatin accessibility redistribution in the regulatory network modulation of critical cilia genes, indicating a novel appendagepa tterning pa thway pr eviously unr ecognized in EVC ciliopa thy pa thogenesis.

Human material and cell culture
Venous blood was collected using standard procedures under ethical approval of the Institutional Ethics Committee of the Central People's Hospital of Zhanjiang (KY-YS-2021-05), with informed consent obtained from the patients and healthy donors. The peripheral blood mononuclear cells (PBMCs) were isolated using Ficoll (Cytiva, 17544202, USA) according to the manufacturer's instructions.

Lentivirus pr epar ation and tr ansfection
The second-generation packaging system plasmids pMD2.G and psPAX2 (Addgene, USA) were used to create lentivirus in 293T cells. Transfection was performed on 5 × 10 6 293T cells growing in a 10 cm culture dish using polyethylenimine (PEI) with 4 g of lentiviral vector, 3 g of psPAX2 and 1 g of pMD2.G (lentiviral vector:psPAX2:pMD2.G = 4:3:1). Cells were collected e v ery 24 h between 24 h and 72 h post-tr ansfection, ultr acentrifugation was used to concentrate the supernatants, and the virus titre was determined through successi v e dilutions.

Antibodies
All antibodies are listed in the Supplementary Data.

A T A C-seq and data processing
Fresh isolated PBMCs and hTERT RPE-1 cells were harvested and resuspended with cold lysis buffer for 10 min on ice. Next, 50 000 nuclei were counted using trypan blue staining and pr epar ed for transposition r eactions according to the manufacturer's instructions for the TruePrep DNA Library Prep Kit V2 for Illumina (Vazyme, TD501, China). Subsequently, nuclei were pelleted and resuspended with transposase at 37 • C for 30 min. DNA fragments were then purified with AMPure XP beads (Beckman, A63882, USA). DNA libraries were constructed after eight cycles of polymerase chain reaction (PCR) amplification using the TruePrep DNA Library Prep Kit and sequenced on the Novaseq PE150 platform to a depth of 4.0 × 10 7 reads.
For ATAC-seq data processing, FastQC (version 0.11.7) was used to evaluate the quality of the sequencing data, and Trimmomatic was explored to remove the adaptor sequences and obtain the clean data, which were subsequently aligned to the hg38 r efer ence genome using Burrows-Wheeler alignment (BWA). Multipl y ma pped reads were removed using SAMtools (version 1.3.1) and the MASCS2 was applied to call significant peaks with a q-value < 0.05. To identify differ entially expr essed genes between healthy donors and patients, the ATAC-seq peaks of each sample were merged to generate a consensus set of unique peaks. The number of peaks among this set was counted using bedtools (version 2.25.0), and CAAs and ciliopa thy-inactiva ted accessibility regions (CIAs) were identified using DESeq2 (version 1.16.0) ( 28 ), with the thresholds |log2FC| > 1 and P < 0.05. Genomic features of peaks were annotated using the ChIPseeker R package. Additionally, Homer software (version 4.6) and CentriMo (version 5.4.1) were applied to discover binding motifs in CAAs and CIAs. Genomewide normalized signal coverage tracks wer e cr eated by bam Coverage in deepTools (version 3.3.0) and were visualized in the Integrati v e Genomics Vie wer (IGV v ersion 2.5.0).

CUT&Tag and data processing
First, 10 000 fresh PBMCs and hTERT RPE-1 cells were counted and processed according to instructions for the Hyperacti v e ® Uni v ersal CUT&Tag Assay Kit for Illumina (Vazyme, TD903, China). Briefly, ConA beads (concanavalin A-conjugated paramagnetic beads) were added to the cell pellet and incuba ted a t room tempera ture for 10 min. Subsequently, the mix ed complex es wer e incuba ted overnight a t 4 • C with primary antibod y. The cells were then washed, followed by incubation with diluted secondary antibody at room temperature for 2 h. Antibodyprotein-DNA complex es wer e then linked by pA / G-Tnp and incubated at room temperature for 1 h. After washing, the eluted complexes were subsequently sheared by adding 5 × Trueprep Tagment Buffer L (5 × TTBL) provided in assay kit. Proteinase K and DNA extract beads were used for obtaining eluted DNA. Sequencing libraries were established after 11 cycles of PCR amplification. Selected sizes were processed using AMPure XP beads and the DNAs were sequenced on the Novaseq PE150 platform to a depth of 2.5 × 10 7 reads. For CUT&Tag data processing, FASTQC was used to measure the data quality distribution, and the clean reads were then aligned to reference genome sequences using the bwa program. Peaks were called using MACS2 peak caller with q-value < 0.05. Genome-wide normalized signal coverage tracks were visualized in the IGV. Peaks that were > 2-fold changed at a false discovery rate (FDR) < 0.1 were considered differential peaks, while peaks that did not have a log2 fold change (FC) significantly different from zero were termed constituti v e peaks.

ChIP-qPCR analyses
PBMCs were cross-linked with 1% formaldehyde solution (Sigma, F8775, USA) for 10 min at room temperature. Glycine was added for 5 min to quench the fixation reaction. Subsequently, the chromatin was sheared into 300-500 bp fragments using Bioruptor ® . The sheared DNA was pr e-clear ed by adding 20 l of salmon sperm DNA / Protein A Sepharose beads (Sigma, GE17-5280, USA) and incuba ted a t 4 • C for 30 min. Next, 3 g of primary antibody was incubated with a 250 g fraction of beads-DNA mixture and incubated overnight at 4 • C. The following day, antibody-protein-DNA complex es wer e pulled down using 30 l of Protein A Sepharose bead slurry. After washing, the eluted complexes were subsequently de-cross-linked at 65 • C for 5 h. Finally, the ChIP (chromatin immunoprecipitation) DNA was purified using phenol-chloroform extraction for further quantitati v e PCR (qPCR) study. The primers used in ChIP-qPCR are listed in the Supplementary Data.

Immunofluorescence and microscopy
In brief, PBMCs and hTERT RPE-1 cells for immunofluor escence wer e seeded on coverslips in 100 mm culture dishes. At 70-80% confluence, cells were washed twice with 1 × icecold phospha te-buf fer ed saline (PBS) and fix ed with 4% paraformaldehyde at room temperature for 20 min. Subsequently, the cells were blocked with 1% bovine serum albumin (BSA) at room temperature for 30 min and incubated overnight at 4 • C diluted in 1% BSA with primary antibody. Subsequently, the cells were washed and incubated with 150 l of secondary antibod y a t room temperature for 2 h. 4 ,6-Diamidino-2-phenylindole (DAPI) stain was added with a final concentration of 2 ng / l, and staining time was limited to 3 min. Finally, slides were sealed with 5 l Fluoromount-G ® mounting medium, and observed under a fluorescence microscope after 5 h.
All the centrosomal and cilia proteins were observed at room temperature using a TH4-200 fluorescence microscope (Ol ympus, Ja pan) equipped with a ×60, 1.42 numerical aperture (NA) Apo oil objecti v e lens (Olympus, Japan) or an A1-SHR LFOV confocal microscope (Nikon, Japan) equipped with a ×60 , 1.42 NA and ×100 , 1.4 NA objec-ti v e lens (Nikon). Images were acquired using DP Controller software (Ol ympus, Ja pan) or NIS-Elements software (Nikon, Japan). For images observed using z-stacking, sections were acquired with 0.5 m distance between zsteps. All images were reconstructed to maximum projections. Images were processed using ImageJ and Photoshop (Adobe, USA).

Cilia measurement analyses
The Pythagor ean theor em (PyT) was used to achie v e the optimal type of measurement of cilia anatomy / structure in 3D space. Knowing the length of two sides in a right-angle triangle, the third side can be calculated using the PyT formula: a 2 + b 2 = c 2 . To calculate the precise cilium length c , the formula can be rewritten as: c = ( a 2 + b 2 ) 1 2 , where the length of the cilium on the maximum intensity projection is a , and the thickness of z-slices is b ( 29 ).

Expression analysis
For western blotting, the total protein was isolated using TRIzol reagent (Invitrogen, 10296028, USA) according to the manufacturer's protocol. The proteins were dissolved in 1% sodium dodecylsulphate (SDS), separated using SDS-polyacrylamide gel electrophoresis (PAGE), and transferred to a nitrocellulose membrane. The membranes wer e pr e-stained using Ponceau S buffer and blocked in 5% BSA for 1 h. The tailored membranes were incubated overnight at 4 • C with the primary antibodies. After washing, the membranes were incubated with a diluted secondary antibody for at least 2 h at room temperature, washed three times with PBST, and finally visualized using the Odyssey infrared imaging system (Odyssey, USA).
For total RN A anal ysis and data processing, the total RNA was extracted using TRIzol reagent according to standard protocols, followed by analysis using an Agilent 2100 BioAnalyzer (Agilent Technologies, USA) to verify RNA integrity. Library construction and high-throughput RNAseq were performed with the Novaseq PE150 platform to a depth of 1.5 × 10 7 reads. The sequencing adaptors were removed from the raw RNA-seq reads, and the obtained clean data were aligned to the human reference genome (hg38) using HISAT2 software (version 2.1.0.). For gene expression quantification, the sequencing reads within each gene were counted, and the counts were normalized using edgeR ( 30 ). Enriched Gene Ontology (GO) terms were identified with MGI Gene Ontology Term Finder using |log2FC| > 1 and P -value < 0.05 as thresholds.

RNA interference (RNAi)
For the knockdown experiment, two verified small interfering RN As (siRN As) of ETS1 ( 31 ) from GenePharma were transfected into hTERT RPE-1 cells. The cells were analysed 48-72 h after siRNA transfection. To obtain siRNA-r esistant ETS1 (r esETS1), silent mutations wer e incorporated into the siRNA-targeted ETS1 sequence. All sequences are provided in the Supplementary Data.

Zebrafish maintenance and microinjections
Zebrafish experiments were performed under ethical approval of the Institutional Ethics Committee of the Central People's Hospital of Zhanjiang (ZJDY2022-13). Adult zebrafish stocks of the AB strain were maintained under standard aquaculture conditions. All embryos were incubated at 28.5 • C.
Morpholinos were designed and synthesized by Gene Tools, LLC (Philomath, USA), including a standard control morpholino and a morpholino targeting the start codon region (AUG) of ets1 transcripts to block translation. Morpholino oligonucleotides of 2 ng / embryo of control or ets1 were separately injected into the yolk at the one-cell stage. Morpholino sequences are provided in the Supplementary Data.
The ets1 coding sequence (CDS) was ligated to the PCS2 + vector at the BamH ι site through Gibson assembl y. ets1 mRN A was in vitro transcribed using the mMES-SAGE mMACHINE SP6 kit (Invitrogen, AM1340, USA). For ov ere xpression and rescue, 200 pg of ets1 mRNA was injected into each embryo at the one-cell stage.

Whole-mount in situ hybridization (WISH)
The cmlc2 RNA probe was in vitro transcribed and digoxigenin labelled using the DIG RNA Labeling Kit (Roche, 11277073910). WISH was performed as reported ( 32 ).

Whole-mount immunofluorescence staining for zebrafish embryo
Embryos were fixed at the 8-somite stage using 4% paraformaldehyde at 4 • C overnight. Anti-acetyl-alphatubulin antibody (1:500, Sigma, MABT868, USA) was used to label the cilia in Kupffer's vesicles (KVs). Immunofluorescence staining was performed as previously described ( 33 ). Before permeation with acetone, embryos were equilibrated in Tris-HCl (0.1 M, pH 9.0) for 5 min at room temperature and then heated at 70 • C for 15 min for antigen retrieval. After staining, embryos were flattened and imaged using the z-stack function of the confocal microscope with a ×60 oil objecti v e.

Behavioural assays
The behavioural assay was performed using a ZebraBox system and video tracking software (ViewPoint Life Sciences, France). At 5 days post-fertilization (dpf), an individual larva was placed in each well of a 96-well plate containing 500 l of embryo medium without light exposure. Plates were sealed with an optical adhesi v e film to pre v ent e vaporation. The mov ement of each larva was recor ded ov er a period of 20 min. A transparent background mode with a detection threshold of 20 was set. Behavioural endpoints were swimming distance (cm) and time (s).

Heart failure assessment
The zebrafish larvae were subjected to video recording under a Zebralab Blood Flow System (ViewPoint Life Sciences, France). To determine differential atrial and ventricular heart rate, heart videos were analysed using Mi-croZebraLab (ViewPoint Life Sciences). The pixel density changes detected by the software correlate with cardiac muscle contractions and chamber filling, which are used to estimate cardiac muscle contractions in beats per minute (bpm).
Blood flow videos were analysed using ZebraBlood (ViewPoint Life Sciences). The pixel density changes detected by the softwar e ar e combined with vessel diameter to determine the flow rate for each frame. Quantitati v e assessments were performed using video-based analysis, from which the linear v elocity, b lood flow and vessel diameter were evaluated.

Statistical analysis
Either a two-tailed unpaired Student's t -test or a oneway analysis of variance (ANOVA) with Dunnett's multiple comparisons test was applied to determine the difference between groups, unless otherwise noted, using GraphPad Prism 5 (version 5.01; GraphPad, USA). The r esults ar e pr esented as the mean ± standard error of the mean (SEM), and P < 0.05 was considered to indicate a statistically significant difference, unless otherwise stated.

Open chromatin landscape undergoes distinct global remodelling in EVC ciliopathy patients
One of the common malformations associated with ciliopathies is post-axial pol ydactyl y. Her e, we r ecruited patients with pol ydactyl y predominantl y of the hands and / or of the feet; based on pr evious r eports and sequencing inf ormation, f our of these patients were clinically diagnosed with the EVC ciliopathy, due to the additional presence of chondrodysplasia (dental abnormalities, short rib), congenital cardiopathy, intellectual disability and variable degrees of ectodermal d ysplasia (d ysplastic nails, hair abnormalities) ( Figure 1 A; Supplementary Data). Since cilia defects ar e r esponsible for the phenotypes in ciliopathies, and since it has been reported that PBMCs can be ciliated ( 34 ), we first examined the cilia morphology of PBMCs isolated from the patients. We found that the cilia in the PBMCs of these patients were longer than those in normal individuals, demonstr ating the structur al cilia defects in the patients, which strengthens the correlation between cilia malfunction in patient-deri v ed cells and the symptoms observed in the patients (Figure 1 B). Additionally, the transcriptional e xpression le v els of se v eral well-recognized cilia genes in the PBMCs, including CEP131 , NEK8 and other cilia genes recorded in the SCGSv2 gold standard list, were substantially elevated in the genome-wide transcriptome profile of EVC ciliopathy patients ( 35 ) (Supplementary Figure S1A-C). Notably, the differentially expressed genes were always enriched in multiple signalling processes --such as body morphogenesis, metabolic processes and immune responses --that have been shown to be closely associated with cilia genes ( 36-39 ) (Supplementary Figure S1D, E). In particular, the genes ov ere xpressed in patients were always enriched in cilia-rela ted pa thwa ys, reinf orcing that a common set of cilia genes were activated during EVC ciliopathy pathogenesis (Figure 1 C).
Chromatin accessibility is a reflection of and corresponds to the transcriptome profile that determines the cellular state ( 40 ). Gi v en the unique alterations in cilia genes and the enrichment of cilia-rela ted pa thways, we conducted ATACseq analyses in isolated PBMCs from both patients and healthy donors (Supplementary Figure S2A). As expected, the ATAC-seq data showed clear nucleosome phasing in the insert size distributions, and the r eads wer e a ppropriatel y enriched at transcription start sites (TSSs) (Supplementary Figure S2B, C). Peaks r epr esenting the chromatin accessibility r egions wer e primarily located in the intergenic regions containing abundant cis -regulatory elements; furthermore, genome annotation re v ealed that a large fraction of these ATAC-seq peaks overlapped with promoter and enhancer regions, consistent with previous findings ( 41 )

Identification of specific chromatin accessible regions in EVC ciliopathy patients
To identify chromatin remodelling that potentially contributes to EVC ciliopathy pathogenesis, we defined > 10 000 differentially accessible regions (DARs) between EVC ciliopa thy pa tients and healthy donors based on quantitati v e peak signals; these DARs included 402 EVC CAAs and 568 EVC CIAs (Figure 2 A; Supplementary Figure  S4A, B). As expected, promoter elements were enriched in both CAAs and CIAs according to chroma tin sta te analysis (Figure 2 B). Furthermore, most CAAs and CIAs were located ∼2 kb from TSSs, indicating the mapping of the potential promoter elements in the peaks mentioned above (Supplementary Figure S4C-E). Using ChIP-qPCR of H3K4me3, randomly selected peaks in CAAs and CIAs were validated based on their promoter signal alterations (Supplementary Figure 4F). By comparing our data with public DNase-seq data from ENCODE, we found that > 50% of the peaks in CIAs overlapped with the chromatin open regions in data from healthy donors, whereas the fractions of overlapped peaks in CAAs w ere substantially low er in all healthy donors, r epr esent-

CAAs regulate adjacent cilia gene expression
Because promoters have the ability to act as enhancers to regulate the expression of nearby and distal genes in cell fate decision ( 26 ), we investigated the impact of EVC ciliopathy-specific DARs on global gene expression. Interestingly, CAA-neighbouring genes were up-regulated, while those adjacent to CIAs were down-regulated (Figure 2 D). Additionally, genes known to be up-regulated in patients, such as NEK8 , showed substantial elevations in both chromatin accessibility and e xpression le v el (Figure 2 E). Verification of the H3K4me3 signal near NEK8 underscored the pre vious e vidence for increased promoter acti vity at CAAs (Figure 2 F).
Next, GO analysis was used to identify the biological functions of DAR-adjacent genes. Notably, CAA-adjacent r egions wer e enriched in morphogenesis-associa ted pa thwa ys, including embry onic digit morphogenesis, cell morphogenesis and morphogenesis of other organs, which are crucial for the de v elopment process (Figure 3 A). Meanwhile, genes neighbouring CIAs wer e pr edominantly enriched in metabolism-associated pathways such as protein ubiquitina tion and kera tiniza tion (Figure 3 B). Consistent with the GO analysis results, CAA-adjacent genes were significantly up-regulated, including the well-known cilia genes CEP131 , CEP41 and CDC20 (Figure 3 C, D; Supplementary Figure S5A, B). Conversely, RNF216 , the key regulator of protein ubiquitination in de v elopment, was flanked by CIAs and showed decr eased expr ession in patients (Figure 3 C; Supplementary Figure S5C). Since w e show ed that the genes neighbouring the CAAs always displayed upr egulated expr ession, and that many of these aberrant elevated CAA-adjacent genes were well-known cilia genes, we wondered whether these cilia genes had a regulatory role in chroma tin accessibility d ynamics. Ther efor e, we quantified the ATAC-seq intensities from the gold standard SCGSv2 list, and found that peaks in patients led to substantial enrichment of cilia genes compared with the signal in healthy donors, while the randomly selected gene set displayed unchanged signals in patients and healthy donors (Figure 3 E).
Collecti v ely, our results suggest that the CAA-harbouring EVC ciliopathy accompanied by multiorgan de v elopmental abnormalities is due to the impaired cellular signalling; moreov er, the ne wly acti vated CAAs are tightly coupled to a robust increase in the expression of adjacent cilia genes, which may serve as a pr er equisite or contributing factor for the aberrant elevated cilia gene expression seen in patients. The relati v e le v els of H3K4me3 for CAAs in the NEK8 gene measured using ChIP-qPCR in PBMCs (lower). H3K4me3 enrichment for the P53 intron served as a negative control, and H3K4me3 enrichment for P53 and RAD51 served as positive controls, as previously described. The enrichment is normalized to 10% input, while IgG is used as a negati v e contr ol. Err or bars r epr esent the means ± SEM for thr ee independent experiments. * P < 0.05, ** P < 0.01, *** P < 0.001, as determined using one-way ANOVA with Dunnett's multiple comparisons test.

Recruitment of ETS1 to CAAs contributes to regulation of the expression of cilia genes
To identify the potential TFs involved in chromatin remodelling at CAAs, thus contributing to the regulation of downstream cilia genes, we carried out motif enrichment analyses using the algorithms HOMER and MEME. The topranked motifs were observed in the consensus binding sites for ETS, bHLH, ZF and the STAT famil y, but onl y the ETS motif was strongly over-r epr esented in CAAs ( Figure  4 A; Supplementary Figure S6A, B). Ther efor e, we carried out footprinting analyses and found considerable TF occupancy around the ETS motif site in EVC ciliopathy patients relati v e to healthy donors, suggesting a protracted ETS-DNA interaction on CAAs during EVC ciliopathy pathogenesis (Figure 4 B). The ETS family of proteins comprises 28 TFs that contain a highly conserved DNA-binding ETS domain ( 44 ) (Supplementary Figure S6C). To identify putati v e TFs that regulate CAAs, we subsequently screened the ETS TFs exhibiting expression changes and found that ETS1 expression --both for mRNA and for protein --was notably higher in EVC ciliopathy patients relati v e to healthy donors ( Figure 4 C; Supplementary Figure S6D, E). To further determine whether ETS1 maintains the chromatin landscapes, we assessed ETS1 occupancy in PBMCs from both EVC ciliopathy patients and healthy donors by using CUT&Tag. Notably, ETS1 displayed a genome-wide ( D ) Distribution of motif scores of sites within ETS1 CUT&Tag peaks in PBMCs, either those that gave more signal in CPY than in the Ctrl (differential) or those that were not significantly different (constituti v e). The maximum scoring ETS full site within each CUT&Tag peak was used. The constituti v e peaks have a higher mean motif score than the differential peaks (Mann-Whitney U-test, P < 1 × 10 −32 ). ( E ) Spearman's correlation of differential ETS1 CUT&Tag signals in accessible regions and differential accessibility in CPY versus Ctrl ( r = 0.399). ( F ) The enrichment of ETS1 CUT&Tag signals at CAAs in both Ctrl and CPY. Tracks ar e centr ed at the peaks and extend ± 3 kb. ( G ) Enrichment plots showing normalized read densities of ETS1 CUT&Tag signals at the TSS for randomly selected genes ( n = 700, left) and SCGSv2 genes (right) in both Ctrl and CPY. Tracks ar e centr ed at the TSS, extending ± 3 kb. ( H ) The relati v e le v els of cilia genes flanking CAAs in Ctrl and CPY measured using ETS1 ChIP-qPCR in PBMCs. ETS1 enrichment for RPL3 -c (coding region) served as the negative control, and ETS1 enrichment for RPL3 -p (promoter region) and RPL13A -p (promoter region) served as positive controls, as previously described. The enrichment of ETS1 is normalized to 10% input. ** P < 0.01, *** P < 0.001, as determined using one-way ANOVA with Dunnett's multiple comparisons test.  Figure S6F). We quantified and compared the motif scores of the CUT&Tag peaks of controls and EVC ciliopathy patients, and found that, on average, peaks that were differentially occupied had lower motif scores than peaks that were constituti v el y occupied, w hich is consistent with the finding that increased ETS1 expression enables binding of lower affinity sites (Figure 4 D). Mor eover, the differ ential accessibility and CUT&Tag signal wer e r elati v ely correla ted, indica ting tha t a high fraction of differential accessibility is attributable to differential ETS1 occupancy (Figure 4 E). In particular, the ETS1 enrichment signal on CAAs was elevated in EVC ciliopathy pa tients rela ti v e to that in healthy individuals, r einfor cing the potential function of ETS1 in modulating the activity of CAAs (Figure 4 F).
Since we have revealed that CAAs are responsible for robust downstream cilia gene expression, ETS1 enrichment among randomly selected genes was compared with the signal around cilia genes in the SCGSv2 gold standard list. Indeed, these cilia genes exhibited pronounced ETS1 enrichment flanking their TSSs relati v e to a randomly selected gene set (Figure 4 G). Moreover, ChIP-qPCR of ETS1 was performed to examine its signal on CAAs flanking cilia genes, which showed that the signal strength of ETS1 was substantially increased on these CAAs (Figure 4 H). IGV mapping of ETS1 occupancy on CEP131 and NEK8 also emphasized the role of ETS1 in recruitment to this subset of CAAs (Figure 4 I). Altogether, these results suggest that the recruitment of ETS1 to CAAs substantially altered the chroma tin sta te transitions, ultima tely causing downstream tr anscriptional alter a tions associa ted with incr eased expr ession of cilia genes.

ETS1 is r equir ed f or cilia f ormation
Previous studies have hypothesized the potential function of genes from the ETS family in regulating the expression of cilia genes during tumour progression ( 45 ). As a member of the ETS family, ETS1 is suspected to be an oncogene involved in various cancers and to have a vital function in haematopoietic cell de v elopment ( 46 , 47 ). To further inv estigate the potential role of ETS1 in EVC ciliopathy, we performed ETS1 ov ere xpression and knockdown experiments in hTERT RPE-1 cells and observed the impact on cilia formation after serum starvation. It was demonstrated that for ced ETS1 expr ession r esulted in se v ere cilia morphology impairments including bulge, truncation ( < 3 m) and elongation ( > 8 m), with a higher percentage of elongated than truncated cells ( Figure 5 A-D), although the increase in the overall percentage of ciliated cells was insignificant (Figure 5 A, E). Conversely, it was found that ETS1 knockdown led to almost the opposite effects of its ov ere xpression. ETS1 knockdown impaired cilia morphology with remar kab ly increased bulged cilia, a higher percentage of cells with truncated cilia than elongated cilia and a slight increase in the overall percentage of cells that underwent cilia formation. These morphology impairments, along with the quantity of ciliated cells, were rescued by treatment with exogenous siRNA-resistant ETS1 ( Figure 5 F-J; Supplementary Figure S6G). Collecti v ely, these r esults r e v ealed the criti-cal and unexpectedly bidirectional regulator role of ETS1 in cilia formation.

Disrupted chromatin accessibility caused by ETS1 silencing underlies the impaired transcriptional program of cilia genes
To determine if prolonged ETS1 expression is r equir ed to sustain an open chroma tin sta te of CAAs, we performed a stable ETS1 knockdown in hTERT RPE-1 cells; using ATAC-seq on siETS1 samples and negati v e controls, we found that ETS1 silencing caused e xtensi v e alterations in chromatin accessibility (Supplementary Figure S7A). The peak annotation suggested that the closure of accessibility in ETS1 knockdown cells occurred primarily in intergenic regions enriched in ETS-family motif sites (Supplementary Figure S7B, C). Notably, ETS1 knockdown substantially r educed accessibility, r e v ealing an ETS1-dri v en disruption of the global open chroma tin sta te, supporting the view that ETS1 is capable of inducing chromatin remodelling by evicting nucleosomes directly from nucleosomebound DNA ( 48 ) (Figure 6 A). In contrast to the elevated CAA peaks mentioned above, the suppr essed expr ession of ETS1 resulted in a robust decrease of CAA activity (Figure 6 B). Meanwhile, considering that nucleosomes and TFs are in competition for DNA-binding sites, we evaluated nucleosome occupancy after ETS1 knockdown by calculating the ATAC-seq data insert sizes. We observed an increase in nucleosome occupancy at ETS motif sites in ETS1silenced cells, leading to less ETS1 occupancy on CAAs (Supplementary Figure S7D). Additionally, the closed chromatin regions in ETS1 knockdown cells were generally anticorrelated with accessibility changes in CAAs, highlighting the function of ETS1 in restructuring chromatin accessibility during EVC ciliopathy pathogenesis (Figure 6 C).
Gi v en that the SCGSv2 genes exhibited substantially decreased densities after ETS1 silencing relati v e to randomly selected genes, we then tested ETS1 occupancy in hTERT RPE-1 cells by using CUT&Tag techniques (Figure 6 D; Supplementary Figure S7E). In ETS1 knockdown cells that alwa ys displa yed ETS1 occupancy, we observed the collapse of CAA activity in regions including the well-described cilia genes KDM3A and SMO , further supporting the role of ETS1 in harbouring CAAs that flank cilia genes ( 49 , 50 ) (Figure 6 E). To verify whether ETS1 alone can sufficiently r egulate the expr ession of cilia genes, we carried out transcriptome profiling in ETS1 knockdown cells and found e xtensi v ely decreased e xpression of cilia genes adjacent to CAAs corresponding to ETS1 silencing (Supplementary Figure S7F, G). Gene set enrichment analysis (GSEA) of the transcriptome profiling re v ealed that many cilia-related signalling pathways, including the HH pathway, were impaired after ETS1 depletion (Supplementary Figure S7H-J). This finding inspired us to investigate whether ETS1 is r equir ed for ciliary signalling transduction. GLI Family Zinc Finger 3 (Gli3), a critical TF in HH signalling, exists in two forms: a full-length activator (Gli3-FL) or a r epr essor (Gli3-R). Gli3-FL activates the HH pathway and modulates HH genes by targeting the Gli1 promoter, whereas Gli3-R --produced by degradation of Gli3-FL --inhibits HH functions. In this study, Gli3-R expression decreased and Gli1 expression substantially increased in Smoothened  (SMO) agonist (SAG)-treated ETS1-ov ere xpressing cells, indica ting tha t forced ETS1 e xpression acti vated the HH signalling pathway (Supplementary Figure S8A). Furthermore, depletion of ETS1 inhibited the HH signalling pathway, and exhibited the opposite response in Gli1 and Gli3 expr ession patterns compar ed with ETS1 over expr ession (Figure 6 F). These results demonstrated that ETS1 plays a critically positi v e role in the HH signalling pathway.
Consistently, SAG-tr eated ETS1-over expr essing cells exhibited activated HH signalling with a notable increase in Gli3 accumulation at the cilia tip ( Supplementary Figure S8B, C). Meanwhile, SAG-induced ciliary SMO localization was a pparentl y enhanced in ETS1-ov ere xpressing cells (Supplementary Figure S8D, E). Conversely, depletion of ETS1 inhibited HH signalling transduction with impaired Gli3 accumulation at the cilia tip and weakened SMO localization along the cilia (Figure 6 G, H; Supplementary Figure S8F, G). Ther efor e, these phenomena further confirmed that ETS1 is a positi v e regulator of the canonical HH signalling pa thway. Altogether, our da ta support a EVC ciliopathy model in which ETS1 is recruited to CAAs to stabilize the open chroma tin sta te, thereby inducing a genome-wide program of expression of cilia genes related to cilia-associated processes such as signalling transduction.

ETS1 is sufficient for cilia development in zebrafish larvae
ETS1-dri v en redistribution of transcriptome profiling in cilia genes in patients provides an important context for understanding organ de v elopment. Here, we e valuated these de v elopmental changes in zebrafish, a powerful model system for the study of cilia structure, and examined the changes in cilia formation and subsequent morphological and physiological effects caused by alterations in ETS1 expression. A morpholino targeting the 'AUG' site of the ets1 transcript was designed and synthesized (Supplementary Figure S9A). At 72 hours post-fertilization (hpf), ets1 morphants displayed various typical cilia defect phenotypes, including body curvature (29.1%) and pericardial oedema (82.6%). Intriguingly, ov ere xpression of ETS1 also led to these cilia defect phenotypes, but with lower frequency (bod y curva ture, 20.8%; pericardial oedema, 43.9%) (Figure 7 A, B). Meanwhile, we performed WISH using a cmlc2 probe to examine heart position distribution. In ets1 morphants, 15.2% and 3.1% of embryos displayed middle and right-sided heart, respecti v ely. In line with the morphogenetic phenotype r esults, ETS1-over expr essing embryos also displayed left-right asymmetry defects; howe v er, only 3.3% and 0.7% of the control embryos displayed middle and right-sided heart, respecti v ely (Figur e 7 C, D). These r esults imply that ets1 plays a vital role in determining left-right patterning. Notably, an incr eased per centage of truncated cilia compared with the control group could be detected not only in ETS1-deficent cells but also in ETS1-ov ere xpressing cells, although with a lesser effect. Thus, we observed cilia defect phenotypes --including body curvature, pericardial oedema and disorder of left-right asymmetry --in both ets1 morphants and ETS1-ov ere xpressing zebrafish. Moreov er, we note that both depletion and ov ere xpression of a protein called nuclear distribution gene C (NudC) tha t regula tes cil-iogenesis can lead to the above-mentioned cilia defect phenotypes in zebrafish ( 51 ), indica ting tha t both e xcessi v e and insufficient expression of cilia-related genes in zebrafish can result in the occurrence of cilia defect phenotypes, and that the pr ecise expr ession of cilia genes is critical for zebrafish de v elopment.
In zebrafish, left-right asymmetry patterning is regulated by the ciliated organ the KV ( 52 ). Hence, we also labelled the cilia in the KV at the 8-somite stage, measured the length of the cilia and discovered that those of ets1 morphants were significantly shorter than those of control embryos (Figure 7 E, F). To validate the specificity of phenotypes in ets1 morphants, we co-injected a mixture of in vitro transcribed ets1 mRNA with ets1 morpholino. Notably, all of the defects caused by ets1 morpholino, including body curvature, pericardial oedema, disruption of left-right asymmetry, and cilia length defects in the KV, were effecti v ely rescued by exogenous ets1 mRNA (Figure 7 B-F). These results rule out morpholino-induced toxicity or secondary effects, and validate the specificity of the phenotypes in ets1 morphants.
Since ETS1 has been shown to have an essential function in angiogenesis, we speculated that heart pumping capacity and blood flow would change after ETS1 r epr ession or ov ere xpression, which may lead to pericardial oedema in zebrafish ( 53 ). Hence, we measured the heart rate in zebrafish larvae over a 10 min period and found that the atrial and ventricular beat rates were robustly reduced in ets1 morphants relati v e to controls, whereas rescuing ets1 in the ets1 morphants partially r estor ed the heart rates toward the le v els of the controls (Figure 7 G, H). As expected, in larvae ov ere xpressing ets1 , heart rates were decreased relati v e to controls (Supplementary Figure S9B, C). These findings suggest that a specific le v el of ets1 gene expression is necessary for heart de v elopment. Car diac morphology is one of the important indicators to evaluate cardiac function, which affects systemic blood circulation. To further address the impact of ETS1 on blood circulation, we quantified the blood flow dynamics by evaluating linear velocity, blood flow and vessel diameter, and discover ed sever e vascular stenosis in ets1 morphants and larvae ov ere xpressing  Figure S9D). This finding explains the heart pumping defects induced by either ETS1 r epr ession or overexpression. Swimming is a highly complex behaviour that involves deeply integrated physiological processes, and organ defects are known to profoundly affect ex er cise ability. Ther efor e, we analysed zebrafish locomotor trajectories to further assess behavioural changes in response to alterations in ETS1 expression. During a 20 min test, we noticed that the ETS1-suppressed and ETS1-ov ere xpressing zebrafish had substantially decreased swimming distance and swimming time relati v e to control zebrafish, while rescue of ets1 expression with mRNA in zebrafish with ets1 knockdown resulted in increased swimming distance and swimming time relati v e to zebrafish with ETS1 repression (Figure 7 J, K; Supplementary Figure S9G, H). In summary, our results provide evidence that ETS1 functions as the important regulator in organ de v elopment, the malfunction of which negati v ely affects the fins and heart, in line with the characteristics of ciliopathies.

DISCUSSION
A myriad of cilia genes involved in cilia formation and function profoundly affect cilia signalling transduction, leading to changes in physiological and de v elopmental functions. Deciphering the cilia gene expression program within comple x networ ks and pathways is crucial for understanding the underlying mechanism of multiorgan-associated ciliopathies. Ne v ertheless, knowledge of the effects of cilia gene regulation on ciliopathies is curr ently limited; ther efore, it is vital to thoroughly characterize the chromatin state-dependent mechanism for the cilia regulatory network. In this study, we investigated the dramatic remodelling of the chromatin state in EVC ciliopathy patients. The distinct CAAs were activ ated b y a single TF, ETS1, leading to further up-regulation of the expression of flanking cilia genes. We also confirmed aberrant ETS1 activation in EVC ciliopa thy pa tients, while suppression of the ETS1-induced open chromatin state collapsed, leading to cilia defects in both cellular and organismal contexts. These data suggest the strong likelihood that cilia genes r equir e ETS1-induced chroma tin sta te altera tions as their mechanism for EVC ciliopa thy pa thogenesis; mor eover, these findings underscor e the powerful combination of genome-wide epigenetic approaches with individual molecular studies to uncover previously unknown variations in rare diseases.
ETS-family TFs have been well characterized in se v eral types of solid tumours. Although their role in ciliopathies has not been studied, a recent study attempted to establish the link between ETS-family TFs and modulation of cilia gene expression at the tumour-micr oenvir onment (TME) interface ( 45 ). Our findings resolve the question of the mechanism by which ETS1 dri v es the redistribution of chromatin accessibility, ther eby dir ecting the r egulatory network of critical cilia genes in EVC ciliopathy patients, coinciding with a previous study showing that ETS1 can maintain the function of cilia by binding to the IFT20 promoter to facilitate its expression ( 54 ). In our study, the ETS1-regulated ciliary genes can be primarily divided into two groups based on their effects on cilia: the ciliary positi v e factors that form the majority of ETS1-regulated genes, and include WDR34 , RSPH1 , CEP41 , NEK8 , CEP131 and FUZ , and the ciliary negati v e factors that include only a few genes responsible for ciliary disassembly, such as CDC20 and KDM3A . Therefor e, it is r easonab le that the ETS1-ov ere xpressing cells e xhibited more elongated cilia than truncated cilia, thus creating a net positi v e effect including the slightly increased cilia formation and activation of HH signalling seen in ETS1ov ere xpressing cells, with a different but not completely opposite effect on cilia in ETS1-deficient cells. Ther efor e, unlike most unidirectional regulation of cilia in single gene defects, ETS1 demonstrates bidirectional regulation in cilia structure and function. Notably and consistently, continuously high ETS1 expression in EVC ciliopathy patients will lead to the aberrant activation of cilia genes; howe v er, lack of ETS1 expression leads to failure to open accessible regions tha t regula te corr ect expr ession of cilia genes, implying the essential role of ETS1 for de v elopment.
Usuall y, pol ydactyl y and bone defects ar e r elated to defecti v e HH signalling. Notab ly, ectopic e xpression of Shh in the anterior limb bud also causes the formation of extra anterior digits ( 55 ). Previous studies have reported that ETS / ETV-family TFs regulate Shh expression at E11.5 during the maintenance / expansion phase of Shh expression, and these studies also indicate that ETS1 activates HH signalling ( 56 , 57 ). Mor eover, a r ecent stud y has found tha t ov ere xpressed ETV2 can displace the chromatin of limb enhancers that ectopically activate Shh signalling and induce pol ydactyl y ( 58 ). Thus, it is reasonable that EVC ciliopathy patients who have enhanced HH signalling induced by ETS1 ov ere xpression e xhibit the symptom of pol ydactyl y.
Additionally, ther e ar e now well-supported experimental distinctions between confirmed and suspected cilia / ciliopathy genes, primary and motile cilia genes, and cilia genes that encode proteins that are acti v e only within cilia and those that encode proteins that are acti v e both within and outside of cilia. In our study, for example, many of these screened cilia genes function within cilia, such as CEP131 and CEP41, while others function both within and outside of cilia, such as FUZ --which governs trafficking of intr aflagellar tr ansport (IFT) pr oteins fr om the cytoplasm to the basal body and then to the ciliary tip ( 59 ) --and KDM3A, which r egulates r ecruitment of IFT proteins into the cilia by modulating actin dynamics through actin cytoskeleton binding ( 49 ). Collecti v ely, these results indicate the comprehensi v e modulation of the transcriptional program of cilia genes by ETS1.
Apart from ETS1, it is well demonstrated that RFX and FOXJ1 TFs directly regulate genes for core ciliary components. RFX TFs play essential roles in the generation of both motile and sensory cilia. The RFX proteins activate core components such as IFT122, IFT172, Dync2li1, Dnah9 and Dnah11, thereby regulating a series of ciliary genes that positi v ely regulate cilia formation and function. Meanw hile, FOXJ1 pro grams motile cilia by activating a network of motile cilia genes that promote cilia formation and function, for different cilia types, in selected cell types and or ganisms. Ho we v er, the TFs of the For khead (FOX) and RFX families either could not be recruited by accessible chromatin regions or were not activated in our TF-CAAs-cilia genes axis (Supplementary Figure S10A-C). In contrast to RFX proteins and FOXJ1, ETS1 activates two groups of cilia genes that can positi v ely and negati v ely modulate cilia formation and function, leading to bidirectional phenotypes very different from those generated by the ciliogenesis 'master regulator' TFs, including RFX proteins and FOXJ1. Notably, manipulation of ETS1 expression in different cell types may produce distinct cilia phenotypes, paralleling species-specific differences in the ability of TFs such as FOXJ1 to regulate cilia.
Recent r esear ch shows that ETS-family TFs can serve as transcriptional activators and / or r epr essors according to gene and conte xt, re v ealing the complicated multidirectional function of ETS-family TFs in cilia regulation ( 45 ).
Whether and to what extent the other ETS TFs function in cilia modulation, and how these family members precisely orchestrate cilia formation and function r equir es further elucidation. In addition to the ETS family, the TFs of other families may also participate in the control of expression of cilia genes, which will probably further valida te the grea t importance of chroma tin d ynamics and epigenetics in ciliopathies ( 14 , 60 ).
Based on previous findings, the ETS family members are oncogenic TFs connected to wide-ranging cellular processes such as self-renewal, DNA damage, metabolism, TME modula tion, d ynamic chroma tin remodelling and epigenetics ( 61 ); furthermore, cilia also take part in these processes. Recent studies have strengthened the genetic and functional links between the DNA damage response (DDR) and ciliary factors including CEP63, CEP152, CEP164, CEP290, MCPH1, NEK8 and PCNT ( 62 ). Additionally, primary cilium, as a spatially localized platform for various signal transduction pathwa ys, pla ys a crucial role in the TME, since paracellular signalling between tumour cells and other cells in the TME strongly af fects initia tion, progression and therapeutic efficacy in tumours ( 63 ). The signal transduction, mediated by primary cilium, also participates in energy and steroid metabolism ( 64 ), which affects tumour cell survival. Recently, the function of ETS-family TFs in modulating cilia gene expression at the TME interface has been partially elucidated ( 45 ). We have demonstrated that ETS1 r egulates the expr ession of cilia genes via dynamic chromatin remodelling epigenetics. Combining our results with the pr eviously r eported close r elationship among the ETS family, cilia and tumour formation, it is plausible that the modulation of ETS1 oncogenic efficacy is accomplished, at least in part, by primary cilia, which provides new insight into the mechanism for the oncogenic efficacy of ETS TFs, which needs further verification.
Moreover, the signalling pathway via primary cilium may also in turn affect the activity of the ETS TFs, such as FLI1. Meanwhile, Wnt / ␤-catenin signalling, primarily mediated by primary cilia, antagonizes the EWS-FLI1-dependent repression of transforming growth factor (TGF)-␤ receptor type 2 in Ewing sarcoma cells ( 65 ). Ther efor e, we speculate that there is a feedback loop between the ETS family and primary cilia and that this feedback loop is the means by which cells accomplish the pr ecise r egulation of cilia structure to ensure cilia function.

DA T A A V AILABILITY
All related sequencing data have been uploaded to NCBI's Gene Expression Omnibus and are accessible through accession number GSE209662.