Comparative analysis of genes downstream of the Hoxd cluster in developing digits and external genitalia

In the course of tetrapod development, both the limbs and the external genitalia first appear as small buds of mesenchymal cells. The mouse forelimb bud begins to emerge from the lateral plate mesoderm during the ninth day of embryonic development (E9) and the genital tubercle (genital bud; genital eminence) becomes visible at E11.0. The subsequent elongation of both structures is dependent on the epithelial-mesenchymal interactions that establish a signaling center at the distal end of the growing bud: the apical ectodermal ridge (AER) in limbs and the distal urethral epithelium (DUE) in the presumptive genitals. During outgrowth, the proliferating mesenchmymal cells of the limb bud are patterned and differentiate into mature limb structures in a proximal to distal sequence (reviewed by Tickle, 2003). Although the patterning of the genitalia has not been studied as extensively, the time-dependent development of proximal and distal segments and the presence of a bony structure (os penis/baculum in males or os clitoris in females) in many mammalian genital organs suggests a common ontogenetic strategy. Therefore,despite the rather different fate of these two structures, their basic underlying morphogenetic processes are similar.

Over the past 15 years, gene expression analyses and functional approaches have largely substantiated the surprising similarity between the distal part of developing limbs (the digits) and the genitalia. Indeed, the majority of genes known to be of interest for developmental growth and patterning in presumptive digits are also transcribed in the emerging genital bud. The initial outgrowth of the limb bud is dependent on the production of retinoic acid by the enzyme ALDH1A2 (RALDH2)(Niederreither et al., 2002),whose expression also suggests a role in the formation of the genital bud. Fibroblast growth factor (FGF) signaling is required for both AER and DUE function (Haraguchi et al.,2000; Sun et al.,2002). In the limb, expression of the sonic hedgehog gene(Shh) from a posterior patch of mesodermal cells called the zone of polarizing activity (ZPA) is required for both anteroposterior and proximodistal patterning (Riddle et al.,1993). In developing genitals, Shh expression in the urethral epithelium, rather than the mesoderm, is required for outgrowth of the bud (Haraguchi et al.,2001; Perriton et al.,2002). The significance of these homologous expression patterns is confirmed by gene inactivation experiments, which often cause phenotypic effects in both structures (e.g. Dolle et al., 1993).

The patterning of the autopod domain of the limb, which corresponds to the hands and feet, is of particular interest in this context. For example, the expression patterns of the genes coding for the transcription factors of the HOX family are comparable in both the genitalia and in distal limbs(Dolle et al., 1991). Specifically, the Hoxa and Hoxd genes related to the Drosophila Abd-Bgene, located at the 5′ extremity of their respective clusters, are necessary for patterning both the digits and genitalia. The most 5′genes from each complex, Hoxa13 and Hoxd13, are the ones that are most strongly expressed at the distal ends of limbs and in the genital bud. The combined inactivation of Hoxa13 and Hoxd13in the same animal leads to the agenesis of both external genitalia and distal limbs (Fromental-Ramain et al.,1996; Kondo et al.,1997). Analyses of compound Hoxa13/Hoxd13mutants revealed that effects of these null mutations are dose-dependent in genitals, in the digestive and urogenital tracts(Warot et al., 1997), and in limbs (Zakany et al., 1997). Likewise, a variety of human congenital syndromes involving mutations in either HOXD13 or HOXA13 are characterized by morphological defects in both digits and genitalia(Goodman, 2002).

These observations raise the question of how the same set of control genes,active in similar developmental processes, ultimately generates morphologically and functionally different structures. One possibility is that Hox transcription factors regulate the expression of different target genes depending on the ontogenetic context. Interestingly, there is little information about the identity of HOX-dependent target genes during limb patterning, and their relationships with either the FGF, WNT, SHH, BMP or retinoid signaling pathways are unclear. Recent studies have begun to address these questions. In one approach, candidate genes were chosen because of their known role in the patterning pathways outlined above. In this way Knosp et al.(Knosp et al., 2004)identified Bmp2 and Bmp7 as being downstream of Hoxa13 in limb development, and a related study showed that Bmp7 and Fgf8 expression is dependent on Hoxa13 in genitals as well (Morgan et al.,2003).

A more comprehensive method is to use microarray technology because in this case no prior knowledge of candidate genes is required. One approach is to express Hox genes in cell lines where they are normally not expressed, and then compare expression profiles with and without Hox expression. In this way, Itga8 was established as a target of Hoxa11 regulation in the developing kidney (Valerius et al.,2002). A perhaps more physiologically relevant method is to compare expression in tissues from wild-type and Hox-mutant mice. With such a technique, Hedlund et al. (Hedlund et al.,2004) identified candidate target genes by comparing embryonic spinal cord tissue from wild-type and Hoxd10-null mice. In a related approach, Pruett et al. (Pruett et al.,2004) used transgenic mice overexpressing the Hoxc13 gene to identify keratin genes of the Krtap16 family as likely direct targets of HOXC13. The use of such approaches is nevertheless complicated by the severity of some Hox mutations and the functional redundancy of these genes. If the mutation is severe enough to radically transform tissues,expression profiling would undoubtedly identify many candidate genes as,effectively, different structures would be being compared. However, the relevance of these candidates would be unclear, because in this case the genes immediately downstream would be difficult to differentiate from the many genes far downstream that are necessary to form the structure. Therefore, we chose not to use a mutant with a drastic phenotype such as the Hoxa13/Hoxd13 double mutant described above. Using this or comparable strains would preclude a comparative approach as external genitalia and distal limbs do not develop in these mutants.

