Ciliary Genes Are Down-Regulated in Bronchial Tissue of Primary Ciliary Dyskinesia Patients

Primary ciliary dyskinesia (PCD) is a rare, genetically heterogeneous disease characterized by recurrent respiratory tract infections, sinusitis, bronchiectasis and male infertility. The pulmonary phenotype in PCD is caused by the impaired motility of cilia in the respiratory epithelium, due to ultrastructural defects of these organelles. We hypothesized that defects of multi-protein ciliary complexes should be reflected by gene expression changes in the respiratory epithelium. We have previously found that large group of genes functionally related to cilia share highly correlated expression pattern in PCD bronchial tissue. Here we performed an explorative analysis of differential gene expression in the bronchial tissue from six PCD patients and nine non-PCD controls, using Illumina HumanRef-12 Whole Genome BeadChips. We observed 1323 genes with at least 2-fold difference in the mean expression level between the two groups (t-test p-value <0.05). Annotation analysis showed that the genes down-regulated in PCD biopsies (602) were significantly enriched for terms related to cilia, whereas the up-regulated genes (721) were significantly enriched for terms related to cell cycle and mitosis. We assembled a list of human genes predicted to encode ciliary proteins, components of outer dynein arms, inner dynein arms, radial spokes, and intraflagellar transport proteins. A significant down-regulation of the expression of genes from all the four groups was observed in PCD, compared to non-PCD biopsies. Our data suggest that a coordinated down-regulation of the ciliome genes plays an important role in the molecular pathomechanism of PCD.


Introduction
Cilia are small cellular projectiles, extending from the cell surface, with which they share the cell membrane. The ciliary axoneme is built on a scaffold of nine peripheral microtubule doublets, associated with many multi-protein complexes, including outer and inner dynein arms (ODA and IDA), nexin links and radial spokes. Cilia act either as sensors (primary cilia) or propel fluid over the epithelia of various organs (motile cilia) [1]; their dysfunction is the underlying cause of many systemic diseases.
Primary ciliary dyskinesia (PCD) is a rare genetic disease characterized by recurrent respiratory tract infections, bronchiectasis and infertility. Pulmonary symptoms occur because of the lack of an efficient mucociliary clearance, caused by the kinetic dysfunction of motile cilia in the respiratory epithelium. Male infertility is caused by the dysmotility of flagella in spermatozoids. Situs inversus, a mirror reversal of body organ positioning, is present in approximately half of the PCD patients because of the immotility of primary cilia of the embryonic node [2,3]. In rare cases, PCD symptoms are a part of syndromic diseases, e.g. Xlinked retinitis pigmentosa [4,5] and oral-facial-digital type I syndrome [6]. The diagnosis of PCD is based on the clinical criteria, the delayed saccharin test [7], and low levels of nasal nitric oxide [8]. These symptoms are supported by the imaging studies of the respiratory cilia. Tissue specimens obtained through bronchoscopy are evaluated in the light microscope for ciliary beating and with the electron microscope for defects of the ciliary ultrastructure. Among a wide spectrum of various ultrastructural ciliary lesions documented in PCD patients [9], the lack of ODA and/or IDA is the most frequently observed defect.
Previously, in search for a better understanding of the molecular basis of functional defects of cilia in PCD, we have performed the whole-genome expression profiling in bronchial biopsies from six PCD patients [36]. Clustering analysis revealed a large group of genes with highly correlated expression pattern in PCD biopsies, but not in controls; based on the series of in silico analyses, we have indicated over 200 new genes potentially involved in the biology of human cilia. In the present study, we further explored the gene expression profiling data to characterize patterns of differential gene expression in bronchial tissue from PCD patients and non-PCD controls. We report that the significant proportion of genes that are down-regulated in PCD encode specific elements of the cilium, suggesting that a coordinated down-regulation of the ciliome genes plays an important role in the molecular pathomechanism of PCD.

Ethics statement
The Institutional Review Board (IRB) of the Medical University of Poznań approved use of bronchial biopsy specimens of patients and controls for genetic studies on PCD. Specimens were collected during routine hospital procedures. Written informed consent was obtained from all subjects.

