SOX9 modulates the expression of key transcription factors required for heart valve development

One-third of cardiovascular birth defects are due to abnormal formation of the heart valves and valve insufficiency can lead to increased disease susceptibility later in life. Heart valves develop at two constrictions in the embryonic heart tube: between the atria and ventricles, known as the atrioventricular canal (AVC), and in the outflow tract (OFT), which forms the base of the aorta and pulmonary trunk (reviewed by de Vlaming et al., 2012). These regions form the mitral and tricuspid valves, and aortic and pulmonary valves, respectively (reviewed by Person et al., 2005). Starting at embryonic day (E) 9.5 in the mouse, epithelial-to-mesenchymal transition (EMT) occurs in the endocardial cells of the AVC, in part by responding to the activation of Notch signalling (Timmerman et al., 2004), and leads to their transformation into migratory cells. These cells invade and populate the extracellular matrix (ECM) that separates the endocardium and myocardium to form the cardiac cushions, the precursors of the valves and septa (reviewed by Person et al., 2005). These mesenchymal cells will proliferate, differentiate and remodel to eventually form the mature, thin, delicate valve leaflets (reviewed by de Vlaming et al., 2012).

Transcription factors (TFs) are essential for the precise coordination of lineage specification and differentiation, and ultimately the control of cell identity. In the heart valve, several TFs contribute to this role, but how these TFs are coordinated to regulate valve formation is not known. SOX (SRY-type box) 9, a high mobility group (HMG)-box TF, is highly expressed in the newly formed mesenchyme of the developing heart valves (Montero et al., 2002) and its expression is directly downstream of Notch signalling (Chang et al., 2014). In mouse, SOX9 has an essential role in heart valve development (Lincoln et al., 2007): loss of SOX9 in early valve formation leads to hypoplastic endocardial cushions, reduced proliferation and altered ECM deposition whereas at later stages of valve development, the loss of SOX9 causes abnormal ECM patterning, loss of cartilage-associated proteins, and thickened valves (Lincoln et al., 2007). Lack of SOX9 at either time point in heart valve development is embryonic lethal. In humans, haploinsufficient mutations in the SOX9 locus cause campomelic dysplasia (Wagner et al., 1994), in which there are skeletal abnormalities such as bowing of the long bones, sex reversal, and multiple organ defects, including pancreas and heart defects, that often lead to neonatal death (reviewed by Jain and Sen, 2014).

Embryonic heart valves have been shown to express a high number of genes associated with chondrogenesis, including extracellular matrix protein genes and a number of common key TFs, such as SOX9, TWIST1 and NFATC1 (Chakraborty et al., 2010a). Similar to heart valves, cartilage-containing tissues such as the limb differentiate from SOX9-positive mesenchyme, suggesting that their transcriptional programmes are analogous. During cartilage formation in developing limb buds, loss of SOX9 before mesenchymal condensation leads to cartilage agenesis and subsequent loss of bone formation (Akiyama et al., 2002). Limb buds with SOX9 deletion following the initial formation of cartilage have mesenchymal condensations that are severely hypoplastic and have defects in chondrocyte differentiation and proliferation (Akiyama et al., 2002). Loss of SOX9 in other organ systems also produces hypoplastic organs and defects in proliferation (Lincoln et al., 2007; Trowe et al., 2010; Rockich et al., 2013).

Despite the essential role of SOX9, its mechanisms of action and its transcriptional targets are not well understood. In adult hair follicle stem cells (HF-SCs), genes bound by SOX9 are required for the cells to maintain stemness via secreted factors in the niche (Kadaja et al., 2014). Recently, the global transcriptional targets of SOX9 in chondrocytes have begun to be explored (Liu and Lefebvre, 2015; Ohba et al., 2015). However, only a limited number of bona fide in vivo transcriptional targets of SOX9 have been identified in the embryonic limb (Bell et al., 1997; Bridgewater et al., 1998) and no direct in vivo SOX9 transcriptional targets are known in embryonic heart valves.

To explore the in vivo functional relevance of SOX9 in valve and limb development, we used chromatin immunoprecipitation coupled with next generation sequencing (ChIP-Seq) to generate genome-wide profiles of putative SOX9 DNA-binding sites in E12.5 mouse embryonic AVC and limb. We identified 2607 and 9092 SOX9 peak regions corresponding to 2453 and 5750 potential target genes in the E12.5 AVC and limb, respectively. SOX9 peaks common across multiple tissues supported a role for SOX9 in regulating cell proliferation and specifically identified AP-1 proteins and other cell cycle regulators as targets of SOX9 interactions. Moreover, in heart valves, we identified a group of essential TFs modulated by SOX9. By comparing these potential gene targets of SOX9 to transcriptional changes identified in RNA-Seq data from Sox9-deficient AVCs, we have demonstrated that SOX9 is a necessary upstream regulator required to modulate the expression levels of factors essential for heart valve formation.