Instead, we analyzed a mutant with an intermediate phenotype in which the genitalia and distal limbs were partially, but not completely, lost. We used a mouse strain deleted for the entire Hoxd cluster(HoxD^Del1-13). When Hoxd genes are deleted, the digits and external genitalia are reduced in size, but the basic pattern of both structures, as well as the relative amount of cell types, apparently remains intact (Zakany and Duboule,1996; Zakany et al.,1997). Because of the common expression patterns of Hoxd genes in genitals and digits, and the related phenotypic effects of their mutations, we hypothesized that genes downstream of Hoxd could be similarly regulated in both structures. Therefore, we simultaneously identified candidate genes in distal forelimbs and genitals by analyzing global gene expression in both tissues, in both the presence and absence of Hoxd genes. The advantage of this comparative strategy is that in addition to identifying specific candidate genes, we can also begin to characterize the general nature of gene expression programs downstream of the Hoxd genes in two distinct structures, and therefore address the effects on gene expression of an evolutionary-conserved Hox expression pattern.

Candidate genes were identified by microarray analysis and subsequently validated by complementary approaches: real-time RT-PCR and in situ hybridization. Only a minority of candidate genes identified in either tissue appeared to be Hoxd dependent in just one of the two structures. This suggests that expression of the same Hox genes in two different developmental contexts modulates at least a subset of overlapping downstream genes. Furthermore many of the identified genes were not previously known to have any role in limb or genital development, and most are not part of any of the established limb or genital patterning pathways.

The HoxD^Del1-13 (previously called Del9) and HoxD^Del1-10 alleles have been previously described(Spitz et al., 2001; Zakany et al., 2004). Hoxa13-mutant mice were provided by P. Chambon and have been previously described (Fromental-Ramain et al., 1996). HoxD^Del1-13 heterozygous mice from a mixed C57BL/6-129/Sv background were crossed to give homozygous mutant and wild-type embryos. Noon on the day of plugging was assumed to be E0.5. For microarray analysis, E12.5 embryos were dissected by cutting the autopod of the forelimbs with fine forceps and by removing the genital bud with a tungsten needle. The precise point of cutting for each tissue is indicated by arrows in Fig. 1A. Samples were immediately placed in RNAlater reagent (QIAGEN) for storage during genotyping.

The RNA samples used for microarray analysis (herein referred to as primary samples) were extracted from pools of limb or genital tissue from five to 10 individual embryos of the same genotype, three replicate pools for each tissue and genotype. RNA was isolated using the RNeasy micro- or mini-kit (QIAGEN),following homogenization and disruption with a POLYTRON device (Kinematica)using the QIAGEN RLT solution. The yield of total RNA was 5-7 μg per pool of genital tissue and 20-45 μg per pool of limb tissue. The quality of all RNA pools was confirmed by analysis on a 2100 Bioanalyzer (Agilent). Secondary samples for quantitative real-time RT-PCR analysis were collected in triplicate from E11.5-E14.5 embryos and extracted as above, except that in most cases a sample from only one individual embryo was used per replicate.

Double-stranded cDNA was synthesized from 2-5 μg total RNA from each pool, according to the GeneChip Expression Analysis Technical Manual(Affymetrix), using the SuperScript Choice system (Invitrogen). cDNA was used to synthesize biotin-labeled cRNA using the BioArray HighYield RNA Transcript Labeling Kit (Enzo). After purification using a QIAGEN RNeasy column, 25 μg of cRNA was fragmented. Each fragmented cRNA (15 μg) was then hybridized to an Affymetrix U74Av2 GeneChip microarray. Hybridization, washing and scanning were performed according to the Affymetrix manual.

Data from the scanned chips were analyzed using Affymetrix MAS 5.0 software(Hubbell et al., 2002; Liu et al., 2002). The significance of differential expression was determined by the number of pairwise comparisons (out of nine total) having a significant change in the same direction. The cutoff for significance was a P-value of 0.0025 by Wilcoxon's Signed Rank test for each pairwise comparison. The fold-change values reported are the means±s.d. of the nine comparisons. The Hoxd genes were not included in the reported results, but detection of their differential expression served as controls for both the tissue samples and detection by the microarrays. Hoxd genes had the greatest fold changes measured, as much as 190-fold greater signal strength for Hoxd13 when comparing wild-type and HoxD^Del1-13 forelimb samples. Candidate genes were prioritized into categories by the number of significant change comparisons, and were categorized secondarily by the fold change. Gene lists were further filtered by eliminating those candidates that were called absent in both genotypes by the MAS 5.0 analysis. In some cases, probe sets were eliminated from the analysis when extensive BLAST searching indicated that they were not specific for an individual gene. The raw data were also analyzed with other publicly available software, including dChip(Li and Hung Wong, 2001) and GCRMA (Wu et al., 2003). The results from these complementary analyses largely confirmed the lists of candidate genes obtained by MAS 5.0. According to MIAME guidelines (Brazma et al.,2003), the complete microarray dataset was deposited in the public data repository of the European Bioinformatics Institute (ArrayExpress) with accession number E-MEXP-257.

