Identification of candidate regulatory elements controlling transcriptome during the formation of interphalangeal joints.

doi:10.21203/rs.3.rs-1413327/v1

PDF

Research Article

Identification of candidate regulatory elements controlling transcriptome during the formation of interphalangeal joints.

https://doi.org/10.21203/rs.3.rs-1413327/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Formation of the synovial joint starts with the visible emergence of interzone, a stripe of densely packed mesenchymal cells located between distal ends of the developing skeletal anlagen. Recently the transcriptome of the early synovial joint was reported. Knowledge about regulatory elements such as enhancers would complement these data and lead to a better understanding of the genetic control of gene transcription at the onset of joint development. Using ChIP-sequencing we have mapped the H3 signatures H3K27ac and H3K4me1 to locate the regulatory elements specific for the interzone and adjacent phalange, respectively. This atlas of candidate enhancers (CEs) was then used to map the association between these respective joint tissue specific CEs and biological processes. Subsequently, an integrative analysis of transcriptomic data and CEs identified new putative regulatory elements of genes expressed in interzone (e.g. GDF5, BMP2 and DACT2) and phalange (e.g. MATN1, HAPLN1 and SNAI1). We also linked such CEs to genes characterized as crucial in synovial joint hypermobility and osteoarthritis, as well as phalange malformations. These analyses show that the enhancer atlas can serve as resource for identifying putative disease-causing genomic regulatory regions in patients with synovial joint dysfunctions and/or phalange disorders, and enhancer-controlled synovial joint and phalange formation.

chondrocytes

enhancers

interzone

joint development

phalange