At E12.5, SOX9 is widely expressed in mesenchyme throughout the developing valves (Akiyama et al., 2004) and in the condensing limb mesenchyme (Wheatley et al., 1996). To examine the regulatory role of SOX9 in these tissues, we generated genome-wide ChIP-Seq profiles of SOX9 DNA-binding sites for E12.5 heart valves (AVC) and E12.5 limbs. Subsequent to mapping and compiling of the sequencing reads, potential binding regions of SOX9, referred to as ‘peaks’, were identified using a false discovery rate (FDR) of 0.01 followed by statistical subtraction of an input DNA control. We identified 2607 and 9092 peaks in the AVC and limb, respectively (Fig. 1A; Tables S2-S4). Peaks of SOX9 interaction in the limb supported previously identified regulatory regions that are bound by SOX9, near to the genes Col2a1, Acan and Col11a2 (Bell et al., 1997; Bridgewater et al., 1998; Sekiya et al., 2000) (Fig. S1C). To identify similarities in the SOX9-initiated transcriptional programmes, we compared the SOX9 AVC and limb ChIP-Seq libraries and found 782 SOX9 peaks that were shared (30.0% of E12.5 AVC peaks) (Fig. 1A), confirming that SOX9 has common regulatory roles in valve and limb development. An additional 1825 and 8310 SOX9 peaks were unique in the E12.5 AVC and limb, respectively, indicating that SOX9 also has tissue-specific binding sites.

De novo motif analysis of the SOX9 ChIP-Seq peaks with SeqPos (Liu et al., 2011) generated a SOX monomer and a SOX dimer position weight matrix (PWM) (Fig. 1C) that is similar to the SOX9 JASPAR motif. In sex determination, SOX9 is known to bind as a monomer to a single DNA-binding motif in the regulatory regions of the loci for steroidogenic factor 1 (Nr5a1) and anti-Müllerian hormone (Amh), whereas genes required for chondrogenesis, such as Col10a1 and Col9a1, feature two SOX9-binding motifs separated by four nucleotides in reverse and complementary orientation to one another where SOX9 binds as a homodimer (Bernard et al., 2003). This is analogous to the SOX9 dimer motif found in our data. To determine how many SOX9 peaks contain a SOX monomer and/or dimer motif, we scanned all SOX9 peaks with Screen Motif (a tool written by Cliff Meyer and Len Taing in Cistrome; Liu et al., 2011) using the PWMs generated by SeqPos. We found that 77% of SOX9 limb peaks and 58% of SOX9 E12.5 AVC peaks contained at least one monomer or dimer motif (Table S5). In many cases, SOX9 peaks contain multiple SOX9 DNA-binding motifs per peak. The SOX dimer sequence was identified in 34% of SOX9 limb peaks compared with 13.5% of SOX9 AVC peaks, suggesting that SOX9 primarily binds as a monomer in E12.5 heart valves.

To determine where SOX9 binds in relation to genes, we associated all SOX9 ChIP-Seq peaks to putative target genes through a ‘yes-no’ process (Fig. S1D; supplementary materials and methods) that categorized where SOX9 peaks are relative to transcriptional start sites (TSSs), intragenic regions or intergenic regions. By weighting associations in this way, SOX9 peaks that localized to intragenic regions were associated with those genes rather than a potentially closer TSS. This mapping system allowed us to associate peaks to genes with fewer potential false positives and ensured that most peaks were only associated with one gene. Approximately 22% and 31% of SOX9 peaks in the limb and AVC, respectively, were located either directly over a TSS or in the 5 kb proximal promoter regions. This is in contrast to the SOX9 family members SOX6, which had only 13.5% of bound sites in TSS and 20 kb upstream regions (An et al., 2011), and SOX3, which had 11% of its bound sites in proximal promoter regions (McAninch and Thomas, 2014). Moreover, a remarkable 63.6% of SOX9 peaks shared between AVC and limb (497 out of 782) were located in the TSS/5 kb promoter regions (Fig. 1D). These data indicate that SOX9 frequently binds to promoters.

We identified 2453 and 5750 genes with associated SOX9 peaks in the E12.5 AVC and limb ChIP-Seq datasets, respectively (Tables S3, S4). Notably, 1605 gene loci that had associated SOX9 peaks were common in both AVC and limb. Of these, 782 SOX9 peaks were overlapping between the two libraries whereas the remaining gene loci had SOX9 peaks associated at different genomic locations in the AVC and limb (Fig. 1B). This suggests that SOX9 interacts with gene loci by using both shared and tissue-specific regulatory elements.

