Hoxc8 is a homeodomain transcription factor that regulates pattern formation, cell migration, and differentiation (Lei et al., 2005; Pearson et al., 2005). Loss- or gain-of-function studies have suggested that Hoxc8 is essential for skeletal pattern formation, hematopoiesis, and cartilage differentiation during embryogenesis (Kruger and Kappen, 2010; Le Mouellic et al,. 1992; Shimamoto et al., 1999; Tiret et al., 1998; Yueh et al., 1998). To further elucidate the underlying mechanisms of developmental defects caused by Hoxc8 mutation, it is necessary to identify the transcriptional target genes of Hoxc8. To date, several target genes have been identified, and several preliminary Hoxc8-mediated regulatory networks have been proposed. Hoxc8 directly downregulates Opn, Zac1, Ncam, and Pedf, and upregulates Cdh11 in mouse embryonic fibroblast cells (Lei et al., 2005; 2006). We also previously identified Pcna, ZNF804A, and Mgl1 as Hoxc8 direct target genes in mouse embryonic fibroblast cells (Chung et al., 2010; Min et al., 2010; Ruthala et al., 2011). In another study, although 15 of 21 target genes were upregulated by Hox genes (Pearson et al., 2005), more than half of the known Hoxc8 target genes were downregulated by Hoxc8. This bias in Hoxc8 regulation patterns might result from use of a single cell type, mouse embryonic fibroblast cells at specific developmental stage in vitro. To more accurately reflect in vivo Hoxc8 transcriptional specificity, more target genes need to be identified in diverse cell types and tissues throughout the stages of embryonic development.

Several characteristics of Hoxc8 make target gene identification difficult. First, it has low DNA binding sequence specificity. Hoxc8 proteins bind as monomers or multimers to specific sequence motifs (TAAT/ATTA, TTAT, and ATAA) in their target genes (Lei et al., 2006). However, other Hox proteins can potentially bind to the same elements, and the DNA-binding specificities are modified through interactions with cofactors, such as Pbx and Meis (Ladam and Sagerstr?m, 2014, Mann, 1995; Mann and Chan, 1996; Mann and Affolter, 1998; Moens and Selleri, 2006;). Second, little is known about Hoxc8 cofactors. Interaction of HOXC8 with PBX1 has been reported in prostate cancer cells (Kikugawa et al., 2006) but in no other tissues. Interestingly, Hoxc8 interacts with Smad1 during osteoblast differentiation (Hullinger et al., 2001, Shi et al., 1999; Yang et al., 2000). This indicates that other unknown factors, even those without a homeodomain like Smads, might interact with Hoxc8 in vivo. Third, Hox members act as both transcriptional activators and repressors, depending on their cellular context (Pearson et al., 2005). Therefore, careful experimental designs, with attention to cell types, developmental stage, time, and location, are necessary to validate candidate genes.

In an attempt to identify more Hoxc8 target genes, we analyzed eight archived microarray datasets generated from diverse cell types and tissues throughout mouse embryonic development. We found that Hoxc8 is often coexpressed with several transcription factors, including known cofactors Pbx1 and Meis, and genes related to extracellular matrix (ECM) or organization. We then randomly selected five genes and assessed their likelihood of being Hoxc8 target genes.

Gene expression data analysis