Single-stranded cDNA templates for real-time RT-PCR analysis were synthesized from the same RNA pools used for the microarray analysis, and from independently derived and extracted secondary samples as described above. Sixty-five to 100 bp amplicons and Taqman probes were designed using Primer Express 2.0 software (Applied Biosystems). Primer pairs were tested and efficiencies were measured using standard curves from serial dilutions of cDNA. Primer pairs having greater than 88% efficiency were judged to be acceptable for subsequent measurements. Results and amplification efficiencies were comparable using either Taqman or Sybr Green. Specificity of Sybr Green reactions was determined by examination of product melting curves as described(Ririe et al., 1997). cDNA was PCR amplified in a 7900HT SDS System (Applied Biosystems) and raw threshold-cycle (Ct) values were obtained from SDS 2.0 software (Applied Biosystems). Relative quantities (RQ) were calculated with the formula RQ=E^-Ct using efficiencies (E) calculated for each run with the DART-PCR algorithm, as described (Peirson et al., 2003). A mean quantity was calculated from triplicate PCR reactions for each sample, and this quantity was normalized to two (for primary samples) or three (for secondary samples) similarly measured quantities of normalization genes as described(Vandesompele et al., 2002). Normalized quantities were averaged for three replicates for each data point and represented as the mean±s.d. The highest normalized relative quantity was arbitrarily designated as a value of 1.0. Fold changes were calculated from the quotient of means of these normalized quantities and reported as ±s.d. The statistical significance of fold-changes was determined by a paired Student's t-test. Taqman PCR was used to quantify the expression of Hoxa11, Papss2, Aldh1a2, Stra6, Msx2, Prrx1,Gfra2 and Epha3. Sybr Green PCR was used to quantify the expression of Sgk, Gdf10, Odz4, Nr2f1, Pcdha, Foxp1, Shox2 and Lisch7. Control genes for normalization, Rps9, Tbp and Tubb4, were amplified either by Sybr Green or Taqman PCR, as appropriate for each run. Primers and probe sequences can be retrieved at: http://www.unige.ch/sciences/biologie/biani/duboule/index_st.htm.

All cDNA fragments used as whole-mount situ hybridization (WISH) probes were prepared by RT-PCR using E12.5 genital or limb bud RNA, with the exception of the Hoxa11 probe (gift of Catherine Fromental-Ramain). Four hundred to 700 bp amplicons were designed with MacVector software(Accelrys) using publicly available GenBank cDNA sequences. Specificity of the primers and chosen cDNA fragment was confirmed by BLAST analysis. Specific primer sequences are listed on the website indicated above. The various cDNAs were cloned into the pGEM-T vector (Promega) and the correct identities of the cloned fragments were confirmed by restriction enzyme analysis. DIG-labeled antisense probes were prepared by in vitro transcription with SP6 or T7 polymerase (Promega). WISH using embryos fixed in 4% paraformaldehyde was performed according to standard procedures. Because sex-specific differences in external genitalia only begin to appear at E16.5(Suzuki et al., 2002), the sex of embryos was not routinely determined. Nevertheless, in some cases the gender of embryos was determined by PCR as described(Lavrovsky et al., 1998) to confirm that the expression patterns were not sex-dependent. All images represent one representative staining of at least two replicates for each condition.

We sought to identify genes regulated by HOXD proteins, while determining whether the comparable expression patterns of Hoxd genes in developing digits and genitalia regulated similar or different downstream developmental pathways. We analyzed developmental gene expression by microarray analysis,quantitative real-time RT-PCR and WISH analysis. Our primary genetic tool was the HoxD^Del1-13 mutant stock, in which all Hoxd genes are deleted (Spitz et al.,2001).

We chose to use samples from E12.5 embryos as a starting point for our analysis, as it is the first day in which the genital bud exists as a structure that can be reliably dissected, and because Hoxd gene expression is reaching high levels in both structures at this time(Fig. 1A). Also E12.5 embryos,unlike younger specimens, have presumptive autopods clearly distinguishable from more proximal domains.