We further delineated common functions of SOX9 by comparing peaks in the AVC and limb with a publicly available SOX9 ChIP-Seq dataset that was generated in mouse HF-SCs (Kadaja et al., 2014) (Fig. 2A). We identified 293 genomic locations with SOX9 peaks in all three ChIP-Seq libraries (Fig. 2A; Fig. S2; Table S6), suggesting that SOX9 occupies analogous regions in mesenchymal tissues and stem cells. Gene Ontology (GO) analysis on the subset of genes associated with these shared SOX9-interacting regions, using the biofunction category (Ingenuity), featured proliferation of cells (Fig. 2B; Table S7). Indeed, 106 of the genes associated with common SOX9 peaks were included in the ‘proliferation of cells’ category (Table S8). Cell cycle regulators associated with SOX9 peaks in all three libraries include Junb, Cops5, Fosl1, Fosl2 and Fos. Additionally, several genes involved in cell proliferation had SOX9 peaks in heart valve and limb but not HF-SCs (Fig. 2C), such as Trp53 and Fgfr2. ChIP-qPCR on E12.5 AVC and limb further validated the SOX9 ChIP-Seq data, demonstrating that these shared regions were enriched for SOX9 occupancy when compared with IgG (Fig. 2C).

Given that SOX9 is required for the development of many different organs and is known to be involved in progenitor cell proliferation (reviewed by Pritchett et al., 2011), we anticipated that shared SOX9 peaks associated with proliferation genes would also be present in other SOX9-expressing tissues. ChIP-qPCR was performed on the E12.5 lung, gut and liver and demonstrated that many of the shared SOX9 peak regions were also enriched in these tissues (Fig. 2D). Together, our data suggest that analogous SOX9 DNA-binding regions are used in the regulation of cell proliferation across both different cell types and many developing tissues.

To identify tissue-specific activities of SOX9, we carried out GO analysis on genes associated with SOX9 non-overlapping peaks from each library and filtered out redundant GO terms between each tissue (Fig. 2B; Table S9). HF-SC-specific SOX9-associated genes were enriched for unique GO terms such as stem cell division, hair follicle morphogenesis and cell fate commitment (Fig. 2B). Limb-specific SOX9-associated genes revealed unique GO terms implicated in mesenchyme development, ECM organization and forelimb morphogenesis (Fig. 2B). SOX9 AVC-specific target genes identified unique GO categories involved in DNA binding, cardiac neural crest cell (NCC) development, and ascending aorta morphogenesis. Identifying genes involved in cardiac NCC in the AVC-specific GO categories supports the similarity between AVC and OFT cushion development. Overall, genes associated with tissue-specific SOX9 peaks strongly reflect the unique characteristics of each tissue and this suggests that SOX9 plays a role in, and is strongly linked with, tissue identity.

To parse out crucial transcriptional targets of SOX9 during heart valve development, we took advantage of a mouse model with specific deletion of Sox9 in which both endocardial cells and newly transformed mesenchymal cells of the developing AVC cushions fail to express SOX9 (Lincoln et al., 2007). These mice have been shown to die during embryogenesis and have heart valve defects including hypoplastic cardiac cushions, reduced mesenchyme proliferation and altered ECM composition. The Tie2-Cre mouse, used for this conditional system, has been previously shown to specifically express in endocardial cells and the resulting mesenchymal cells and not cardiomyocytes, epicardium or distal OFT mesenchyme (Kisanuki et al., 2001; Snarr et al., 2008). To ensure deletion of Sox9 in the E12.5 AVC, immunofluorescence for SOX9 was performed on Sox9^fl/fl (wild type; WT) and Sox9^fl/fl;Tie2-Cre (Sox9 cKO) embryonic hearts (Fig. S3A,B; Fig. S4). As expected, deletion of Sox9 leads to death in the majority of embryos by E13.5 due to severely hypoplastic and malformed AV valves (Fig. S3A,B; Fig. S4). Additionally, Sox9 cKO hearts had thinning of the ventricular walls (Fig. S3A, asterisk) and both the atrial and ventricular septums had not fused with the AVC in the Sox9 cKO compared with WT (Fig. S3A, arrowheads). Illustrating the specificity of Sox9 deletion with the Tie2-Cre lineage, SOX9 expression was lost only in the AVC whereas epicardial cells that descend from a different lineage still express SOX9 (Fig. S3A,B; Fig. S4). A decrease in proliferation in Sox9 cKO AVCs and AVC explant cultures was also noted (Fig. S5), as seen in previous work (Lincoln et al., 2007).