To identify genes with expression patterns similar to those of Hoxc8 during mouse embryonic development, we used gene expression profiles from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). We analyzed eight datasets (Table 1) and assessed the similarity of each gene expression pattern to that of Hoxc8 based on Pearson’s correlation coefficient (r-value > 0.75) and t-test (p-value < 0.001). For the t-test, we grouped samples based on differential Hoxc8 expression, as shown in Table 1. We selected genes with r-values greater than 0.75 in at least four out of the eight datasets. Genes with Pearson’s correlation coefficients greater than 0.75 in at least four out of the eight datasets were subjected to gene set analysis with the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (http://david.abcc.ncifcrf.gov/home.jsp) program. Gene Ontology (GO) annotations with an enrichment score > 1.0, p-value < 0.01, and false discovery rate (FDR) ≤ 5% were considered biologically significant. We used the rVista sequence analyzer (http://zpicture.dcode.org/) to search for Hox consensus binding elements present in the 5 kb promoter, first exon, and first intron of putative target genes. The gene network was analyzed with Ingenuity Pathway Analysis (IPA) (http://www.ingenuity.com/products/ipa).

Animal preparation

To obtain E14.5 embryos, male and female ICR mice were caged together for mating at around 6 pm. The next morning, when vaginal plugs were present, was defined as 0.5 days post-coitum (dpc) or as an E0.5 embryo. After 14 days, the pregnant female mice were sacrificed, and the E14.5 embryos were extracted. The maternal and extra-embryonic tissues, cervical region, internal organs, tail bud, and limbs were removed. The embryos were divided into three parts: brain, trunk anterior (somites 12?23), and trunk posterior (somites 24?41). Total RNA was isolated from each part. This study was carried out in strict accordance with the recommendations in the Guide for the Institutional Animal Care and Use Committee of Yonsei University College of Medicine. The protocol for obtaining embryonic samples was approved by the Committee on Animal Research at Yonsei University College of Medicine (permit number 2013-0174-1).

Cell culture and transfection

Three cell lines, MC3T3-E1, NIH3T3, and C3H10T1/2, were cultured in Dulbecco’s modified Eagles medium (WelGENE Inc., Korea) supplemented with 10% fetal bovine serum (FBS; WelGENE Inc., Korea) and 100 μg/ml penicillin-streptomycin (WelGENE Inc., Korea) at 37°C inside a 5% CO₂ and 95% humidified incubator. The construction of the pcDNA3.1-Hoxc8 plasmid, harboring the murine Hoxc8 gene, has been previously described (Kwon et al., 2003). pcDNA3.1-Hoxc8 or pcDNA3.1 empty vector was transfected into MC3T3-E1 cells using Lipofectamine 2000™ reagent (Invitrogen, USA), as indicated by the manufacturer. To establish stable cell lines expressing Hoxc8, NIH3T3 and C3H10T1/2 cells were transfected with Lipofectamine 2000™ reagent and pcDNA3.1-Hoxc8 or control empty vector and selected in culture media containing 500 μg/ml G418 antibiotic (Invitrogen, USA). The media was changed every three days. Cells were subcultured when the cells reached 90% confluence.

MC3T3-E1 cells were seeded into 12-well dishes at 1 × 10⁵ cells/well. After an overnight incubation, the media was replaced with fresh media supplemented with 10 ng/ml of TGF-β2 (R&D Systems, Inc., USA). Cells were harvested at time points 0 h, 6 h, and 12 h after treatment with TGF-β2, and total RNA was isolated using Trizol reagent (Invitrogen, USA).

Knockdown of Hoxc8 with short interfering RNA (siRNA)

The sequences of the Hoxc8 and control siRNAs are as follows: Hoxc8 sense, 5′-AGA CGC CUC CAA AUU CUA UTT-3′, Hoxc8 antisense, 5′-AUA GAA UUU GGA GGC GUC UTT-3′, control sense, 5′-AUG AAC GUG AAU UGC UCA ATT-3′, and control antisense, 5′-UUG AGC AAU UCA CGU UCA UTT-3′ (Samchully Pharm Co., Ltd., Korea). MC3T3-E1 cells were seeded into 12-well dishes (1 × 10⁵ cells/well) and incubated overnight. Then, 100 nM siRNA (final concentration) was transfected into the cells using HiPerfect transfection reagent (Qiagen, Germany), according to the manufacturer’s instructions. Cells were harvested 60 h after transfection, followed by isolation of total RNA. All experiments were performed in triplicate and representative examples are shown.

RNA isolation and semi-quantitative RT-PCR

Total RNA was isolated using Trizol reagent, according to the manufacturer’s instructions (Invitrogen, USA). Two micrograms (μg) of total RNA was reverse-transcribed with ImProm-II™ Reverse Transcriptase (Promega, USA) and poly (dT)₂₀, according to the manufacturer’s instructions. Semi-quantitative RT-PCR was performed with hTaq DNA polymerase (Solgent, Korea) using the following thermo cycling conditions: initial denaturation for 5 min at 95°C, followed by 30?33 cycles of 94°C for 30 s (denaturation), 58°C for 30 s (annealing), and 72°C for 30 s (polymerization). We determined the highest PCR cycle numbers at which PCR products increased linearly and were detectable on agarose gel. Primer sequences are listed in Supplementary Table S1. The PCR products were analyzed on a 1.5% agarose gel containing ethidium bromide.

Hoxc8 was differentially expressed in several gene expression datasets

To select archived gene expression datasets in which Hoxc8 was differentially expressed, the GEO database was analyzed. We selected datasets by focusing on two criteria: datasets associated with studies on mouse development and cell differentiation, as well as those with greater than two-fold differential Hoxc8 expression under the given experimental conditions, tissue types, or cell types. Eight datasets (GDS2843, GDS2743, GDS1500, GDS2123, GDS2209, GDS2699, GDS2044, and GDS2421) met our criteria (Table 1). For example, Hoxc8 was differentially expressed during mammary gland development (GDS2843; Fig. 1). In addition, when different cell types were compared, Hoxc8 was overexpressed in white adipocytes compared to brown adipocytes (GDS2743), spinal cord compared to dorsal root ganglion (GDS2209), mesenchymal cells compared to epithelial cells (GDS2699), and in skin fibroblasts among fetal corneal, skin, and tendon fibroblasts (GDS1500). Next, we examined the datasets in which the known Hoxc8 target genes Opn and Fzd2 were concordantly expressed with Hoxc8 (Table 1). Representative data from the GDS2843 dataset, which originated from the developing mammary gland, are shown in Fig. 1. Therein, the gene expression levels of Opn and Fzd2 were negatively and positively correlated with that of Hoxc8, respectively. In most datasets, the expression patterns of Fzd2 were highly correlated with Hoxc8 (Table 1) except in the GDS1500 and GDS2209 datasets.

Enriched functions of Hoxc8 coexpressed genes

To search for genes coexpressed with Hoxc8 in each dataset, we assessed the similarity of gene expression patterns between each gene and Hoxc8 using Pearson’s correlation coefficient (r-value > 0.75) and t-test (p-value < 0.001). A total of 567 genes that had r-values greater than 0.75 in at least four of eight datasets were selected (Fig. 2). Among these, 23 genes were found to have r-values greater than 0.75 in six of the eight datasets (Table 2). Among these 23 genes, Ncam1 and Fzd2 are already known to be Hoxc8 target genes, and Meis1 is a known Hox co-factor. As the most of the known target genes are associated with ECM or cell adhesion, many of the identified genes were found to be associated with ECM (Isir, immunoglobulin superfamily containing leucine-rich repeat; Timp3; Sparc; Fbln1; Tgfbi, transforming growth factor beta-induced, 68kDa) or cell adhesion (Thbs2, thrombospondin 2; Fyn; Mylk, myosin light chain kinase).

To gain insight into the functional characteristics of Hoxc8 coexpressed genes, we performed gene set analysis for 567 genes, using the DAVID program. Genes associated with the ECM were most significantly enriched (Fig. 2 and Supplementary Table S2). In addition to the 23 genes mentioned above, collagen family proteins (Col27a1, Col3a1, Col1a2, Col1a1, Lox, Col11a1, Col5a2, Col5a1, and Col4a5) and proteolysis metallopeptidases (Adam19, Adam9, and Adamts4, -5, and -12) were highly enriched. These genes are closely connected with the Hoxc8 subnetwork (Supplementary Fig. S1). Genes associated with the development of the skeletal system were also highly enriched. Of note, these included six other Hox genes: Hoxa5, Hoxa7, Hoxc6, Hoxc9, Hoxd8, and Hoxd9. These genes are neighboring genes localized next to Hoxc8 and its paralogues located in other clusters, with the exception of the Hoxb cluster. We also performed gene set analysis for genes that were negatively correlated with Hoxc8 at least four times in eight datasets (Supplementary Table S3). Most significantly, the enriched gene set was associated with the generation of precursor metabolites and energy, especially the TCA cycle and localization in the mitochondria.

Interestingly, 17.4% of the coexpressed genes were transcriptional regulators, such as transcription factors involved in a pattern formation, Hox cofactors, and chromatin remodeling factors. In addition to Meis1, other well-known Hox cofactors, including Meis2, Pbx1, and Pbx3, belonging to the TALE (Three Amino acid Loop Extension) family, were concordantly expressed in four datasets. Zhx1, Prrx1, Zeb2, Pknox1, and Meox2 also have a homeobox domain, as do members of the TALE family. Zhx1 and Prrx1 are members of the zinc-finger and homeobox protein families, respectively, both of which act as repressors. However, whether they function as Hox cofactors is unknown. Although the tissues or cells used in the gene expression datasets are heterogeneous, gene set analysis of the coexpressed genes seemed to summarize Hoxc8 phenotypes.

Selected genes were endogenously coexpressed with Hoxc8

To assess the possibility that the 567 selected genes are coexpressed with Hoxc8, we selected five genes and compared their expression patterns with that of Hoxc8 in mouse embryos and cell lines. We selected the five following genes showing different coexpression frequencies against Hoxc8: Tgfbi (transforming growth factor, beta induced), Ptpn13 (protein tyrosine phosphatase non-receptor type 13), Prkd1 (protein kinase D1), Adam19 (a disintegrin and metallopeptidase domain 19), and Aldh1a3 (aldehyde dehydrogenase 1a3). Tgfbi was the most frequently coexpressed (6 of 8 datasets) gene, while Aldh1a3 was coexpressed in only three out of eight datasets. Prkd1 and Ptpn13 were coexpressed in five, and Adam19 was coexpressed in four of the eight datasets (Supplementary Table S4). These proteins have been implicated in various biological systems, such as ECM and adhesion (Tgfbi and Adam19), phosphate metabolic processes (Ptpn13 and Prkd1), and retinoic acid biosynthetic processes (Aldh1a3).

First, the endogenous expression patterns of these five genes were analyzed in mouse embryos (Fig. 3). E14.5 embryos were dissected into three parts, brain, trunk anterior, and trunk posterior tissues, and then semi-quantitative RT-PCR was performed after isolating total RNA from each tissue. In agreement with previous reports (Kwon et al., 2005; Min et al., 2012; 2013), Hoxc8 was strongly expressed in the trunk anterior and moderately expressed in the trunk posterior; it was not expressed in the brain (Fig. 3). All five of the examined genes, along with Fzd2 as a control, exhibited similar expression patterns to those of Hoxc8, although Prkd1 and Ptpn13 were also weakly expressed in brain tissue (Fig. 3).

Since Hoxc8 has been reported to be induced in preosteoblastic cells through the TGF-β signaling pathway (Li et al., 2006; Yang et al., 2000), we tested whether the putative target genes are also concordantly induced along with Hoxc8 in the MC3T3-E1 preosteoblastic cell line. After treating the cells with TGF-β2, they were harvested at multiple time points (0 h, 6 h, and 12 h), RNA was extracted, and semi-quantitative RT-PCR was performed. As expected, expression of Hoxc8 gradually increased following TGF-β2 treatment (Fig. 4). Similarly, four genes of interest, Adam19, Ptpn13, Prkd1, and Tgfbi, were also induced (Fig. 4), although Aldh1a3 was not detected at any time points (data not shown), probably due to gene silencing in MC3T3-E1 cells. Expression of the selected genes, except Aldh1a3, was confirmed to be positively correlated with Hoxc8 in mouse embryonic tissues and cell lines, in which the TGF-β signaling pathway was activated.

Hoxc8 positively regulates expression of the selected genes in vitro

To further confirm the correlation between Hoxc8 and expression of the selected genes, overexpression and knockdown conditions of Hoxc8 were generated in MC3T3-E1 preosteoblastic cells using a Hoxc8 expression vector (pCDNA3.1-Hoxc8) and siRNA against Hoxc8, respectively. When Hoxc8 was transiently overexpressed through introduction of the Hoxc8 expression vector into MC3T3-E1 cells (Fig. 5A), the expressions of Adam19, Ptpn13, Prkd1, and Tgfbi were concordantly increased. Likewise, the expression levels of the four genes decreased when Hoxc8 was knocked down (Fig. 5A). We also stably transfected the Hoxc8 expression vector into C3H10T1/2 and NIH3T3 cell lines, which lack any endogenous expression of Hoxc8. The basal expression levels of Adam19, Ptpn13, Prkd1, and Tgfbi were increased in the Hoxc8-expressing cell lines (Fig. 5B). Exceptionally, Aldh1a3 was downregulated upon overexpression of Hoxc8 in C3H10T1/2 cells.

Microarray technology, in combination with chromatin immunoprecipitation (ChIP), has greatly enhanced the discovery of transcription factor target genes (Lei et al., 2005; 2006). Nonetheless, identification of Hox target genes, including Hoxc8, has yet to reach a point that would facilitate delineating their individual roles in specifying the identity of body segments, cell differentiation, migration, and proliferation (Hueber and Lohmann, 2008). Low DNA binding sequence specificity, context-dependent activation or repression, and unknown cofactors present extensive challenges for identification of Hox target genes. These potential pitfalls spurred us to undertake alternative methods for identifying Hox target genes. Genes that show similar expression patterns in multiple independent microarray datasets are considered to be highly functionally correlated (Lee et al., 2004; Price and Rieffel, 2004). Therefore, we hypothesized that genes coexpressed with Hoxc8 might suggest plausible roles for Hoxc8, and that they might include cofactors or downstream target genes. In this study, 567 genes were found to be coexpressed with Hoxc8 in at least four out of eight datasets; 23 genes were determined to be coexpressed in six datasets. Among these, genes associated with ECM and cell adhesion were most significantly enriched (Fig. 2 and Supplementary Table S2), irrespective of the tissue types studied. Since the functions of known Hoxc8 target genes are biased toward cell adhesion, it is not surprising that most significantly enriched genes are associated with ECM or cell adhesion. During vertebrate development, Hox genes are expressed during gastrulation, when epithelial mesenchymal transition is initiated at the midline of the embryos, the so-called primitive streak, in order to produce a new class of cells (mesodermal cells) between the epiblasts and hypoblasts. Eventually, Hox genes regulate morphogenesis, during which ECM and cell adhesion molecules play critical roles. Consequently, it is possible that ECM and/or cell adhesion molecules could be major targets of Hox proteins. Among the ECM genes and/or cell adhesion molecules identified (Supplementary Table S2), half of them are ECM structural constituents, like collagen, fibrillin, fibulin, and laminin. Interactions between ECM and cell adhesion molecules regulate cell migration and differentiation in cooperation with growth factors and hormones (Daley and Yamada, 2013). Enrichment of ECM genes has been reported as a significant feature of each dataset analyzed (LeDouxe et al., 2006; Li et al., 2007; Timmons et al., 2007). For example, Hoxc8 was overexpressed in white adipocytes compared to brown adipocytes (GDS2743) and in intestinal mesenchymal cells compared to epithelial cells (GDS2699). The transcriptomes of white adipocytes and mesenchymal cells were specifically enriched with ECM genes (Li et al., 2007; Timmons et al., 2007). Therefore, direct or indirect regulation of ECM genes might be a good signature of Hoxc8 activation.

In addition to ECM genes, significant enrichment of transcriptional regulators is also noteworthy. These genes include known Hox cofactors, transcription factors, and chromatin remodeling factors (Supplementary Table S2). The known Hox cofactors, such as the Pbx (Pbx1 and Pbx3) and Meis (Meis1 and Meis2) classes of TALE (Three Amino acid Loop Extension) homeodomain proteins (Moens and Selleri, 2006) suggest that Hoxc8 transactivates downstream target genes by interacting with these cofactors, as other Hox proteins do to overcome poor sequence specificity. Previously, our group identified Pcna, which harbors both Pbx1 and Hoxc8 binding sites in the promoter region, as a Hoxc8 target gene (Min et al., 2010). Other interesting genes analyzed here are those with homeodomains, such as Zhx1 and Prrx1, which function as transcriptional repressors. Zhx1 was found to interact with DNA methyltransferase (DNMT) 3B (Kim et al., 2007). Interestingly, Prrx1^?/? mice showed limb bud and skeletal losses (Mann, 1995), and adipogenesis was inhibited (Du et al., 2013). Given that mutations in cofactors result in similar loss-of-function phenotypes of Hox (Moens and Selleri, 2006), Prrx1 might interact with Hoxc8 as a cofactor to repress genes during morphogenesis and adipo-genesis. We also identified many other transcription factors that play important roles in embryonic development as Hoxc8 companion genes. This finding was not observed in a previous in vitro cell line experiment (Lei et al., 2006). However, about 50% of known Hox target genes have been reported to be transcription factors in Drosophila (Hueber and Lohmann, 2008). Thus, it is possible that Hoxc8 transactivates or represses other transcription factors, which in turn amplify signaling and construct Hoxc8-mediated regulatory networks.

To validate the genes analyzed in silico, five genes (Tgfbi, Ptpn13, Prkd1, Adam19, and Aldh1a3) among the 567 were selected at random, and their concordant expression patterns were analyzed with Hoxc8. We confirmed that their expression levels were positively correlated with that of Hoxc8 in vivo in embryonic tissues and in vitro in cell lines. As an in vitro model, we used a TGF-β-induced system in MC3T3-E1 cells and overexpression/ knockdown in different cell lines. Therein, we showed that four of the selected genes (Tgfbi, Ptpn13, Prkd1, Adam19), except Aldh1a3, are coexpressed with Hoxc8 and positively regulated by Hoxc8. This allowed us to make sure that Aldh1a3 is, as we expected (coexpressed with Hoxc8 in 3 out of 8 datasets), less relevant to Hoxc8 network, whereas the other four genes, which showed coexpression in at least four of the eight datasets, sufficiently met the requirements to support our in silico methods. Evolutionally conserved Hoxc8 binding sites were found to be present in the promoter regions of the four genes (Supplementary Table S5). Practically, enrichment of Hoxc8 on Hoxc8 binding sites in Adam19 and Ptpn13 was confirmed by ChIP assays (Supplementary Fig. S2), suggesting that the list of genes reported by us include direct target genes. Based on the results, we assume that these genes can be directly or indirectly regulated by Hoxc8 and functionally linked. Especially, Adam19 and Tgfbi, whose functions are associated with ECM and cell adhesion, are likely to be Hoxc8 target genes. Adam19 is a membrane-anchored glycoprotein that plays key roles in embryo implantation, neurogenesis, cardiovascular morphogenesis, and the release of proteins, such as epidermal growth factor receptor ligands and osteoprotegerin ligand, a protein important in osteoclast differentiation and mammary gland development (Fata et al., 2000; Kong et al.,1999; Qi et al., 2009). Tgfbi is a secreted protein that is induced by treatment with transforming growth factor-β and inhibits cell attachment (Skonier et al., 1994). It is strongly expressed in the mesenchyme of numerous collagen-rich tissues throughout all stages of murine development (Ferguson et al., 2003; Schorderet et al., 2000), implicating its role in tissue morphogenesis. Together, our collective results suggest that the listed Hoxc8 companion genes are a worthy resource for exploration of Hoxc8 target genes and regulatory networks.