Triplicate wild-type and HoxD^Del1-13 homozygote pools of E12.5 RNA were analyzed using Affymetrix microarrays that contain probe sets representing approximately 12,000 cDNAs. The resulting data sets allowed a total of nine comparisons between pools to be made for each type of tissue(Fig. 1B). Using the criteria described in Table 1, the best candidates were chosen by the number of comparisons in which differential expression was statistically significant, and by the magnitude of the fold change. Surprisingly, of the strongest candidates, in which at least seven out of nine comparisons showed significant differential expression, only four genes showed a greater than twofold change in expression(Table 1; first row). When the criteria were relaxed to include genes with a fold-change cutoff of 1.4-fold,only 13 more genes were identified with seven significant comparisons(Table 1; second row). A third,even less stringent category, yielded 28 more candidates in which the threshold for fold-change remained at 1.4-fold, but the number of significant differential comparisons was reduced to 6. Further relaxation of the criteria(Table 1; fourth row) expanded the list to hundreds of genes with less than six statistically significant comparisons. Strikingly, none of the candidate genes identified in the two tissues were found to be differentially expressed in both distal forelimbs and genitals.

The small fold changes observed made it crucial to further verify the potential differentially expressed genes by RT-PCR. We chose to validate a subset of 25 genes, including all four from the most stringent category, most from the second category (10 out of 13), and a sample of 11 genes from the two least restrictive categories. We first quantified expression of these 25 genes in our primary E12.5 samples that were used for the microarray analysis. Of these 25 genes, 16 showed a differential expression by real-time RT-PCR that was similar to that observed by microarray analysis(Fig. 1C). Our threshold for defining a gene as `validated' was a change of at least 1.4-fold by RT-PCR, in the same direction as was detected by the microarray analysis. The precise fold changes measured for validated genes are indicated in Tables 2 and 3. The percentage of candidates confirmed from the two most stringent categories in Table 1 was 78% (11 out of 14). By contrast, only 5 out of 11 (45%) genes were confirmed from the two least stringent categories, thus supporting the use of these categories to identify the most likely candidates. Nonetheless, apparent false-positives from the microarray analysis were found in each category, emphasizing the need to verify any specific gene before being assured of its differential expression.

Results comparing samples from E12.5 indicated that there were no common Hoxd-target candidate genes in distal limbs and genitalia. However, because the genital eminence starts to develop 1 to 2 days later than the forelimb, we analyzed and compared samples obtained from different stages of development,reasoning that shared molecular mechanisms could be implemented at different times in both structures. For this reason, and to further confirm the differential expression seen in our primary samples, we used real-time RT-PCR to analyze secondary samples from embryos at days E11.5 to E14.5 to generate a developmental time course of expression in genitals and distal forelimbs for each of the 16 candidate genes and two control genes. If a candidate gene were involved in digit or genital bud formation, we expected that its differential expression would be detected during this time period because digit morphogenesis is essentially complete by E14.5 and the genital bud is a well-developed structure at this time.

Because the same secondary RNA samples were to be used to characterize all candidate genes, we first assayed for any systematic differences amongst samples. Accordingly, we chose two genes, Prrx1 and Msx2,which showed no significant difference in expression by microarray analysis. Furthermore, these particular genes served as robust controls for either random or systematic effects, as both are expressed in the same mesenchymal tissues as the Hoxd genes and function in limb patterning(Leussink et al., 1995; Satokata et al., 2000). The resulting real-time RT-PCR expression profiles for these control genes indicated no differential expression in either forelimb or genital buds (see Fig. S1 in the supplementary material), supporting the use of these samples in identifying truly differentially expressed genes.

We confirmed differential expression at E12.5 for 14 out of the 16 genes in the secondary samples, giving a final list of 14 validated candidate genes(Tables 2 and 3). Therefore, despite the relatively small fold changes detected, multiple rounds of verification of these genes rigorously confirmed their differential expression. Furthermore,among the 11 candidates identified in distal limbs, seven also had highly significant (P<0.01) differential expression in genital buds on at least one day of development (Table 3; last column). With the exception of Gdf10, the most significant differential expression in genitals was found in embryos older than 12.5 days (the differential expression of Gdf10 in E12.5 genitals was apparently not detected by the microarrays because of its weak signal). Similarly, of the three candidates initially identified in genitals,two also had highly significant differential expression in limbs(Table 2; last column) at stages other than E12.5. Ultimately, nine out of the 14 candidate genes were found to be differentially expressed in both presumptive digits and genitals.

After quantifying variations in the RNA levels of these candidate target genes over developmental time, we investigated their spatial expression patterns. We designed WISH probes for each of the 14 candidate genes and stained wild-type and mutant embryos at the appropriate stages. The three genes initially identified as differentially expressed in genitals(Table 2) were among the most clearly confirmed. Firstly, although Hoxa11 was barely detectable in wild-type E12.5 genital buds by WISH, it was consistently upregulated in both the genital bud and trunk of the HoxD^Del1-13 mutant specimens (Fig. 2A). This observation confirmed the stable 2- to 3-fold elevation in Hoxa11expression measured in HoxD^Del1-13 genital buds by real time RT-PCR (Fig. 2B). By contrast, only a slight 1.3- to 1.4-fold increase in Hoxa11expression was measured in HoxD^Del1-13 limbs.

Contrastingly, the Sgk (serum glucocorticoid-regulated kinase) and Gfra2 (glial cell line-derived neurotrophic factor family receptor alpha 2) genes clearly showed significant downregulation in both HoxD^Del1-13 limbs and genitals. The downregulation in the expression of Sgk was initially detected by microarray in E12.5 genital buds. Upon WISH staining, differential expression was also clearly apparent in forelimb buds (Fig. 3A). Sgk expression first appears in a distal domain of the forelimb bud at E11. This domain was not observed in three out of three stained HoxD^Del1-13 E11-E11.5 embryos. The differential expression at E11.5 was quantified as being 2.3-fold by real-time RT-PCR(Fig. 3C). At E12.5, some expression did appear in HoxD^Del1-13 forelimbs(Fig. 3A; second panel), but only proximally, whereas a more distal interdigital domain developed in wild-type forelimbs. The proximal domain in the center of the HoxD^Del1-13 E12.5 autopod may explain why differential expression was not detected in E12.5 limbs by microarray analysis, or by real-time RT-PCR (Fig. 3C). At E13.5, the differential expression was again significant, as the strong interdigital domain seen in wild-type limbs was never scored in HoxD^Del1-13 limbs. Differential expression of Sgkwas no longer detected at E14.5, but by this time the interdigital domain of expression in wild-type forelimbs had disappeared (data not shown). A similar distal domain of Sgk expression was seen in wild-type genital buds from E11.5 (data not shown) through to E13.5, which was not observed in HoxD^Del1-13 genitals(Fig. 3B).

Similar to Sgk, Gfra2 expression was visibly reduced in both HoxD^Del1-13 mutant forelimbs and genital buds(Fig. 4). In contrast to Sgk, however, differential expression of Gfra2 was apparent first in genitalia and later in limbs, as seen both by WISH and real-time RT-PCR analysis. Robust expression of Gfra2 never developed in HoxD^Del1-13 forelimbs, at least not through E14.5, whereas a clear signal was observed in wild-type limbs at E13.5, and then only in the ventral mesenchyme of the autopod (Fig. 4A). The differential expression was first detected in genitals at E11.5 (not shown) and continued through E14.5(Fig. 4B). Although Gfra2 expression did not disappear in HoxD^Del1-13genitalia, it was clearly reduced and restricted to a more distal domain.

The expression patterns of all of the candidate genes in Table 3 were determined(Fig. 5, see also Figs S2-S4 in the supplementary material) for the developmental stages in which they showed significant differences by real-time RT-PCR. Odz4 and Epha3RNA levels were increased in both HoxD^Del1-13 limbs and genitals (Fig. 5A,B),suggesting that they could be targets for Hoxd repression (as are Shox2 and Pcdha, see Fig. S4C,D). The increase in Odz4 and Epha3 expression in mutant forelimbs was clearly detected by WISH. At E12.5, the domain of Odz4 expression increased both in its distal extent and intensity, and this increase was still visible in E13.5 and E14.5 embryos (not shown). Similarly, the domain of Epha3 expression extended slightly more distally, and was reinforced towards the posterior and proximal parts of mutant digits at E12.5(Fig. 5B). An increase of Odz4 and Epha3 in genitals was difficult to visualize by WISH, but was highly significant (P<0.01) by real-time RT-PCR in both cases (Fig. 5A,B).

In contrast to Odz4 and Epha3, Gdf10 is a candidate for Hoxd upregulation (as are the six other genes initially identified in limbs listed in Table 3). Because Gdf10 is expressed broadly in both forelimbs and genital buds from both genotypes, its loss of expression in the HoxD^Del1-13mice was difficult to document by WISH(Fig. 5C), although the difference was highly significant by real-time RT-PCR in both tissues. For other genes, such as Papss2 (see Fig. S2A in the supplementary material) and Aldh1a2 (see Fig. S3A in the supplementary material),the in situ staining confirmed the differential expression as quantified by microarray and real-time PCR, but their expression patterns were so dynamic that an interpretation of their role downstream of Hoxd was difficult. The patterns of the remaining six genes confirmed to varying degrees the differential expression (see Figs S3, S4 in the supplementary material).

In order to further validate Epha3 and Sgk as candidate target genes, we made use of additional mouse strains with specific Hox gain or loss of functions. In the case of Epha3, its expression in the most anterior part of the developing hand-plate was virtually exclusive from Hoxd gene expression in digits, suggesting a repression of the former by the latter products. We tested this possibility by using a mouse strain showing a gain of function of Hoxd13 and Hoxd12 in distal limbs(HoxD^Del1-10) (Zakany et al., 2004). In HoxD^Del1-10 mice, the expression of Hoxd13 in presumptive digits is extended anteriorly,thus overlapping with the Epha3 domain(Fig. 6A). In the HoxD^Del1-10 mutant limbs, expression of Epha3 was no longer detected, concomitantly with a strong gain of Hoxd13expression in the same domain (Fig. 6B). This result confirmed the repressive effect of Hoxd gene products over Epha3 transcription, in agreement with the observed gain of expression of Epha3 in limbs lacking all Hoxd function(Fig. 5B). In this latter case, Epha3 expression did not extend throughout developing digits, which normally express posterior Hoxd genes (Fig. 5B), suggesting that Hoxd gene products cannot be the only proteins preventing Epha3 expression in developing distal limbs.

Functional redundancy between Hoxa and Hoxd genes in limb and genital development is well established (e.g. Kondo et al., 1997; Fromental-Ramain et al., 1996; Davis et al., 1995). Hoxa13, in particular, is strongly expressed during digit development, in a domain largely overlapping with that of posterior Hoxd genes. In addition, the most distal digit domain where Epha3expression was not gained in the absence of Hoxd gene functions precisely coincides with the Hoxa13 transcript domain(Fig. 5B). Therefore, we tested whether the Epha3 gain of expression would be enlarged upon inactivation of the Hoxa13 gene in combination with the absence of Hoxd genes. Epha3 expression was gained in mutant forelimbs in a Hox-dose dependent fashion (Fig. 7). When one copy of Hoxa13 was inactivated in combination with a Hoxd-homozygous deletion, a much larger distal and posterior gain in Epha3 expression was scored than in the Hoxd mutant alone (Fig. 7C, compare with Fig. 5B). When the remaining copy of Hoxa13 was also inactivated so that the only Hox products present in the autopod were from the Hoxa11 gene, the expression domain of Epha3 was further gained to the posterior margin of the forelimb bud (Fig. 7D). Therefore, Epha3 appears to be downstream of both the Hoxa and Hoxd genes.

The gain of expression of posterior Hoxd genes in HoxD^Del1-10 mutant mice was also used to monitor expression of the Sgk gene, a candidate for transcriptional upregulation by Hoxd gene products. Unlike Epha3, Sgk expression in wild-type developing digits perfectly overlaps with that of Hoxd13(Fig. 6; compare A with C). In particular, Sgk, like Hoxd13, was not detected in the most anterior aspect of the developing distal limb(Fig. 6). In mutant HoxD^Del1-10 limbs, however, Sgk expression was detected concomitantly with the gain of expression of Hoxd13 in this anterior domain (Fig. 6A). This confirms the results obtained with both microarray and RT-PCR analyses, which placed Sgk downstream of HOXD function.

Although HOX proteins were amongst the first transcription factors whose genes were isolated in higher eukaryotes, very few Hox target genes have been described so far. This apparent paradox finds its roots in the functional history of this gene family. In the course of evolution, different Hox genes were co-opted to achieve a variety of functions in different organs or structures. Therefore the functional redundancy varies quite dramatically in different structures and reflects the phylogeny of this gene family(Duboule and Wilkins, 1998). For instance, a late co-option of a few Hox genes in a defined organ will result in a low level of redundancy, with only some genes being functional. By contrast, the ancestral functional domain of this gene family, i.e. the developing trunk, will display a maximum of redundancy, making functional approaches difficult.

An intermediate situation is observed in both limbs and genitalia, where mostly two clusters (Hoxa and Hoxd) are involved in patterning. Consequently,the level of functional redundancy observed in these structures is probably not as important as that in the developing trunk, and removing the functions of two orthologous genes leads to drastic alterations. Therefore, the search for target genes should, in principle, be easier in developing limbs and genitalia than in the developing spinal cord or sclerotomes. This is, however,not necessarily the case, as removing Hox gene functions progressively deletes the corresponding limb and genital structures, rather than inducing distinct morphological alterations.

We tried to overcome this problem by using an intermediate genetic condition in which the limbs and genitalia, while being affected in their overall size and patterning, still display fairly good morphologies due to the function of the remaining Hoxa genes. Because limbs and genitals have comparable early developmental phases, we wondered whether their different morphological fates could be accounted for by distinct transcriptional outputs following the activation of the same Hox genes, i.e. whether the same HOX proteins would trigger the activation of different genetic programs.

We identified and validated 14 genes as potential HOXD targets. Six of these (Hoxa11, Sgk, Gfra2, Epha3, Odz4 and Gdf10) are especially strong candidates. Interestingly, with the exception of Hoxa11, these genes are similarly regulated in limbs and genitals,which further indicates that these two structures display closely related developmental strategies. Initially, these genes were not found to be regulated simultaneously in both developing buds, until the analyses were extended to various time-points. This revealed that the same variations in gene regulation were observed in both structures, but usually slightly later in genitalia than in limb buds (Gfra2 is an exception to this trend,see below). This observation nicely fits the developmental delay that exists between these buds, as the genital eminence emerges with a one to two day delay with respect to forelimb budding.

Strikingly, except for Hoxa11, the function of the five other potential target genes has not yet been fully explored in limb or genital development. In fact, none of the 14 candidate genes are members of the classical FGF, BMP, WNT or SHH signaling pathways. Only the retinoid pathway(Aldh1a2, Stra6) is represented among the candidates. This observation supports the conclusion that, at the developmental stages we analyzed, the Hoxd genes act downstream or independently of most of the previously described limb and genital patterning pathways.

Post et al. (Post et al., 1999) previously reported that the expression domain of Hoxa11 extends more distally into the autopod domain when Hoxa13 is mutated. We have found a similar upregulation of Hoxa11 in genitals devoid of Hoxd gene function. Therefore Hoxa11 is a common downstream target of Hox proteins. In our study, Hoxa11 did not increase significantly in forelimbs, probably because of the presence of the Hoxa13 product. As Hoxa13 is also expressed in developing genitals, it is unclear why we observe a clear increase in Hoxa11 transcripts in the absence of Hoxd genes. Nonetheless, the combined results from our study and the findings of Post et al. (Post et al., 1999) identify Hoxa11 as similarly regulated by Hox proteins in limbs and genitals. In both cases, one may wonder whether such an upregulation of Hoxa11 could both weaken the phenotype and,accordingly, reduce the changes seen in the expression of candidate target genes through functional compensation.

Despite this clear effect on Hoxa11 regulation and a few other reported cases of auto-regulation (Popperl et al., 1995; Popperl and Featherstone, 1992), cross-regulatory and autoregulatory interactions amongst Hox genes and their products does not appear to be the rule, particularly for those posterior Hoxd genes involved in limb and genitalia development. This is demonstrated well by the inability of Hoxd transgenes to be faithfully expressed in distal limbs and genitalia whenever integrated outside the Hoxd cluster, even in the presence of the full Hox gene complement (see van der Hoeven et al.,1996).

The Sgk gene is certainly amongst the more unlikely candidates identified in this work. No developmental role for this kinase has been reported, even though the Sgk orthologous gene in C. eleganshas been identified as a critical determinant of life span and stress response(Hertweck et al., 2004). A previous study reported tissue-specific expression of Sgk during development, but an analysis of limbs beyond E10.5 was not described, thus overlooking the expression phase starting at E11.5(Lee et al., 2001). The mouse genome encodes two other Sgk isoforms (Sgk2 and Sgk3), which is likely to account for the virtual absence of a phenotype in mice null for the gene (Wulff et al., 2002).

Among protein kinases, the SGK family is most closely related to the AKT protein kinases. A developmental role for AKT kinases was not fully evident until null mutations for Akt1 and Akt2 were combined(Peng et al., 2003). Likewise,the developmental role of the Sgk genes will probably require detailed combined analysis of all three forms. Because both AKT and SGK kinases act downstream of PDK1 (3-phosphoinositide-dependent protein kinase 1), SGK could participate in the recently reported involvement of AKT in the regulation of apoptosis in the limb(Kawakami et al., 2003). Biochemical studies have shown that SGK has an anti-apoptotic function(reviewed by Lang and Cohen,2001), at least partially through its phosphorylation and inactivation of the pro-apoptotic transcription factor FKHRL1(Brunet et al., 2001; Mikosz et al., 2001). Intriguingly, the domain of Sgk expression we report, especially at E13.5 (Fig. 3), is adjacent to the domain of interdigital cell death occurring more distally(Chen and Zhao, 1998). Likewise, a domain of apoptosis has been reported in the distal genital bud(Haraguchi et al., 2001; Suzuki et al., 2003),immediately adjacent to the domain of Sgk expression. SGK could be involved in regulating the domains of apoptosis or proliferation in both structures and its misregulation could contribute to the smaller size of appendages seen in the HoxD^Del1-13 mutant.

Gfra2 codes for a receptor for the neutrophic factor neurturin,which signals through the RET receptor tyrosine kinase(Buj-Bello et al., 1997). A null mutation in this gene caused defective parasympathetic innervation of the gut and penis (Laurikainen et al.,2000; Rossi et al.,1999) but, as for Sgk, the existence of closely related family members makes full assessment of the developmental role difficult. We report here that the Gfra2 gene is highly expressed from the onset of genital budding, an expression that is markedly diminished in the absence of Hoxd genes. Although the expression pattern correlates with the role in the innervation of the genitalia, the breadth of Gfra2 expression throughout the genital bud suggests an additional developmental role as well. The expression outside of purely neuronal tissue has been previously noted,but as yet no function has been assigned(Klein et al., 1997).

The Gfra2 expression pattern we report in limbs is significantly more restricted than that in genitals, but is equally Hoxd dependent. In contrast to most of the other candidates, differential expression of Gfra2 appears later in limbs than in genitals. The only function for Gfra2 that has been reported in limbs is a postnatal requirement for the innervation of sweat glands on the ventral surface of the paws(Hiltunen and Airaksinen,2004). Although this function is likely to be required too late in development to be assigned to the expression we see at E13.5-14.5, it suggests a possible role in the innervation of this domain of the limb.

Ephrins and their receptors have a well-established role in neuronal pathfinding, cell migration and cell adhesion (reviewed by Poliakov et al., 2004). Their role in limb development has only recently begun to be explored, but the redundancy and overlapping expression domains of the ephrin and Eph genes(eight and 13 family members, respectively) makes this task particularly challenging. Compagni et al. showed that females heterozygous for a null mutation of ephrin B1 have digit duplications and bifurcations(Compagni et al., 2003). Stadler et al. initially established a link between Hox gene expression and ephrin signaling by showing that Epha7 expression is markedly lower in Hoxa13-null forelimbs (Stadler et al., 2001). Another study showed that overexpression of ephrin A2, one of the ligands of EPHA3, caused digit bifurcations and fusions in chick limbs (Wada et al.,2003).

In addition, previous studies have reported that ephrin receptor (Eph)genes are the direct targets of Hox proteins. HOXA9 has been shown to directly bind to and regulate the expression of the Ephb4 gene in cultured endothelial cells (Bruhl et al.,2004). Similarly, when co-expressed in vitro along with PBX1,HOXA1 and HOXB1 bind to and activate transcription from an enhancer sequence that is known to direct rhombomere-specific expression of Epha2(Chen and Ruley, 1998). Taken together these studies clearly establish a link between Hox genes and ephrin signaling.

Our data show that Hoxd and Hoxa genes act in combination to downregulate Epha3 expression in developing digits. However, the biological relevance of this repression is unclear in part because Epha3-null mice have no abnormal phenotype (Vaidya et al., 2003). Our microarray data indicate that at least 12 of the ephrin and Eph genes are expressed in developing limbs and genitals. Therefore, defining a precise role for Eph and ephrin genes in limb development will require comprehensive studies of all of these genes. Regardless of its function, the Epha3 gene can serve as a model for Hox regulation, as demonstrated by its complementary expression response in three kinds of Hox mutant stocks.

Like the Hox genes themselves, all of the best candidate genes identified have multiple paralogs in the mouse genome, suggesting that functional redundancy may prevent the rapid elucidation of the biological role of the downstream genes. This complex situation is not unexpected considering the global function of Hox genes during vertebrate development, which is to modulate the fate of a given morphological module, rather than triggering the activation of novel genetic pathways. For example, the difference between a cervical and a lumbar vertebra is most likely to be due to subtle modulations of the same genetic determinants, rather than to the function of distinct pathways. In the case of both the limbs and genitalia, Hox genes participate in both the elaboration and the specification of the structures, and therefore affect the entire process rather than some specific parts of it.

In this view, the difficulty to assign clear target genes to the Hox proteins, at least during trunk, limb and genital development, is not surprising, but may be a precise indication of how the system works. This`horizontal' regulatory strategy, whereby subtle variations in the amounts of related proteins impact upon the balance between a large set of products is in marked contrast with the situation found in arthropods, where Hox genes seem to be part of more `vertical' regulatory processes. In the former case,various thresholds of target production, or combinations thereof, may induce a structure to produce a given morphology rather than another, related one. In the latter case, the presence or absence of a Hox product may trigger a chain of events leading to the choice of a given genetic pathway.