To identify downstream mRNA changes resulting from deletion of Sox9 that might be responsible for the valve defects, we compared the transcriptome of the E12.5 AVC in WT with that of Sox9 cKO. Because the efficiency of Sox9 deletion in the AVC by Tie2-Cre is variable, RNA was isolated from individual AVCs from E12.5 WT and Sox9 cKO, and qRT-PCR verified the loss of Sox9 transcript prior to pooling of two to three AVCs for RNA-Seq libraries. Duplicate libraries were created for each genotype (Fig. S6). Upon comparison, 634 genes were downregulated in the Sox9 cKO at least 1.5-fold and 610 genes were upregulated in the Sox9 cKO at least 1.5-fold (Tables S10, S11). Interestingly, downregulated genes in the Sox9 cKO AVC were categorized by GO as involved in cartilage development, mesenchyme differentiation and EMT (Fig. 3A). Upregulated genes in the Sox9 cKO included functions such as response to oxidative stress, gas transport and positive regulation of heart rate (Fig. 3B).

To further investigate genes that might be directly regulated by SOX9, we focused on the 145 genes that had both altered expression in AVCs lacking Sox9 and had an E12.5 AVC SOX9 peak associated. Of these 145 genes, for simplicity referred to as SOX9 target genes, ∼60% were downregulated and 40% were upregulated (Tables 1, 2; Table S12). Analysis of the up- and downregulated SOX9 target genes by Ingenuity Pathway Analysis highlighted enrichment in functions associated with transcription and cardiogenesis. Moreover, several of these genes have been linked to abnormal heart morphogenesis (Fig. 3C). Notable target TFs categorized as having roles in heart development that were downregulated in Sox9 cKO include Sox4, Hand2, Twist1, Foxp4, Mecom (also known as Evi1) and Pitx2; upregulated target genes also contained several TFs, such as Bhlhe40, Ddit3 and Junb (Tables 1, 2; Table S12). In addition, several ECM components, such as periostin and elastin, were reduced (Table 1; Table S12). Thus, our data suggest that SOX9 can act as a transcriptional activator and repressor. Moreover, our data support a direct role for SOX9 in regulating a network of TFs and ECM components that function during heart valve development.

Comparison of E12.5 AVC SOX9 ChIP-Seq and RNA-Seq data from Sox9 deficient valves suggests that SOX9 lies upstream of a network of TFs in developing heart valves. We confirmed reduced mRNA expression in Sox9 cKO AVC by qRT-PCR of the following TFs: Sox4, Mecom, Twist1, Pitx2, Hand2 and Nfia (Fig. 4A). In all cases, mRNA expression of these TFs was substantially reduced in Sox9 cKO AVCs. Two additional TFs associated with SOX9 peaks that are important in heart valve development, Lef1 and Tbx20, were also significantly reduced in Sox9 cKO AVC as shown by qRT-PCR (Fig. 4A) but were below the 1.5-fold change cut-off used for the RNA-Seq analysis (Tbx20=1.33, Lef1=1.29; Table S11). Of note, the two most downregulated (Btn1a1 and Prelp) and two of the most upregulated genes (Bhlhle40 and Fos) in the Sox9 cKO AVC also associated with SOX9 peaks were confirmed by qRT-PCR (Fig. 4B). Interestingly, little is known about the exact role of these factors during heart valve development.

To confirm that SOX9 can act as a transcriptional activator and repressor, SOX9 peak regions associated with Mecom, Nfia and Junb were cloned into luciferase reporter vectors and, together with SOX9 ectopic expression (pcDNA3-SOX9), luciferase activity was assayed in HEK293T cells (Fig. 4C). Mecom- and Nfia-associated SOX9 peaks were upstream of the genes in putative enhancers, whereas the SOX9 peak for Junb was located within the promoter. Mecom has been shown to be expressed in the developing heart valves (Hoyt et al., 1997; Bard-Chapeau et al., 2014) but its exact role in the AV valves is not well understood. NFI factors have been shown to work together with SOX9 to regulate downstream target genes in other systems (Nagy et al., 2011; Kang et al., 2012) and Junb has a well-established role in regulating cell cycle. SOX9 ectopic expression activated the Mecom and Nfia enhancers by 2.45- and 2.37-fold, respectively, whereas SOX9 inhibited the activity of the Junb promoter by 1.5-fold (Fig. 4C). These data support the proposal that SOX9 can function to both activate and repress transcription.