Bronchial biopsies and subjects
Material for this study, obtained from six unrelated PCD patients and nine unrelated control individuals, has been described in detail in Geremek et al. [36]. Briefly, clinical evaluation of the PCD patients was performed at the Institute of Tuberculosis and Lung Diseases in Rabka, Poland, by the experienced physician (AP). The primary bronchopulmonary symptoms in the patients were: sinusitis, nasal polyps, bronchiectasis, and recurrent infections of the upper respiratory tract. PCD status was also indicated by the positive results of routine diagnostic tests: the delayed saccharine test and lack of the ciliary motility in light microscopy imaging. In four patients (#1, 2, 3, 4), PCD diagnosis was confirmed by the low level of NO (25-188 ppb), measured in the nasal cavity (chemiluminescence analyzer, the threshold value of 200 ppb) [8] (Table S1). Three of these patients (#1, 3, 4) had ODA/IDA defect. One patient (#6) for whom no NO data were available had IDA defect and another one (#5) in whom NO level was in the normal range had ODA/IDA defect. Two of the patients (#5, 6) had situs inversus. The mutational status was known for one patient (patient #4), who was a compound heterozygote for two mutations in DNAH5. Non-PCD controls were individuals referred to the Institute of Tuberculosis and Lung Diseases in Rabka for various reasons unrelated to PCD; they presented no symptoms of acute pulmonary disease, no bronchoscopic signs of the disturbance of mucociliary transport, and had normal ciliary beating in the light microscopy.
Biopsies were obtained during a diagnostic bronchoscopy, using the same protocol in PCD and non-PCD individuals. Specimens were stored in RNAlater buffer and used for RNA isolation.

Gene expression analysis
The quality and concentration of RNA were determined using the 2100 Bioanalyzer (Agilent, Amstelveen, The Netherlands) with the Agilent RNA 6000 Nano Kit. Anti-sense RNA was synthesized, amplified and purified using the Ambion Illumina TotalPrep Amplification Kit (Ambion, USA) according to the manufacturer's protocol. Complementary RNA was hybridized to Illumina HumanRef-12 Whole Genome BeadChips and scanned on the Illumina BeadArray Reader. Data was handled through the Illumina BeadStudio Gene Expression module v3.2. The expression data has been submitted to GEO under accession number GSE25186.

Data-analysis
Initial steps of data preprocessing, quantile normalization, background subtraction and quality control were performed using the BeadStudio software. The data were exported to GeneSpring 9. Expression values with a threshold value below 5 were raised to 5, and the data were scaled by a base 2 logarithm.
Ensembl Biomart module was used to obtain Ensemble gene IDs of the differentially expressed genes that were subject to annotation enrichment analysis and the Ciliome database (www. ciliome.com) [37] search. Fold-change analysis for all the probes was performed using GeneSpring 9 software. Genes that had at least 2-fold expression change (with the t-test p-value ,0.05) in PCD patients as compared to the controls were subjected to gene annotation enrichment analysis using DAVID (Database for Annotation, Visualization, and Integrated Discovery) [38]. DAVID applies a modified Fisher exact test, to establish if the proportion of genes falling into an annotation category differs for a particular group of genes and the background group of genes; the p-values were corrected for multiple testing using the FDR method. The resulting annotation terms were then clustered into functional groups.
PubMed was used as references for known ciliary genes encoding components of radial spokes and intraflagellar transport proteins. Outer and inner dynein arms genes were taken from Pazour et al [39]. The R package was used for statistical calculations.
The log2-transformed mean expression values for the probes representing protein components of the specific ultrastructural ciliary elements were compared in the patients and controls, and the Wilcoxon Mann-Whitney test was used to obtain p-values for individual comparisons. In case of genes represented by multiple probes, mean values were used for calculations. The unweighted Z-method of combining probabilities from independent tests was used to assess if the distribution of p-values deviated from normal [40]. The Z-transform test converts one-tailed p-values from each of the independent tests into standard normal deviates; the sum of these standard normal deviates, divided by the square root of the number of tests, has a standard normal distribution if the common null hypothesis is true.

Whole genome analysis (fold change and gene enrichment)
We identified 1323 genes (1396 probes), which displayed at least 2-fold difference in the expression level in PCD compared to non-PCD bronchial tissue (t-test p-value ,0.05); 721of these genes (755 probes) were over-expressed, whereas 602 genes (641 probes) showed a reduced expression in PCD (Table S2).
To investigate the biological meaning of these differentially expressed genes, we performed annotation enrichment analysis, separately for the up-regulated and down-regulated genes, using DAVID software [38]. Annotations of the genes with the reduced expression in PCD patients were significantly enriched for terms related to microtubules, motility, cilia, and dyneins (Table 1). In the group of up-regulated genes, the highest scoring cluster was related to cell cycle and mitosis (Table 2). To investigate how many of the down-and up-regulated genes were functionally linked to cilia, the two groups (genes from Table S2) were used to query the Ciliome database [37]. This search revealed that 19% of the down-regulated and 11% of the up-regulated genes were related to cilia. Under a theoretical assumption that approximately 10% of the human genes are related to cilia (the Ciliome database contains 2024 entries), the proportion of ciliary genes found in the down-regulated group was significantly higher than in the whole set of human genes (binominal p-value 1.69e-08), while the list of up-regulated genes was not significantly enriched for genes related to cilia (p-value 0.26).

Known ciliary genes
To search for the evidence of differential expression of the genes related to specific elements of the axonemal ultrastructure, we   Table 3. Differential expression of the genes encoding outer dynein arms, inner dynein arms, radial spokes and intraflagellar transport proteins in the group of PCD cases compared to the non-PCD controls. assembled a list of human proteins represented on the microarray and predicted to be a part of the known ultrastructural ciliary components (ODA, IDA, radial spokes, intraflagellar transportrelated proteins). We used the data from Pazour et al. [39], who applied phylogenetic analysis to assign human heavy dynein chains to ultrastructural components of outer and dyneins arms of flagella, and from Hom et al. [41], who compiled data on the nomenclature of dyneins. The assignment of the remaining genes was based on the PubMed search. The resulting list of 37 human genes encoding proteins predicted to be a part of known ciliary components included: 8 ODA genes, 10 IDA genes, 8 radial spoke genes and 11 intraflagellar transport genes ( Table 3). The log 2-transformed mean expression levels in PCD and non-PCD bronchial tissues, the fold change, the log2 of the fold change and the corresponding Wilcoxon Mann-Whitney test p-values are presented in Table 3. All the genes had a lower expression in PCD than in controls. Statistically significant evidence of the down-regulation in PCD was also observed when the combined p-values were calculated, using an unweighted Z-method [33], in four subgroups of the genes encoding ciliary components: ODA (p = 3.5e-08), IDA (p = 3.19e-12), radial spokes (p = 1.06e-4) and intraflagellar transport proteins (p = 1.27e-07), respectively ( Table 3).
The results were not dependent on individual PCD patients. Removing of each of the PCD subject did not significantly change the result of gene annotation analysis. Evidence of the downregulation in four subgroups of the genes encoding ciliary components remained also significant.

Discussion
A combination of the fold change with the significance measured by an uncorrected t-test (which takes into consideration the variance between the samples), has been shown to be a reliable and reproducible measure of differential expression [42].The gene expression profiling using the whole-genome Illumina panel, performed in bronchial biopsies from six PCD patients and nine non-PCD controls, revealed 1396 probes (representing 1323 genes) with at least 2 fold expression level difference between the PCD and non-PCD groups (t-test p-value ,0.05). Functional annotation of the differentially expressed genes revealed that those down-regulated in PCD were related to microtubules, motility and cilia, and the up-regulated ones to the cell cycle and mitosis, i.e. to the processes that are related to microtubules by mitotic spindle and centrioles.
To further characterize the expression changes in PCD, we analyzed expression levels of a number of cilia-related genes represented on the chip. We assembled the list of 37 human genes representing: three subgroups, related to the ultrastructural components of motile cilia (outer dynein arms, inner dynein arms, and radial spokes), plus one functional subgroup (intraflagellar transport proteins). In all these subgroups, the reduced expression in bronchial biopsies from PCD patients was observed. The largest number of probes exhibiting significantly decreased expression was observed in the group of genes encoding components of inner and outer dynein arms.
Of the twenty-six genes known to be mutated in PCD, CCDC39, CCDC164 and SPAG1 were not represented on the array and TXNDC3 was not stably expressed. Of the remaining twenty-two stably expressed genes (Table 4), ten were significantly downregulated (.2-fold change and p,0.05). The majority of the down-regulated PCD genes encoded components of cilia ultrastructure: ODA (DNAH11, DNAH5, DNAI1, DNAI2, DNAL1, CCDC114); central pair appendix (HYDIN). The remaining down- Table 3. Cont.   Table 4. Differential expression of the genes mutated in PCD and two genes mutated in syndromic disorders associated with PCD symptoms (RPGR and ODF1). regulated PCD genes encoded cytoplasmic proteins required for the assembling of dynein arms (LRCC50 alias DNAAF1, ZMYND10). In twelve of the PCD genes the expression level did not significantly differ between PCD and non-PCD biopsies; five of these genes (KTU alias DNAAF2, CCDC103, DYX1C1, c21orf59 and C19ORF52 alias DNNAAF3) encoded cytoplasmic assembly proteins, two (CCDC40 and CCDC65) encoded N-DRC (nexindenein regulatory complex) proteins, one (HEATR2) encoded protein of an unknown location and three encoded radial spokes proteins (RSPH4A, RSPH1 and RSPH9), one (ARMC4) encoded ODA docking component; only the expression of DNAAF3 and CCDC40 was slightly higher in PCD than in healthy controls. ODF1 associated with the syndromic form of PCD, was expressed at the significantly lower lever in PCD patients. DYX1C1, LRCC6, RPGR, RSPH1 had the fold change .2 but the p-values were between 0.05 and 0.1. Given the lack of linear relation between RNA and protein, direct conclusions about a possible correlation between the RNA level and ultrastructural ciliary defects cannot be made. Large multi-protein complexes of the dynein arms are pre-assembled in the cytoplasm, before their intraflagellar transport to the ciliary axonemes. Immunohistochemical staining of the ciliated epithelium from PCD patients with identified mutations has shown that a defect/absence of one of the ciliary proteins may impair the delivery or assembly of other, non-mutated members of the axonemal substructures [43]. This relationship has been shown e.g. for mutated DNAH5 or DNAI1 and DNAH9 [43], mutated KTU and components of ODA/IDA [16], mutated CCDC40 or CCDC39 and components of IDA [19,20], mutated LRRC6 and components of ODA/IDA [23]. Our data suggest that this impairment of the ciliary ultrastructure assembly may in fact originate at the transcription level; we hypothesize that a regulatory system may exist, which -in the presence of a mutated ciliary gene -down-regulates the expression of at least part of the ciliary genome. Such a regulation would play an important role in the molecular pathomechanism of PCD. During ciliogenesis the expression of ciliary genes is regulated by common transcription factors. The FOXJ1 transcription factor and RFX (regulatory factor X) family of transcription factors have been shown to regulate the expression of many cilia related genes. We have found FOXJ1 (FC = 2,086838) and RFX2 (FC = 2,412601) among the down-regulated genes indicating that the putative regulatory event acts upstream of these two transcription factors [44].
In our study, the strongest evidence of a significant downregulation of the cilia-related genes in PCD bronchial biopsies concerned the genes predicted to encode elements of the inner and outer dynein arms and to a lesser extent, elements of radial spokes and proteins involved in ODA assembly or in the intraflagellar transport. The small study group provides only the proof of the concept and further studies are warranted. For example, to estimate the specificity of a putative regulatory system it will be important to repeat similar analyses using groups of patients carrying mutations in the same gene. In our small group of patients, only one of the patients had a known mutational status (two complementary heterozygous mutations in DNAH5). DNAH5 is the strongest contributor to the genetic background of PCD, but one cannot exclude the possibility that there were mutations in other genes among the patients we examined. One may envision that the homogeneous group of patients with the mutations in the specific genes would cause different groups of the ciliary genes to undergo down-regulation.
Apart from the contribution to the understanding of the mechanisms of ciliary assembly, such studies may contribute to finding further candidate genes in PCD. To date, a number of candidate disease genes in several genomic regions reported to be linked to PCD remain to be identified; the failure of many linkage studies may be ascribed to the noise in linkage data, and difficulties in the functional prioritization of the candidate genes. In this context, observation of a down-regulation of the particular gene's expression in PCD might provide clues for selecting candidate genes from the linkage peaks.

Supporting Information
Table S1 Electron microscopy, NO measurements and situs status in the PCD patients studied. (DOCX)