The vertebrate case will be difficult to solve with our current analytical tools, given the difficulty of identifying global correlations out of multiple parameters, rather than punctual downstream effectors. The analysis of HOX protein-binding sites may help to some extent, but in vitro studies have shown that Hox proteins by themselves bind DNA rather nonspecifically, using core-binding sequences of only four nucleotides(Phelan and Featherstone,1997). The binding specificity apparently comes from Hox proteins forming complexes with other proteins, including those of the Pbx and Meis gene families, each with their own binding specificities. Indeed recent work in Drosophila has shown that a complex of at least five protein subunits represses the well-characterized Hox target gene Distalless(Gebelein et al., 2004). While some of the binding sites for higher order complexes on mammalian DNA have been determined (reviewed by Mann and Affolter, 1998), it is as yet difficult to reconcile these in vitro studies with the in vivo evidence. In particular, we do not know which complexes function in limbs and genitals, especially since MEIS1, which has been reported to be the strongest in vitro binding partner for AbdB-related HOX proteins (Shen et al., 1997), is not expressed in distal limbs(Mercader et al., 1999). Therefore, although some of the genes identified here could be direct Hox targets, demonstrating this conclusively will require future studies of the Hox transcription complexes in limbs and genitals. When appropriate antibodies become available, chromatin-immunoprecipitation experiments will be invaluable for conclusively demonstrating Hoxd interactions with specific regulatory regions. Our study provides a list of genes that can serve as substrates for these future studies.

We thank J. Zákány for mutant mice and consultations. The assistance of P. Descombes, M. Docquier, O. Schaad and D. Chollet (NCCR Genomics Platform) was essential for our microarray and real-time PCR analysis. In addition we thank P. Descombes and M. McDonald for critical reading of the manuscript, sharing reagents and help with data analysis; M. Kmita for discussions and generous sharing of mutant embryos; T. Montavon for sharing reagents; and P. Chambon for providing Hoxa13-mutant mice. This work was supported by funds from the Canton de Genève, the Claraz and Louis Jeantet foundations, the Swiss National Research Fund, the National Pole of Research (NCCR) `Frontiers in Genetics' and the European community(EUMORPHIA). J.C. was supported by an NRSA from the National Institutes of Health (USA).