The putative enhancer for Mecom contains both a SOX9 dimer motif and several monomer motifs. To show that these SOX9 motifs are required for SOX9 regulation of this enhancer, we generated two mutant versions: one with a mutated SOX9 dimer motif (DM) and one with the dimer and monomer motifs mutated (Fig. S7). Deletion of the SOX9 dimer motif resulted in a substantial loss of SOX9-dependent luciferase activity from the Mecom enhancer, and mutation of all potential SOX9 motifs led to a complete loss of SOX9-dependent activity (Fig. 4D). These data suggest that the SOX9 dimer motif in the Mecom enhancer is essential for the majority of SOX9-dependent induction.

TWIST1 is highly expressed during the formation of the cardiac cushion mesenchyme and plays important roles in valve mesenchyme proliferation and differentiation (Shelton and Yutzey, 2008; Chakraborty et al., 2010b). Owing to the key role of TWIST1 in heart valve development, we focused on characterizing the transcript levels of Twist1 in the Sox9 cKO valves. qRT-PCR confirmed that the Twist1 transcript is not only reduced at E12.5, as discussed, but also reduced at an earlier time point (E10.5) in the Sox9 cKO AVC (Fig. S8). Using in situ hybridization, we found that Twist1 mRNA was specifically expressed in the valve mesenchyme in the E12.5 AVC (Fig. 4E; Fig. S9). Although precise quantitative analysis by in situ hybridization is not possible, Twist1 mRNA showed a similar level of reduction as the qRT-PCR in Sox9 cKO mutant valves at E12.5 (Fig. 4A). Despite the in vivo reduction of Twist1 expression in Sox9-deficient AVCs, the Twist1-associated enhancer is highly active without SOX9 overexpression using luciferase assays (data not shown). Motif analysis on the Twist1-associated enhancer identified a centrally located NF-Y motif and recent work suggests that NF-Y might recruit SOX9 (Shi et al., 2015).

Overall, the identification of SOX9 ChIP-Seq peaks in the E12.5 AVC, and alterations in mRNA levels between WT and Sox9 cKO AVCs support a model in which SOX9 functions to both activate and repress transcription during heart valve development. SOX9 plays an important role in regulating heart valve-specific regulatory networks by directly activating transcription factors such as Mecom and Nfia. However, SOX9 might also indirectly regulate factors like Twist1 via interactions with NF-Y or additional binding partners.

During heart valve development, SOX9 is highly expressed in the mesenchyme of the cardiac cushions of the AVC and OFT (Akiyama et al., 2004). Despite the important role of SOX9 in heart valve development, little is known about its mechanisms or the genes regulated by SOX9. Our study has suggested a dual role for SOX9 in the AVC by regulating proliferation of the mesenchyme cells and by modulating key transcription factors that function during heart valve development. To date, our study is the only one to examine SOX9 genomic interactions in embryonic heart and limb tissues.

In this study, we identified 2607 SOX9 ChIP-Seq peaks in the AVC of E12.5 embryonic heart. Owing to reported similarities in gene expression between valve formation and chondrogenesis (reviewed by Lincoln et al., 2006), and the key role of SOX9 in chondrogenesis (Akiyama et al., 2002), we generated and identified 9092 SOX9 ChIP-Seq peaks in E12.5 limb. Recent studies have similarly identified several thousand ChIP-Seq peaks for SOX9, including >27,000 in neonatal rib chrondrocytes (Ohba et al., 2015) and 8177 in HF-SCs (Kadaja et al., 2014). The relatively lower number of SOX9 peaks observed in the E12.5 AVC is potentially due to the limited amount of chromatin obtained from embryonic heart valves. Alternatively, it might reflect a more restricted role for SOX9 or less heterogeneity in the AVC at this time point. Remarkably, the overlap of SOX9 peaks between AVC, limb bud and HF-SCs is <2% of the total peaks, suggesting that SOX9-interacting regions are numerous and dynamic between tissues. Thus, these studies combine to illustrate the broad diversity of genomic regions that can interact with SOX9, either directly with the DNA or indirectly though other factors.

Although many SOX9 peaks do not overlap, we identified nearly 300 context-independent SOX9 peaks that are located at analogous genomic regions in the three different tissues examined. Many of these shared SOX9 peaks are associated with genes that function in proliferation, which supports a common function of SOX9 in maintaining a proliferative state during embryonic development (Fig. 5). It has long been known that SOX9 is linked with proliferation. However, a direct mechanistic connection and the transcriptional targets of SOX9 involved in proliferation have remained elusive. In AVC, limb and HF-SCs, SOX9 peaks are associated with three Fos transforming protein family members, Fos, Fosl1 and Fosl2. FOS and JUN family members heterodimerize to form the AP-1 complex, which is known to regulate cell proliferation and survival, in part via cyclin D1 expression (Shaulian, 2010). In mesenchymal stem cells, a stable SOX9 knockdown caused reduced proliferation, delayed S-phase progression, and increased cyclin D1 protein stability (Stockl et al., 2013). Of note, Junb also has an associated SOX9 peak. Although highly context dependent, JUNB is best known to inhibit cell growth by antagonizing JUN activity (Shaulian, 2010). We show that Fos and Junb expression is upregulated in SOX9-deficient AVCs, suggesting that their increased expression could contribute to tissue hypoplasia.