Synovial joints, organs present at the articular ends of long bones, are essential for mobility in bony vertebrates. They comprise of articular cartilage, synovium, ligaments and the synovium capsule¹ . Due to their function, the joints are frequently exposed to mechanical stress and thus prone to injuries. Congenital malformations and a number of diseases affect joint structure, thereby causing a decrease of joint functionality. For instance, misexpression of PITX1 caused by enhancer adoption results in dysplastic elbow joints in Liebenberg syndrome (OMIM #186550)², a homozygous mutation in IMPAD1 leads to chondrodysplasia with joint dislocations (OMIM #614078)³, and loss of EXOC6B causes joint dislocations and defects in joint mobility, characteristic for patients with spondylo-epimetaphyseal dysplasia with joint laxity, type 3 (OMIM #618395)⁴. Osteoarthritis (OA) is the most prevalent synovial joint disease affecting adults⁵. Typical in OA is the progressive degeneration of articular cartilage and accompanying subchondral bone sclerosis, joint space narrowing, osteophyte formation and the variable degree of synovium inflammation⁶, eventually causing joint destruction. The latter frequently needs intervention by joint replacement⁷.

A comprehensive understanding of gene regulatory networks (GRNs) orchestrating synovial joint formation will contribute to the understanding of both healthy and pathological processes taking place in this organ. Such knowledge will help the development of novel cell and/or gene-based strategies for treatment of injured articular cartilage, particularly these based on the concepts of developmental engineering⁸. The first morphologically distinguishable event in joint development is the formation of interzones, with distinct progenitor cells giving rise to the majority of articular tissues⁹. The condensing cartilage anlage, at the locations of future joints, undergoes several rounds of cell proliferation¹⁰. Nonetheless, the discovery that the influx of cells from outside of the interzone contributes to overall increased interzone cell density, points to an important mechanism in interzone formation^11,12. At the same time, SOX9 (a member of Sry family of transcription factors) expression becomes repressed, arresting the chondrogenic program and allowing the interzone to form¹³. The final structure is comprised of two layers of cells, named the outer and intermediate interzone layers, with differentially expressed genes (DEGs), including COL2A1 (encoding α(II)-collagen) and MATN1 (Matrilin-1), which have higher expression in the outer interzone, and GDF5 (Growth and Differentiation Factor-5, a ligand of the BMP subgroup of the TGFβ family) with higher mRNA level in the intermediate interzone¹⁴. These layers will contact the ends of future bones, while an inner cell layer of yet to be defined function is also present.

Cells within the cartilage anlage change their phenotype progressively from round to columnar, pre- and eventually hyper-trophic chondrocytes, contributing to longitudal cartilage and bone growth (Green et al., 2015). Acknowledged molecular markers for the round chondrocytes include the aforementioned MATN1 and COL2A1¹⁵. Additionally, chondrogenic induction and differentiation is accompanied by the expression of RUNX2 (encoding a RUNT family transcription factor (TF)) followed by the expression of COL10A1 (α(X)-collagen). The latter is a specific marker gene for hypertrophic chondrocytes, and some of its enhancers have been elucidated¹⁶. It has been suggested that the round chondrocytes may contribute to the articular cartilage, but their contribution to the interzone structures remains unclear and may depend on restricted exposure to BMP and/or WNT signals¹¹. Thus, while new knowledge is emerging regarding early stages of joint formation, including interaction with adjacent cartilage, the insight into gene expression control within cells of the joint interzone remains incomplete. Indeed, while numerous reports have described the transcriptome in synovial joint formation^14,17, including at single-cell level¹⁸, relatively little is still known about the activity of enhancers during that process.

Enhancers regulate gene transcription mainly in cis, within the topologically-associating domains (TADs), where they promote intra-TAD control of transcription of loci by making TFs and co-factors bridge between their bound distal enhancer sites and the promoter-proximal region of the appropriate target gene(s), hence achieving physical proximity^19-21. Enhancers are associated with histone modification signatures, such as H3K27ac and H3K4me1, chromatin accessibility, high-affinity binding of co-factors (e.g. histone-acetylation containing P300), TFs and RNA-Pol2²². Defects in enhancer function have been linked to limb malformations as well, for example in Hass-type polysyndactyly (OMIM #186200)²³, split hand/foot malformation (OMIM #183600)²⁴, Leri-Weill dyschondrosteosis (OMIM #127300) and Laurin-Sandrow syndrome (OMIM #135750)²⁵.

Here, we focus on identifying regulatory elements (REs) that are active in the joint interzone and adjacent phalange, respectively, permitting us to develop a candidate enhancer (CE) atlas. In combination with the available transcriptomes, such an atlas is expected to help elucidating the molecular mechanisms that control joint interzone formation and/or cause joint disease. We opted for microsurgical dissection of interzone as opposed to using Gdf5-positive (+) cell selection procedures, because not all cells during early stages of joint formation are convincingly Gdf5+^12,18. We identified unique interzone/phalange CEs that are conserved between chicken, mouse and human, and functionally annotated these CEs, followed by integrative analysis of cell-type specific CEs and DEGs. We also associated the CEs with synovial joint and phalange abnormalities, and a higher risk of OA.

Microdissection of joint interzones and phalanges

The interzones and the adjacent proximal part of phalange were dissected from the third digit of the hindlimb of chick embryos (stage HH32, when the interzone was distinguishable under the microscope; Fig. 1a). Next, the RNA-sequencing was performed on the separated tissues (see Materials and Methods).

Following the bioinformatic analysis, a list of DEGs was compiled. We confirmed that the interzone samples had increased mRNA levels of GDF5, ENPP2, COL3A1 and ERG, each already known to be expressed in joint interzones. Phalange samples had significantly higher expression of COL2A1, MATN1, SNAI1 and RUNX2 (Fig. 1b; Supplementary Table 1). Notably, COL10A1 mRNA expression was not detected in phalanges (Supplementary Table 2), supporting the notion that the phalanges were collected at the early developmental stage, prior to hypertrophy. We further validated differential expression of selected interzone markers (i.e. GDF5, ENPP2, ERG) using RT-qPCR (Fig. 1c), which confirmed our RNA-seq data. These results show that the we successfully dissected tissues of interest.

An enhancer atlas of joint interzone and phalange identifies candidate REs involved in the regulation of cell identity

In the separated joint and phalange we also mapped the global H3K27ac and H3K4me1 signatures by chromatin immunoprecipitation followed by next-generation sequencing (ChIP-seq). The unsupervised clustering analysis revealed that interzone and phalange have distinct profiles for both H3K27ac and H3K4me1 (Fig. 2a,b). This mapping of the regions enriched for these H3 modifications enabled us to generate a joint/phalange CE atlas: regions enriched for both H3K27ac and H3K4me1 were denoted strongly active enhancers, for H3K27ac active enhancers, and for H3K4me1 poised enhancers. Since enhancers are often evolutionarily conserved²⁶, we decided to select regions conserved among chick, mouse and human. In doing so, we identified 14,217 strongly active, 5,479 active and 11,913 poised enhancers in the interzone, and 14,224, 6,041 and 12,997, respectively, in the phalange (Supplementary Table 3; for their frequency and similar ratios in both samples, see Fig. 3a).

For functional annotation of these REs we used GREAT, which extracts gene ontology (GO) terms linked to biological processes²⁷. The denoted strongly active enhancers associated with cartilage and skeletal development (Fig. 3b), whereas active and poised enhancers mostly linked to general cell functions or processes not specific for skeletal development (Supplementary Fig. 1). Our results also matched with Cheung et al.²⁸ who retrieved GO terms that associate with cartilage and skeletal development of strongly active enhancers of human chondrocytes, and more general GOs of active and poised enhancers.

This agreement with data of Cheung et al.²⁸ prompted us to proceed exclusively with strongly active CEs. These CEs are typically located >5 kb away from the respective transcription start site (TSS) (Fig. 3c). Among them, we confirmed previously characterized enhancers of well-studied loci (Fig. 4a), specifically those expressed in interzone (e.g., GDF5)²⁹ or chondrocytes (e.g., phalangeal IHH, SOX9, ACAN)^30-32. In parallel, and further confirming the validity of our in silico selection approach, we extended the analysis by using Vista Enhancer Browser dataset³³, leading to identification of 257 enhancers (Supplementary Table 4), which have been functionally validated during embryogenesis (for illustration of 6 of these, see Fig. 4b).

To further characterize strongly active CEs of developing interzone and phalange, we selected the mutually exclusive strongly active enhancers, yielding 3,406 CEs (out of the aforementioned 14,217 in total) unique for interzone and 3,407 (out of 14,224) for phalange (Supplementary Table 5). GREAT linked many of such interzone-specific CEs to mesenchymal cell differentiation, and regulation of transmembrane receptor protein serine/threonine kinase signaling (Supplementary Fig. 2, two top panels). In contrast, CEs exclusive for phalange retrieved GO terms including chondrocyte differentiation and endochondral bone morphogenesis (Supplementary Fig. 2, two bottom panels).

Next, we investigated whether the change of enhancer state from strongly active to poised would be relevant to the regulation of tissue-specific genes. Indeed, 2,111 changes occur with strongly active CEs (out of the aforementioned identified 14,217 in total) in interzone and are found poised in phalange; in comparison, 1,502 changes occur with strongly active CEs (out of 14,224 in total) in phalange and are found poised in interzone (Supplementary Table 6). Strikingly, the interzone strongly active enhancers that are poised in phalange were found to associate with regulation of transmembrane receptor protein serine/threonine kinase signaling (Supplementary Fig. 2), which is consistent with our functional annotation of interzone-specific CEs. In contrast, phalange strongly active enhancers that are poised in interzone linked to positive regulation of cartilage differentiation (Supplementary Fig. 2).

Integrative analysis of DEGs and CEs

To correlate the transcribed genes with CEs, we superimposed our RNA-seq and ChIP-seq data. First, DEG analysis revealed 116 upregulated genes in interzone, and 61 such genes in phalange (log2FC > 0.5, p.adj < 0.05) (Supplementary Table 2). Pathway enrichment analysis of these showed that genes upregulated in interzone again linked to transmembrane receptor protein serine/threonine kinase signaling (Fig. 5a), in line with our preceding annotation of interzone-exclusive CEs (see above). In particular RASL11B, LTBP1, TGFB2, GDF5, FSTL1, BMP2, DACT2, CCN3, BMP6, CILP, INHBB and BMPR2 were found upregulated in interzone as compared to phalange (Fig. 5b). Similarly, analysis of genes upregulated in phalange linked these to chondrocyte differentiation and also endochondral bone morphogenesis (Fig. 5c), which is consistent with functional annotation of phalange-specific CEs. The genes involved it these two processes are RUNX2, COL2A1, TRPV4, COL27A1, MATN1, COMP and CYTL1 (Fig. 5d).

Next, to associate the CEs to all DEGs significantly upregulated in interzone, we extracted strongly active enhancers from genomic regions located within -/+ 1 Mb from the TSS of the annotated genes, resulting in the identification of 857 interzone-specific CEs (Supplementary Table 7). Examples of enhancer analyses are presented in Fig. 6a. For phalange-upregulated DEGs we identified 547 phalange-specific CEs (again within -/+ 1 Mb from the TSS of the genes; Supplementary Table 8), with examples given in Fig. 6b. Collectively, the integrative analysis of transcriptome data with CEs assignment, and considering interzone vs. phalange signatures, showed that the DEGs involved in celltype-specific processes are regulated by cell-specific CEs.

CEs regulate skeletal malformation and disease-relevant genes, and are associated with a higher risk of OA

Mutations in genes and REs have been linked to various limb malformations and skeletal defects^21,23-25. We applied two types of analysis to screen for CEs that link to molecular etiology of limb disorders in general. First, we assigned our strongly active CEs to the proximal genes (including relevant respective marker genes and DEGs) and tested whether these genes have been associated previously with limb phenotypes, either in patients (including in syndromes) or mouse models. Analysis of interzone/phalange specific strongly active CEs showed that these are indeed involved in the regulation of genes linked to joint and phalange abnormalities (Supplementary Table 9a-o).

The interzone-specific, strongly active CEs particularly associate with defective joint mobility in humans (Fig. 7a). For instance, we identified such putative enhancers of OTX2 and TGFB2, which are genes that have been linked to joint laxity (OMIM #610125 and #614816, respectively); candidate CEs of FLNB, a gene associated with joint dislocation and carpal fusion (OMIM #150250 and #272460, respectively); we also predicted enhancers of COL5A1, a gene linked to joint hypermobility (OMIM #130000) (Supplementary Table 9a). In mouse, the interzone-specific CEs associate with abnormal joint morphology and fused joints (Fig. 7b; see also Supplementary Table 9b-c). Phalange-specific CEs have in humans been linked to aplasia/hypoplasia of the phalanges, short phalanges, and abnormality of the phalanges of the toe (Fig. 7c). For example, we identified putative enhancers of BMPR1B and IHH (Supplementary Table 9d), both being genes associated with brachydactyly type A (OMIM #112500), and candidate enhancers of RUNX2, a gene linked to cleidocranial dysplasia, with brachydactyly (OMIM #119600) (Supplementary Table 9f). The phalange-specific CEs have also been linked to abnormal chondrocyte and cartilage morphology, chondrodystrophy, abnormal bone ossification and short limbs in mouse (Fig. 7d; see also Supplementary Table 9j-o).

Next, we screened the GWAS catalog (NHGRI-EBI³⁴), which resulted in identification of 3,263 single-nucleotide polymorphisms (SNPs) within CEs, with 232 of these linking to skeleton related traits (Supplementary Table 10). For instance, we identified single-nucleotide variations (SNVs) within CEs that have been associated with OA relevant genes, such as for ALDH1A2 (rs4775006, P-value 8x10^-10)^35,36 and WWP2 (rs34195470, p-value 3x10^-13)^5,37. Additionally, we identified SNVs associated with increased risk of OA, which are located within CEs that map to LRIG3 (rs79056043, P-value 1x10^-9), CRADD (rs7953280, P-value 5x10^-12) and ROCR (rs8067763, P-value 2x10^-9)^5,38. Altogether, the analysis of the tissue-specific CEs showed the association with either synovial joint or phalange congenital abnormalities, which affect their function, as well as the identified regulatory elements of genes relevant for joint degenerative disorders, in particular OA.

Transcriptome analysis of joint interzone gained attention in recent years^14,18,39 but unlike in other fields, relatively little is known about the cis-regulatory elements involved in the establishment of interzone during limb development. Here, we carried out an integrative analysis of transcriptomic and epigenetic data, and subsequently generated a ChIP-seq based CE atlas, for separated interzone and phalange, respectively. For this, we optimized a fast-dissection protocol for careful collection of joint interzone samples. Using both RNA-seq and validation RT-qPCR with selected markers we showed that such collected interzones have significantly higher expression of GDF5, ENPP and ERG as compared to adjacent phalange. At the same time, we validated the dissection protocol for collecting interzone cells for genome-wide experiments that require high numbers of input cells. This optimization of sample collection and separation permitted for the identification of CEs and correlation of the latter with digit/joint formation. Functional annotation of the CEs illustrates that regions enriched in both H3K27ac and H3K4me1 are associated with genes important for cartilage/skeletal development, unlike the CEs enriched in one of the two H3 marks only. We also identify CEs unique for interzone/phalange linked to specific biological processes, which correlate with pathway enrichment analysis of DEGs.

One of the pathways enriched in upregulated genes and CEs in interzones is transmembrane receptor protein serine/threonine kinase signaling, i.e. TGFβ/BMP family signaling. This signaling system has been well-characterized in the process of chondrogenic differentiation^40-43, but is still not well characterized in joint interzones. Suppression of BMP activity in the interzone region is essential for normal joint development^11,44. In contrast, GDF5 is expressed at high level in the interzone region, and GDF5 null mutation results in joint defects⁴⁵. In our study we confirm in an alternative way previously described enhancers of GDF5 and identify novel CEs of GDF5. Another example of a BMP-upregulated gene in the interzone region, consistent with reported in situ RNA-hybridization, is BMP2⁴⁶. BMP2 is involved in joint maturation; its genetic inactivation in synovial joint forming cells results in changes in extracellular matrix and also shape of the meniscus⁴⁷. Interestingly, within the genomic regulatory landscape of BMP2, we identify several CEs likely to be active in the interzone, but this will require further investigation in future projects. We also document a change of enhancer state between strongly active and poised enhancers, which enabled us to identify CEs associated with unique biological processes. These processes (some, but not all) correlate with analysis of interzone/phalange specific CEs and pathway enrichment analysis of certain DEGs.

When considering changes in gene expression as causal for limb malformation, it becomes necessary to also include characterization of enhancers which may drive misexpression of disease-causing genes. We show that many CEs associate with genes important for normal development as well as etiology of both synovial joints and phalanges. For instance, in our study we identified CEs of GDF5 (as shown by Chen et al.²⁹) to be active in the interzone region. Mutations in GDF5 lead to joint malformations, and in a genomic region containing GDF5 these experimentally tested enhancers have been linked to higher risk of OA^29,48-50. We also characterize interzone-specific CEs of the BMP2 locus. Importantly, BMP2 conditional knock-out mice develop progressive OA in the knees⁴⁷. Another set of gene-linked enhancers in joint abnormalities are CEs for FSTL1, associated with rheumatoid arthritis⁵¹. We also identify CEs linked to genes (e.g., OTX2, TGFB2, COL5A1) associated with defects in joint mobility.

We used the GWAS catalog NHGRI-EBI to identify SNVs located within CEs and that are associated with higher risk of OA. This yielded CEs linked to the genes previously described in OA, such as ALDH1A2 and WWP2. In addition to interzone CEs, we also characterized several phalange CEs linked to the genes important in chondrocyte-related disorders. For instance, we identified CEs associated with RUNX2 and IHH, both also pivotal genes in the molecular etiology of limb malformation, including brachydactyly^52,53. Altogether, this illustrates that our CE atlas provides information on association of CEs with already existing gene-to-disease correlations. This will be helpful in studies of possible variations in genomic regions in patients without mutations in the protein-coding genes.

We also discovered that BMPR2 is co-expressed with GDF5 in interzones. Therefore, it is possible that BMP signaling in the interzone region is prevented at the intracellular level, or BMP ligands play a dual function, mutually exclusive in the interzone cells and chondrocytes. Both GDF5 and BMPR2 are upregulated in chondrocyte dedifferentiation in vitro⁵⁴, which could support the hypothesis that interzone originates from dedifferentiated chondrocytes⁹. WNT signaling was shown to be crucial in the formation of joint interzone. We identified that DACT2 (encoding a WNT inhibitor) is significantly upregulated in interzone as compared to phalange, which is consistent with reported in situ RNA-hybridization⁵⁵. This shows the complexity of WNT signaling regulation during joint formation. Within the genomic regulatory landscape of DACT2 we identified several CEs likely to be active in the interzone. However, the more detailed characterization of these DACT2 enhancers in the interzone requires additional studies.

Our work identified several genes from the TGFβ family ligands or system components upregulated in interzone, for instance TGFB2, LTBP1 and INHBB. Several studies showed antagonistic action of TGFβ/Activin-Nodal pathways on BMP signaling in several cell types^56,57. Therefore, it is tempting to suggest that the upregulation of TGFβ family system components in interzone may have effects on attenuation and/or inhibition of BMP signaling within prospective sites of synovial joint formation.

Tissue collection

All vertebrate animal experiments (chick early embryos) were carried out in accordance with the relevant guidelines as applied and approved by the Ethical Committee at the Medical University in Lublin, where this work was performed, and also comply with the European regulations (directive 2010/63/EU). The tissue collection was performed on the chicken embryos until 7.5 days post fertilization which is exempt from the Ethical Committee Approval. Chick White Leghorn embryonated eggs were incubated at 38.5ºC in fixed humidity for 7.5 days, followed by evaluation of developmental stage based on the Hamilton Hamburger classification (HH), using Zeiss Stereo Discovery V8 microscope equipped with 0.63x Plan Apo S Objective Lens. Selected embryos at HH32 were sacrificed for tissue microdissection. The joint interzones and adjacent phalange samples were microdissected from hindlimb digit-3 using Dumont No.5 forceps (tip dimensions: 0.005 x 0.025 mm)

RNA extraction

Total RNA was extracted and prepared using Syngen Tissue RNA Kit followed by DNA digestion with QIAGEN RNase-free DNase Set. The RNA quality was validated on 1% agarose gel (for RT-qPCR) or on Agilent 2100 Bioanalyzer system with RNA 6000 Nano Assay (for RNA-seq). All latter samples had a RIN value > 9.0.

RT-qPCR

Three independently extracted RNAs from both interzone and phalange were reverse transcribed to cDNA using Invitrogen™ SuperScript™ IV Reverse Transcriptase and Oligo(dT) primer. The qPCR was performed using PowerUp™ SYBR™ Green Master Mix II on LightCycler® 480 Instrument II. Gene expression was normalized to expression of GAPDH. Statistics were computed using Mann-Whitney-Wilcoxon Test. The list of used primers is given in Supplementary Table 11.

RNA-seq library preparation

The 6 samples (3 of interzone, 3 of phalange- independent biological replicates) were prepared with the Smart-seq2 method (Picelli et al., 2013). In brief, poly(A)-RNA was reverse transcribed using oligo(dT) primers. Template switching by reverse transcriptase was achieved by using an LNA-containing TSO oligonucleotide. The reverse-transcribed cDNA was pre-amplified with primers for 18 cycles, followed by clean-up. Tagmentation was performed on 500 pg of the pre-amplified cDNA with Tn5 followed by gap repair. The tagmented library was extended with Illumina adaptor sequences by PCR for 14 cycles and purified. The resulting sequencing library was measured on Bioanalyzer and equimolar loaded onto a flowcell and sequenced according to the Illumina TruSeq v3 protocol on the HiSeq2500 with a single- read 50 bp and dual 9 bp indices.

RNA-seq data analysis

The fastq files were checked for quality using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) followed by removal of adapters using Trimmomatic⁵⁸. Further, reads were mapped to the Gallus gallus 6.0 reference genome using STAR with default parameters⁵⁹. Gene expression values were called using featureCounts with Ensemble release 104 annotation⁶⁰. The differential data analysis has been performed with DEseq2⁶¹, and heatmaps have been created in R environment for statistical computing.

Chromatin-immunoprecipitation (ChIP)

For each sample, either 100 interzones or phalanges were dissected from the 3rd hindlimb digit and pooled together. Further, tissues were dissociated for 3 hours at 37ºC using 2.4% Collagenase-II (Gibco™) resuspended in DMEM/high-glucose medium containing 10% fetal bovine serum (FBS). Cells were passed through a 40-μm cell strainer (BD Falcon) and then counted using a hemocytometer. 10⁶ cells were cross-linked using 1% formaldehyde at room temperature (RT, 24ºC) for 9 min. Fixation was quenched with ice-cold 0.125 M glycine for 5 min at 4ºC. To remove excess formaldehyde, two rounds of centrifugation of the cells followed by resuspension in ice-cold PBS were carried out. Next, the cell nuclei were isolated using ice-cold nuclei extraction buffer (NEB) containing 10mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA, 1% IGEPAL® CA-630 and complete protease inhibitors (Roche). Subsequently, the nuclei were resuspended in SDS-containing lysis buffer (50 mM Tris-HCl pH 8.0, 10 mM EDTA and 1% SDS), and chromatin was sheard obtaining the average size of 150 bp in AFA Fiber Pre-Slit Snap-Cap (130 µl) microtube using a S220 Focused-ultrasonicator.

For ChIP, 500 ng of sonicated chromatin was immunoprecipitated with 7.5 μg of anti-H3K27ac antibody (Active Motif, Cat. No. 39133) or anti-H3K4me1 (61781). Input sample was collected prior to immunoprecipitation reaction. Chromatin pre-cleaning incubation with protein-A and protein-G agarose beads (Millipore) was carried out in immunoprecipitation buffer (50 mM Tris-HCl pH 8.0, 0.15 M NaCl, 1 mM EDTA pH8.0, 1% Triton X-100 and 0.1% sodium deoxycholate) for 3 hours at 4ºC while rotating. In parallel, the antibodies were incubated with previously blocked A- and G-agarose beads (Millipore) also for 3 hours at 4ºC, again while rotating.

After pre-cleaning, the chromatin was mixed with the pre-bound antibodies with A- and G-agarose beads and incubated overnight (O/N) at 4ºC, rotating. The next day, multiple rounds of washes of the beads were conducted. Each wash was carried for 10 min at 4^ºC while rotating. The beads were washed once with RIPA-150 buffer (50 mM Tris-HCl pH8.0, 0.15M NaCl, 1mM EDTA pH8.0, 0.1% SDS, 1% Triton X-100, 0.1% sodium deoxycholate), twice with RIPA-500 (50 mM Tris-HCl pH8.0, 0.5 M NaCl, 1 mM EDTA pH8.0, 0.1% SDS, 1% Triton X-100, 0.1% sodium deoxycholate), once with RIPA-LiCl (50 mM Tris-HCl pH8.0, 1mM EDTA pH8.0, 1% Nonidet-P40, 0.7% sodium deoxycholate, 0.5M LiCl) and twice in TE buffer (10 mM Tris-HCl pH8.0, 1mM EDTA pH8.0). Subsequently, the chromatin was eluted with 200 μl of fresh elution buffer (1% SDS and 0.1 M NaHCO₃) followed by addition of 100 μl of TE buffer and 25 μl of 5 M NaCl prior to reverse-crosslinking at 65^ºC for 16 hours. The next day, chromatin was incubated with 2 μl of Proteinase-K (10 mg/ml) for 1 hour at 56^ºC, and 2 μl of RNaseA (10 mg/ml) for 45 min at 37ºC, and DNA was further purified using QIAquick PCR Purification Kit (QIAGEN). The size distribution of immunoprecipitated fragments was evaluated using Agilent 2100 Bioanalyzer system with High Sensitivity DNA. Additionally, the DNA-concentration of input and immunoprecipitated samples was measured on Qubit 2.0 Fluorometer (Invitrogen).

ChIP-sequencing

ChIP-seq libraries were prepared using QIAseq Ultra Low Input Library Kit (QIAGEN, Hilden, Germany). Briefly, DNA was end-repaired, adenosines were added to the 3′ ends of dsDNA and adapters were ligated (adapters from NEB, Ipswich, MA, USA). Following the adapter ligation, uracil was digested by USER enzyme from NEB (Ipswich, MA, USA) in a loop structure of the adapter. Adapters containing DNA fragments were amplified by PCR using NEB starters (Ipswich MA, USA). Library quality evaluation was done with Agilent 2100 Bioanalyzer using the Agilent DNA High Sensitivity chip (Agilent Technologies, Ltd.) Quantification and quality evaluation of obtained samples were done using Nanodrop spectrophotometer (Thermo Scientific, NanoDrop products, Wilmington, USA), Quantus fluorometer (Promega Corporation, Madison, USA) and 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA). Mean library size was 300 bp. Libraries were run in the rapid run flow cell and were single-end sequenced (65 bp) on HiSeq 1500 (Illumina, San Diego, CA 92122 USA).

ChIP-seq data analysis, CEs identification and CE annotation

The quality of raw fastq files were validated using FastQC and adapters were removed using Trimmomatic. Next, reads were mapped to the Gallus gallus 6.0 reference genome using Bowtie 2 with default parameters⁶² and PCR-duplicates were marked and removed using Picard (http://broadinstitute.github.io/picard/). The peaks were called using MACS2 with significance level threshold FDR <0.05, and normalization to input sample⁶³. Further, biological replicates were merged using BEDTools⁶⁴.

The unsupervised clustering of H3K27ac and H3K4me1 peaks was performed using DiffBind (https://bioconductor.org/packages/DiffBind/) with normalized IP samples to input. The experiment-specific lists containing anomalously enriched regions were generated using the GreyListChIP and further removed from datasets.

CEs were identified using in-house script. Briefly, the promoter regions (1 kb -/+ from TSS) were filtered out from H3K27ac and K3K4me1 dataset and nearby peaks (< 1 kb) were merged with GenomicRanges::reduce(min.gapwidth=1000)⁶⁵. Further, the conserved CEs were selected using BEDTools intersect followed by merging nearby genomic intervals with GenomicRanges::reduce(min.gapwidth=1000). The consensus and cell-specific CEs were identified using BEDTools::intersect -f 0.9 -r and BEDTools::intersect -v, respectively. The functional interpretation of CEs was performed using GREAT, with the genomic regions previously lifted to hg38 genome using liftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver)

The enrichment tracks for H3K4me1 and H3K27ac ChIP-seq data were generated using deepTools⁶⁶. Specifically, bamCoverage with reads per kilobase per million mapped reads (RPKM) per bin normalization was used. The enrichment tracks were visualized by loading to UCSC Genome Browser. The tracks visualized in UCSC Genome Browser were merged for biological replicates using transparent method (https://genome.ucsc.edu/cgi-bin/hgCollection).

Data availability

The RNAseq and ChIPseq data are available under the gene expression omnibus (GEO) accession number GSE198819.

Acknowledgments

We thank the members of the Center for Biomics-Genomics at the Erasmus University Medical Center for technical support. We also thank all the members of the Cell Biology Department for fruitful discussion. This work was supported by Polish National Science Centre (UMO-2015/19/B/NZ4/03184) and primary funding to the Department of Cell Biology at Erasmus University Medical Center.

Author contribution

Conceptualization: D.H., K.N., P.T., F.G.; Methodology: A.K., B.G., B.W., K.N., W.I.; Investigation in the lab: K.N.; Computational investigation: A.O., K.N., R.B., Software: A.O., K.N.; Resources: D.H and P.T; Data Curation: K.N, P.T; Writing Original draft: K.N; Writing, Review & editing: D.H., P.T.; Supervision: D.H and P.T; Project Administration: D.H and P.T; Funding Acquisition: D.H and P.T.

Khan, I. M.et al. The development of synovial joints. Curr Top Dev Biol 79, 1-36 (2007).
Spielmann, M.et al. Homeotic arm-to-leg transformation associated with genomic rearrangements at the PITX1 locus. Am J Hum Genet 91, 629-635 (2012).
Vissers, L. E.et al. Chondrodysplasia and abnormal joint development associated with mutations in IMPAD1, encoding the Golgi-resident nucleotide phosphatase, gPAPP. Am J Hum Genet 88, 608-615 (2011).
Girisha, K. M.et al. A novel multiple joint dislocation syndrome associated with a homozygous nonsense variant in the EXOC6B gene. Eur J Hum Genet 24, 1206-1210 (2016).
Boer, C. G.et al. Deciphering osteoarthritis genetics across 826,690 individuals from 9 populations. Cell 184, 4784-4818 e4717 (2021).
Chen, D.et al. Osteoarthritis: toward a comprehensive understanding of pathological mechanism. Bone Res 5, 16044 (2017).
Dieppe, P., Lim, K. & Lohmander, S. Who should have knee joint replacement surgery for osteoarthritis? Int J Rheum Dis 14, 175-180 (2011).
Roelofs, A. J., Rocke, J. P. & De Bari, C. Cell-based approaches to joint surface repair: a research perspective. Osteoarthritis Cartilage 21, 892-900 (2013).
Decker, R. S., Koyama, E. & Pacifici, M. Genesis and morphogenesis of limb synovial joints and articular cartilage. Matrix Biol 39, 5-10 (2014).
Spagnoli, A.et al. TGF-beta signaling is essential for joint morphogenesis. J Cell Biol 177, 1105-1117 (2007).
Ray, A., Singh, P. N., Sohaskey, M. L., Harland, R. M. & Bandyopadhyay, A. Precise spatial restriction of BMP signaling is essential for articular cartilage differentiation. Development 142, 1169-1179 (2015).
Shwartz, Y., Viukov, S., Krief, S. & Zelzer, E. Joint Development Involves a Continuous Influx of Gdf5-Positive Cells. Cell Rep 15, 2577-2587 (2016).
Luyten, F. P., Tylzanowski, P. & Lories, R. J. Wnt signaling and osteoarthritis. Bone 44, 522-527 (2009).
Jenner, F.et al. Differential gene expression of the intermediate and outer interzone layers of developing articular cartilage in murine embryos. Stem Cells Dev 23, 1883-1898 (2014).
Hyde, G., Dover, S., Aszodi, A., Wallis, G. A. & Boot-Handford, R. P. Lineage tracing using matrilin-1 gene expression reveals that articular chondrocytes exist as the joint interzone forms. Dev Biol 304, 825-833 (2007).
Gu, J.et al. Identification and characterization of the novel Col10a1 regulatory mechanism during chondrocyte hypertrophic differentiation. Cell Death Dis 5, e1469 (2014).
Pazin, D. E., Gamer, L. W., Capelo, L. P., Cox, K. A. & Rosen, V. Gene signature of the embryonic meniscus. J Orthop Res 32, 46-53 (2014).
Bian, Q.et al. A single cell transcriptional atlas of early synovial joint development. Development 147 (2020).
Fudenberg, G.et al. Formation of Chromosomal Domains by Loop Extrusion. Cell Rep 15, 2038-2049 (2016).
Garcia-Gonzalez, E., Escamilla-Del-Arenal, M., Arzate-Mejia, R. & Recillas-Targa, F. Chromatin remodeling effects on enhancer activity. Cell Mol Life Sci 73, 2897-2910 (2016).
Lupianez, D. G.et al. Disruptions of topological chromatin domains cause pathogenic rewiring of gene-enhancer interactions. Cell 161, 1012-1025 (2015).
Calo, E. & Wysocka, J. Modification of enhancer chromatin: what, how, and why? Mol Cell 49, 825-837 (2013).
Lohan, S.et al. Microduplications encompassing the Sonic hedgehog limb enhancer ZRS are associated with Haas-type polysyndactyly and Laurin-Sandrow syndrome. Clin Genet 86, 318-325 (2014).
Tayebi, N.et al. Deletions of exons with regulatory activity at the DYNC1I1 locus are associated with split-hand/split-foot malformation: array CGH screening of 134 unrelated families. Orphanet J Rare Dis 9, 108 (2014).
Bunyan, D. J.et al. Duplications upstream and downstream of SHOX identified as novel causes of Leri-Weill dyschondrosteosis or idiopathic short stature. Am J Med Genet A 170A, 949-957 (2016).
Fish, A., Chen, L. & Capra, J. A. Gene Regulatory Enhancers with Evolutionarily Conserved Activity Are More Pleiotropic than Those with Species-Specific Activity. Genome Biol Evol 9, 2615-2625 (2017).
McLean, C. Y.et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol 28, 495-501 (2010).
Cheung, K.et al. Histone ChIP-Seq identifies differential enhancer usage during chondrogenesis as critical for defining cell-type specificity. FASEB J 34, 5317-5331 (2020).
Chen, H.et al. Heads, Shoulders, Elbows, Knees, and Toes: Modular Gdf5 Enhancers Control Different Joints in the Vertebrate Skeleton. PLoS Genet 12, e1006454 (2016).
Hu, G., Codina, M. & Fisher, S. Multiple enhancers associated with ACAN suggest highly redundant transcriptional regulation in cartilage. Matrix Biol 31, 328-337 (2012).
Will, A. J.et al. Composition and dosage of a multipartite enhancer cluster control developmental expression of Ihh (Indian hedgehog). Nat Genet 49, 1539-1545 (2017).
Yao, B.et al. The SOX9 upstream region prone to chromosomal aberrations causing campomelic dysplasia contains multiple cartilage enhancers. Nucleic Acids Res 43, 5394-5408 (2015).
Visel, A., Minovitsky, S., Dubchak, I. & Pennacchio, L. A. VISTA Enhancer Browser--a database of tissue-specific human enhancers. Nucleic Acids Res 35, D88-92 (2007).
Buniello, A.et al. The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019. Nucleic Acids Res 47, D1005-D1012 (2019).
Shepherd, C.et al. Functional Characterization of the Osteoarthritis Genetic Risk Residing at ALDH1A2 Identifies rs12915901 as a Key Target Variant. Arthritis Rheumatol 70, 1577-1587 (2018).
Styrkarsdottir, U.et al. Severe osteoarthritis of the hand associates with common variants within the ALDH1A2 gene and with rare variants at 1p31. Nat Genet 46, 498-502 (2014).
Mokuda, S.et al. Wwp2 maintains cartilage homeostasis through regulation of Adamts5. Nat Commun 10, 2429 (2019).
Tachmazidou, I.et al. Identification of new therapeutic targets for osteoarthritis through genome-wide analyses of UK Biobank data. Nat Genet 51, 230-236 (2019).
Feregrino, C., Sacher, F., Parnas, O. & Tschopp, P. A single-cell transcriptomic atlas of the developing chicken limb. BMC Genomics 20, 401 (2019).
Keller, B.et al. Interaction of TGFbeta and BMP signaling pathways during chondrogenesis. PLoS One 6, e16421 (2011).
Kobayashi, T., Lyons, K. M., McMahon, A. P. & Kronenberg, H. M. BMP signaling stimulates cellular differentiation at multiple steps during cartilage development. Proc Natl Acad Sci U S A 102, 18023-18027 (2005).
Tsumaki, N.et al. Bone morphogenetic protein signals are required for cartilage formation and differently regulate joint development during skeletogenesis. J Bone Miner Res 17, 898-906 (2002).
Wang, W., Rigueur, D. & Lyons, K. M. TGFbeta signaling in cartilage development and maintenance. Birth Defects Res C Embryo Today 102, 37-51 (2014).
Brunet, L. J., McMahon, J. A., McMahon, A. P. & Harland, R. M. Noggin, cartilage morphogenesis, and joint formation in the mammalian skeleton. Science 280, 1455-1457 (1998).
Storm, E. E. & Kingsley, D. M. Joint patterning defects caused by single and double mutations in members of the bone morphogenetic protein (BMP) family. Development 122, 3969-3979 (1996).
Macias, D.et al. Role of BMP-2 and OP-1 (BMP-7) in programmed cell death and skeletogenesis during chick limb development. Development 124, 1109-1117 (1997).
Gamer, L. W.et al. The Role of Bmp2 in the Maturation and Maintenance of the Murine Knee Joint. J Bone Miner Res 33, 1708-1717 (2018).
Capellini, T. D.et al. Ancient selection for derived alleles at a GDF5 enhancer influencing human growth and osteoarthritis risk. Nat Genet 49, 1202-1210 (2017).
Egli, R. J.et al. Functional analysis of the osteoarthritis susceptibility-associated GDF5 regulatory polymorphism. Arthritis Rheum 60, 2055-2064 (2009).
Miyamoto, Y.et al. A functional polymorphism in the 5' UTR of GDF5 is associated with susceptibility to osteoarthritis. Nat Genet 39, 529-533 (2007).
Li, D.et al. Follistatin-like protein 1 is elevated in systemic autoimmune diseases and correlated with disease activity in patients with rheumatoid arthritis. Arthritis Res Ther 13, R17 (2011).
Hordyjewska, E., Jaruga, A., Kandzierski, G. & Tylzanowski, P. Novel Mutation of the RUNX2 Gene in Patients with Cleidocranial Dysplasia. Mol Syndromol 8, 253-260 (2017).
Kirkpatrick, T. J.et al. Identification of a mutation in the Indian Hedgehog (IHH) gene causing brachydactyly type A1 and evidence for a third locus. J Med Genet 40, 42-44 (2003).
Schlegel, W.et al. Dedifferentiation of human articular chondrocytes is associated with alterations in expression patterns of GDF-5 and its receptors. J Cell Mol Med 13, 3398-3404 (2009).
Sensiate, L. A.et al. Dact gene expression profiles suggest a role for this gene family in integrating Wnt and TGF-beta signaling pathways during chicken limb development. Dev Dyn 243, 428-439 (2014).
Chen, P., Carrington, J. L., Hammonds, R. G. & Reddi, A. H. Stimulation of chondrogenesis in limb bud mesoderm cells by recombinant human bone morphogenetic protein 2B (BMP-2B) and modulation by transforming growth factor beta 1 and beta 2. Exp Cell Res 195, 509-515 (1991).
Kawamura, I.et al. SnoN suppresses maturation of chondrocytes by mediating signal cross-talk between transforming growth factor-beta and bone morphogenetic protein pathways. J Biol Chem 287, 29101-29113 (2012).
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114-2120 (2014).
Dobin, A.et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 (2013).
Liao, Y., Smyth, G. K. & Shi, W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923-930 (2014).
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014).
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat Methods 9, 357-359 (2012).
Zhang, Y.et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol 9, R137 (2008).
Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841-842 (2010).
Lawrence, M.et al. Software for computing and annotating genomic ranges. PLoS Comput Biol 9, e1003118 (2013).
Ramirez, F.et al. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 44, W160-165 (2016).

No competing interests reported.

PDF

Editorial decision: Major revision
19 Apr, 2022
Reviews received at journal
13 Apr, 2022
Reviewers agreed at journal
31 Mar, 2022
Reviewers invited by journal
30 Mar, 2022
Editor assigned by journal
26 Mar, 2022
Editor invited by journal
18 Mar, 2022
Submission checks completed at journal
18 Mar, 2022
First submitted to journal
02 Mar, 2022

You are reading this latest preprint version

Identification of candidate regulatory elements controlling transcriptome during the formation of interphalangeal joints.

Status:

Version 1

Abstract

Figures

Introduction

Results

Microdissection of joint interzones and phalanges

An enhancer atlas of joint interzone and phalange identifies candidate REs involved in the regulation of cell identity

Integrative analysis of DEGs and CEs

CEs regulate skeletal malformation and disease-relevant genes, and are associated with a higher risk of OA

Discussion

Materials And Methods

Tissue collection

RNA extraction

RT-qPCR

RNA-seq library preparation

RNA-seq data analysis

Chromatin-immunoprecipitation (ChIP)

ChIP-sequencing

ChIP-seq data analysis, CEs identification and CE annotation

Declarations

References

Competing Interests

Supplementary Files

Status:

Version 1