In addition to AP-1 complex factors, several genes with roles in cell proliferation have SOX9 peaks in their promoters or potential regulatory regions, including Cops5, Srpk2, Akt2, Eed, Hdac1, Hdac2, p53 (Trp53) and Prkaca (PKA). COPS5 (COP9 signalosome subunit 5) associates with JUN proteins to increase binding specificity and can degrade the cell cycle inhibitor p27Kip1 (Cdkn1b) (Claret et al., 1996). Loss of Cops5 in embryonic limb results in shortened limbs due to impaired chondrogenesis and Sox9 levels are decreased in mutant long bones (Bashur et al., 2014) suggesting a potential feedback loop between SOX9 and COPS5. SRPK2 (SRSF protein kinase 2) can promote proliferation and cell cycle progression by enhancing cyclin D1 levels (Jang et al., 2009), and AKT2 regulates progression of cell cycle via phosphorylation of its targets including cyclin-dependent kinase inhibitors and maintaining protein stability of MYC and D-type cyclins via GSK3β (Xu et al., 2012). EED (embryonic ectoderm development), HDAC1 and HDAC2 (histone deacetylase 1 and 2) are epigenetic regulators associated with cell proliferation (Bracken et al., 2003; Kelly and Cowley, 2013). p53 activates DNA repair and arrests proliferation by pausing the cell cycle. If the DNA damage is severe, p53 initiates cell death (Giono and Manfredi, 2006). PKA (protein kinase A) is induced by cyclic adenosine monophosphate and regulates cellular growth and proliferation through a variety of mechanisms (Stork and Schmitt, 2002). Of note, PKA phosphorylates SOX9 and increases its activity during chondrogenesis (Huang et al., 2000). The activities of p53 and PKA in cell proliferation might be specific to mesenchyme as SOX9 peaks are found in AVC and limb but not HF-SCs.

Taken together, our data suggest a model in which SOX9 promotes proliferation across multiple cell types during development, including heart valves, by interacting with promoters or enhancers of genes encoding proliferative factors such as AP-1 proteins, kinases, and histone modifiers. As many of these genes are misregulated in a diverse range of diseases, identification of the nature of SOX9 interactions with the regulatory regions of these genes could help to elucidate how regulation of these genes is controlled.

Previous work demonstrated that Sox9 deletion using Tie2-Cre is embryonic lethal with hypoplastic cardiac cushions, decreased mesenchyme proliferation and altered ECM (Lincoln et al., 2007). Despite this crucial role of SOX9 in valve development, a very limited number of genes show drastic changes in gene expression in the Sox9 cKO AVC. Surprisingly, <6% of genes with associated SOX9 peaks showed altered mRNA expression. Approximately 60% of the SOX9 target genes were downregulated. Although in vitro promoter analysis supported both a positive and negative role of SOX9 in gene regulation, the induction or repression of gene expression by SOX9 was at most threefold. Together, these data suggest that SOX9 might function as a modulator of gene expression rather than being required to produce large changes in gene expression. SOX9 has recently been suggested to act as a pioneer factor and has been linked with super enhancers (Adam et al., 2015), implying that SOX9 plays a role in chromatin dynamics. A recent model for SOX9 function in chondrocytes proposes two types of SOX9 interactions: Class I SOX9 engagement regions are TSS biased and indirectly bind SOX9 via the transcriptional machinery; Class II regions are more distally located, are directly bound by SOX9 and are enriched for skeletal enhancers (Ohba et al., 2015). Of note, most of the SOX9 peaks associated with genes that have altered gene expression in the AVC are located outside of the promoter (data not shown).

Several of these SOX9 target genes are known to be important for valve formation, including ECM-related genes and TFs. Downregulated TFs with associated SOX9 peaks included Lef1, Pitx2 and Hand2. Loss of Lef1 (via TBX20 deletion), Pitx2 or Hand2 is associated with valve defects (Liu et al., 2002; Holler et al., 2010; Cai et al., 2013). TFs with the most reduced expression in Sox9 cKO valves were Twist1, Sox4 and Mecom. Similar to SOX9, these three TFs are highly expressed in the cardiac cushions and mutations in each of these factors cause major valve defects, resulting in embryonic fatality (Ya et al., 1998; Vincentz et al., 2008; Bard-Chapeau et al., 2014).

SOX proteins are known to regulate TFs that will function as their future co-factors (Kamachi and Kondoh, 2013). SOX9 is known to regulate and cooperate with SOX5/6 to regulate target genes in the developing limb (Han and Lefebvre, 2008; Liu and Lefebvre, 2015) and both have associated SOX9 peaks in our limb ChIP-Seq data. Similarly, SOX9 might activate SOX4 in the heart to help co-regulate valve-specific genes. Motif analysis revealed EVI1 (protein product derived from Mecom) as another potential co-factor for SOX9 and comparison of EVI1 peaks in cancer cells (Bard-Chapeau et al., 2012) with SOX9 AVC peaks identified hundreds of overlapping target genes (P.A.H., V.C.G. and R.C., unpublished).

We have shown that SOX9 can modulate the levels of Twist1 expression in the developing AVC. In the absence of SOX9, Twist1 mRNA expression was reduced by approximately threefold in the valve mesenchyme. TWIST1 can induce proliferation and migration of valve mesenchyme during early valve formation (Shelton and Yutzey, 2008; Chakraborty et al., 2010b) and, following EMT, TWIST1 plays a role in regulating differentiation of the AVC mesenchyme (Vrljicak et al., 2012). When TWIST1 persists at later stages of valve development, it leads to increased mesenchyme proliferation, increased TBX20 expression, and primitive ECM (Chakraborty et al., 2010b). TWIST1 directly regulates Tbx20 (Lee and Yutzey, 2011), but, of note, Twist1-null hearts do not show a difference in Tbx20 levels in the AVC compared with WT (Vincentz et al., 2008). Interestingly, we found that TBX20 was downregulated in the Sox9 cKO AVC and that SOX9 occupied regulatory regions associated with TBX20. These data suggest that TWIST1 and SOX9 might cooperate to regulate Tbx20 in developing heart valves. Furthermore, SOX9 is downregulated in Twist1-null AVC (Vrljicak et al., 2012), suggesting the existence of a feedback loop between the two factors.

Taken together, our data reveal that SOX9 occupies regulatory regions and impacts the expression of key TFs that are vitally important for heart valve development (Fig. 5). Intriguingly, many of these TFs have been suggested to regulate each other's expression. For example, TWIST1 has been shown to regulate Tbx20 (Lee and Yutzey, 2011) and TBX20 in turn has been shown to regulate LEF1 in the valve endocardial cells (Cai et al., 2013). EVI1 has been shown to regulate SOX4 and alters its expression (Bard-Chapeau et al., 2014), and EVI1 and SOX4 can collaborate together in myeloid leukemia (Boyd et al., 2006). This demonstrates that complex interactions occur at multiple levels in transcriptional regulation by SOX9 and suggests that the essential TFs are regulated in numerous ways to ensure proper valve formation. Our work suggests that SOX9 and its transcriptional target TFs form a gene regulatory network to drive valve morphogenesis. In humans, aberrant expression of SOX9 and its transcriptional targets have been associated with adult heart valve disease (Hulin et al., 2013). Thus, understanding the SOX9-initiated transcription networks in heart valve development could provide additional insights into adult heart valve disease.

The Animal Care Committee at the University of British Columbia approved all animal procedures. C57BL/6J and ICR mice were used for ChIP-Seq and ChIP-qPCR validation, respectively. Sox9^fl/fl (B6.129S7-Sox9tm2Crm/J) and Tie2-Cre [B6.Cg-Tg(Tek-cre)12Flv/J] mice (Jackson Laboratories) were bred as described (Lincoln et al., 2007) and genotyped (Table S1). See supplementary materials and methods for additional information.

Three independent ChIPs generated from embryonic tissues (AVC ChIPs from 45, 45 and 74 embryos; limb ChIPs from 12, 11 and 16 limbs) were pooled for SOX9 ChIP-Seq. See supplementary materials and methods for details. ChIPs were performed as described (Vrljicak et al., 2012) with modifications (see supplementary materials and methods). Protein A/G beads (Pierce) with 3 μg of rabbit polyclonal anti-mouse SOX9 antibody (Millipore AB5535) or 3 μg of rabbit polyclonal IgG antibody (Santa Cruz sc-2027) were used. Unbound sonicated DNA was sequenced as input. ChIP-Seq libraries were constructed and sequenced at Canada's Michael Smith Genome Sciences Centre (Vancouver, BC). The Burrows-Wheeler Aligner (Li and Durbin, 2009) aligned reads to mm9. FindPeaks3.1 (Fejes et al., 2008) identified regions of enrichment or ‘peaks’. An FDR was applied to threshold ChIP-Seq data (Robertson et al., 2008) and false positives were limited using input. A local z-score was calculated between peak height and control coverage and peaks below the threshold were filtered out (Fig. S1A,B). Peaks that passed filtering, z-score, peak height, and FDR-based peak height cut-offs were retained for analysis. Wig/BED files were analyzed using UCSC Genome browser (Kent et al., 2002), Galaxy (Blankenberg et al., 2010) and Cistrome (Liu et al., 2011). Peaks were associated with genes through a ‘yes-no’ process as depicted by the flowchart (Fig. S1D) to reduce the number of mis-associated peaks (see supplementary materials and methods for details of peak-to-gene associations). Data have been deposited in Gene Expression Omnibus under accession number GSE73225.

Individual AVCs were dissected and put into Trizol (Thermo Fisher Scientific) and genotyped from unused embryonic tissue. Sox9 cKO hearts were only taken for analysis if the heart was functional. cDNA was synthesized to confirm loss of Sox9 in the AVCs by qRT-PCR. Two to three AVC RNA samples for each genotype were pooled. Duplicate RNA-Seq libraries were generated and sequenced on Illumina Mi-Seq or Next-Seq for each genotype. Reads were aligned with Tophat2 (Kim et al., 2013) to mm9. Fragments per kilobase of exon per million reads (FPKMs) were calculated using Cufflinks (Trapnell et al., 2010) and gene FPKMs were an average of duplicate libraries for each genotype. EdgeR (Robinson et al., 2010) was used for normalization and differential expression was determined by log2 fold change between WT and Sox9 cKO gene values. For downstream analyses, genes had to be ≥0.5 FPKM and ≥±0.58 log2 fold change. GO analysis was performed using GOrilla (Eden et al., 2009) and Ingenuity Pathway Analysis (Qiagen). For further details, see supplementary materials and methods.

cDNA was synthesized with Transcriptor First Strand cDNA Synthesis Kit (Roche). HiFi Taq polymerase (Thermo Fisher Scientific) was used for genotyping. Genomic DNA was isolated from embryonic tissue using KAPA express extract. ChIP-qPCR and qRT-PCR were performed with FastStart Universal SYBR Master (Roche) on the ABI 7900HT Fast Real-Time PCR System. See Table S1 for primers. Taqman assays (Thermo Fisher Scientific) used Perfecta qPCR Fastmix (Quanta Biosciences). Actb and Gapdh were used for relative quantification for qRT-PCR and Taqman assays, respectively. ChIP-qPCR fold enrichment was calculated by 2ΔCt difference between IgG ChIP and SOX9 ChIP).

Hearts were fixed in 4% paraformaldehyde (Sigma) overnight, subjected to a sucrose gradient, embedded in OCT (Sakura) and cryosectioned at 8 μm thickness. IF was performed as described (Chang et al., 2011). The primary antibody was rabbit anti-mouse SOX9 (Millipore AB5535; 1:600) and the secondary antibody was anti-rabbit Alexa Fluor 488 or 594 (Thermo Fisher Scientific). Rabbit anti-phospho histone H3 (Abcam) was used for proliferation assays. 4′,6-Diamidino-2-phenylindole, dihydrochloride (Sigma) was used to stain nuclei. Images were captured on a Zeiss Axioplan 2 microscope or TCS SP5 Leica confocal microscope. In situ hybridization was performed as described (Hou et al., 2007). Images were captured on a Zeiss Axioplan 2 microscope.

HEK293T cells were maintained in DMEM (Stemcell) with 10% foetal bovine serum (Thermo Fisher Scientific) and transfected at 60% confluency using polyethylenimine (Polysciences) (Baker et al., 1997). Luciferase assays were performed 2 days post-transfection using the Dual-Luciferase Reporter System (Promega). Regions containing SOX9 peaks were PCR amplified (Table S1) and cloned into either a modified pGL3 promoter luciferase plasmid (Promega), containing an E1b promoter (pGL3-E1Bp; Benchabane and Wrana, 2003), for enhancers) or pGL4-basic (Promega) promoter-less vector (for promoters). Luciferase activity was tested in the presence of a pcDNA3-SOX9 expression vector (Lefebvre et al., 1997) (0.1 µg/well, 24-well plate) or pcDNA3 backbone (0.1 µg/well). Firefly luciferase activity was normalized to Renilla luciferase activity for each sample. Enhancer and promoter firefly luciferase activity is shown to the empty vector. Mecom enhancer mutants were generated using gBlocks by Integrated DNA Technologies (IDT) (Fig. S7).

We would like to thank Dr Aly Karsan for all of his excellent ideas and suggestions throughout this project. We would also like to thank Dr Amanda Kotzer for all of her help in managing this project.