Digital gene expression tag profiling of bat digits provides robust candidates contributing to wing formation

doi:10.1186/1471-2164-11-619

BMC Genomics. 2010; 11: 619.

Published online 2010 November 6. doi: 10.1186/1471-2164-11-619

PMCID: PMC3017863

Digital gene expression tag profiling of bat digits provides robust candidates contributing to wing formation

Zhe Wang,^1,² Dong Dong,^1,³ Binghua Ru,¹ Rebecca L Young,² Naijian Han,⁴ Tingting Guo,¹ and Shuyi Zhang¹

¹Institute of Molecular Ecology and Evolution, iAIR, East China Normal University, Shanghai 200062, PR China

²Yale Systems Biology Institute, Department of Ecology and Evolutionary Biology, Yale University, West Haven, CT 06516, USA

³Donnelly Centre for Cellular and Biomolecular Research, University of Toronto, 160 College Street, Toronto, ON, M5S 3E1, Canada

⁴Institute of Zoology, Chinese Academy of Sciences, Beijing 100101, PR China

Corresponding author.

Zhe Wang: zhe.email/at/gmail.com; Dong Dong: ddseesea/at/gmail.com; Binghua Ru: waoe0.1/at/163.com; Rebecca L Young: rebecca.youngbrim/at/yale.edu; Naijian Han: hannaijian/at/163.com; Tingting Guo: hugt06171003/at/163.com; Shuyi Zhang: syzhang/at/bio.ecnu.edu.cn

Received February 5, 2010; Accepted November 6, 2010.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

Abstract

Background

As the only truly flying mammals, bats use their unique wing - consisting of four elongated digits (digits II-V) connected by membranes - to power their flight. In addition to the elongated digits II-V, the forelimb contains one shorter digit (digit I) that is morphologically similar to the hindlimb digits. Here, we capitalized on the morphological variation among the bat forelimb digits to investigate the molecular mechanisms underlying digit elongation and wing formation. Using next generation sequencing technology, we performed digital gene expression tag profiling (DGE-tag profiling) of developing digits in a pooled sample of two Myotis ricketti and validated our sequencing results using real-time quantitative PCR (RT-qPCR) of gene expression in the developing digits of two Hipposideros armiger.

Results

Among hundreds of genes exhibiting significant differences in expression between the short and long digits, we highlight 14 genes most related to digit elongation. These genes include two Tbx genes (Tbx3 and Tbx15), five BMP pathway genes (Bmp3, RGMB, Smad1, Smad4 and Nog), four Homeobox genes (Hoxd8, Hoxd9, Hoxa1 and Satb1), and three other genes (Twist1, Tmeff2 and Enpp2) related to digit malformations or cell proliferation. In addition, our results suggest that Tbx4 and Pitx2 contribute to the morphological similarity and five genes (Acta1, Tnnc2, Atp2a1, Hrc and Myoz1) contribute to the functional similarity between the thumb and hindlimb digits.

Conclusions

Results of this study not only implicate many developmental genes as robust candidates underlying digit elongation and wing formation in bats, but also provide a better understanding of the genes involved in autopodial development in general.

Background

Bats are a speciose group which constitute approximately 20% of extant mammals, and are the only mammalian group that have powered flight [1]. Unlike the bird wing, made up of an elongated radius and ulna and modified wrist and hand bones that support the flight feathers, bat wings are primarily shaped by elongated digits (digits II-V) that support the wing membranes or dactylopatagia (Figure (Figure1A)1A) [2]. This digit elongation is extreme. For their body length, bats have the longest relative digit lengths with the longest digit (digit III) exceeding the length of the body. However, not all of the forelimb digits are elongated, digit I - the thumb - maintains a short morphology similar to the hindlimb digits and is used to grasp surface or move on the ground (Figure (Figure1).1). This unique morphology of the bat forelimb digits has sparked interest and motivated investigation by many developmental and evolutionary biologists. Despite these efforts, the mechanisms underlying development of these two distinct digit morphologies (i.e., short and long) of the bat hand remains unresolved. In this study, we capitalize on variation in digit development morphology within the forelimb and between the fore- and hindlimb in bats to identify candidate genes associated with development of distinct digit morphologies. Using DGE-tag profiling (described below) of the thumb, the elongated forelimb digits and the hindlimb digits we identified 14 candidate genes likely associated with digit elongation and seven candidate genes likely associated with the morphological and functional similarity between the thumb and hindlimb digits.

Figure 1

Limbs of the adult and fetal bats (modified from our previous study [18]). (A) Left limbs of adult Myotis ricketti. DI, DII, DIII, DIV and DV represent digits I-V of the forelimb; (B, C) Left limbs of Miniopterus schreibersii fuliginosus in the Fetal (more ...)

Elongation of the posterior digits (digits II-V) in the bat forelimb is initiated at an early stage when digit anlagen initially appear, embryonic Stage 15 [3]. By embryonic Stage 17, the elongation becomes obvious with lengthening of forelimb digits II-V greatly exceeding that of the thumb and hindlimb digits [3,4]. This rapid elongation of the posterior forelimb digits is thought to results from modification in expression of key developmental genes [5,6]. In fact, several experimental and empirical studies have found differences in expression of limb and digit associated genes between mouse and bat autopodia. First, a mouse transgenic experiment that replaced an enhancer of the mouse Prx1 - a transcription factor that regulates long bone growth in mouse - with the orthologous putative enhancer sequence from a bat resulted first in an increase in Prx1 expression and second in a small, but significant elongation of the mouse forelimb digits [7,8]. Second, Shh, a key factor in limb development, exhibits a unique second phase of signaling late in development of the bat wing that is not present in the development mouse autopodia [9]. Next, other developmental genes including Fgf8 and Bmp2/4 are up regulated in bat forelimb digits compared to bat hindlimb or mouse digits [10,11]. Finally, the distribution Hoxd13 expression is more posteriorly restricted in the bat forelimb relative to mouse fore- and hindlimbs and bat hindlimbs during late embryonic stages [12,13]. While these results indicate that differences in expression of Hoxd13, Fgf8, Bmp2/4 and Prx1 during forelimb digit development may contribute to bat digit elongation, the complex developmental changes required to generate the bat wing likely require numerous molecular changes [6]. Thus, further examinations of gene expression during development of the bat digits are critical to identify the genes associated with digit elongation in the posterior forelimb digits and those responsible for the morphological and functional similarities between the short thumb and hindlimb digits.

Here, using a DGE-tag profiling approach, we compared patterns of gene expression among three cDNA libraries formed from distinct regions of the fetal autopodia (the thumb, the four posterior forelimb digits and the all five hindlimb digits) in Rickett's big-footed bat, Myotis ricketti. The method of DGE-tag profiling combines the SAGE principle with the massively parallel sequencing technology to study genome-wide gene expression profiles [14]. DGE-tag profiling has three main advantages compared to the conventional microarray technology. First, it is high-throughput sequencing. As a result, data from a single lane of a flow cell (consisting of eight lanes allowing for eight samples) is enough to analyze of genome wide gene expression including weakly expressed genes which cannot be assessed using microarray [15-17]. Next, this type of sequencing data is highly replicable and accurate such that technical replicates have little variation [16]. Finally, sequencing can be performed without a priori sequence information, which is required in microarrays [17]. These features of DGE-tag profiling aided our identification of genes associated with digit morphology in bats. Overall, our examination of differentially expressed genes revealed candidate genes associated with both digit elongation as well as morphological and functional similarities between the thumb and hindlimb digits. Finally, our results revealed genes associated with digit position in both limbs.

Results

Specimen features

All fetuses obtained were in the "Fetal Stage" which is the last prenatal stage according to bat embryonic staging systems [3]. At this stage, the overall appearance of the embryo is similar to the former stages, but the body size as well as the digit length increases rapidly [3,18]. The genital tubercle had become a vagina or a penis allowing determination of specimen sex. All the specimens used in this study had relatively weak pigmentation indicating an early Fetal Stage as overall body pigmentation increases as fetal development progresses [3,18]. Claws were sharp and keratinized in the thumb and hindlimb digits. Crown rump length (CRL) and forearm length (FA) for both species were recorded in Table Table11.

Table 1

Specimen id and measurements

Statistics of tag sequencing

In total, 4377261, 4650734 and 7514943 reads (not including the adapter tags) were sequenced from Library Hand DI, Library Hand DII-V and Library Foot, respectively. Library Hand DI contains 375701 distinct tags, Library Hand DII-V contains 483043, and Library Foot contains 532850. These distinct tags and their genomic frequency as well as the raw data were deposited in NCBI Gene Expression Omnibus (GEO) database with the accession number [GEO: GSE20038]. Most of the distinct tags (over 62%) appeared only once (frequency = 1); however, a small number of distinct tags (less than 1.7%) with frequency higher than 100 made up over 63% of all the tags in all three libraries (Additional file 1, Table S1). When sequencing depths reach 1,000,000 total tags, the number of distinct tags discovered (especially those with frequency >1) drops dramatically in all the three libraries (Figure (Figure2).2). From that point, increasing sequence depth results in a slow and stable accumulation of new distinct tags indicating that the sequencing has reached saturation (Figure (Figure22).

Figure 2

Saturation measured as the number of new distinct tags added with increasing sequencing depth in each library Hand DI, Hand DII-V and Foot. FRQ indicates the frequency of distinct tags in the library.

Differences in expression of genes among the three libraries are indicated in Figure Figure3.3. In total, 7080, 8205 and 8967 human genes and 4742, 5666 and 6379 mouse genes were mapped using the tags from the three libraries Hand DI, Hand DII-V and Foot. In total, 38.5% human genes (10490 genes) and 24.6% mouse genes (7843 genes) were mapped using the tags from the three libraries (Additional file 2, Table S2 and additional file 3, Table S3). Some tags mapped to human genes but not mouse genes or vice versa.

Figure 3

Pairwise comparison of distinct tag expression between all three libraries (Hand DI, Hand DII-V and Foot) using DEGseq software. Genes exhibiting significant differences in expression are provided in red (p < 0.001), genes with similar expression (more ...)

DGE-tag profiling of marker gene expression

The expression patterns of five marker genes (Tbx5, Pitx1, Meis2, Hoxd10 and Krt17) in the three DGE-tag profiling libraries are consistent with previously reported patterns of expression at embryonic earlier stages (reviewed in the Discussion; Table Table2).2). We found that expression of Tbx5 was significantly higher in forelimb digits than the hindlimb digits (p < 0.001) and Pitx1 expression was significantly lower in forelimb than in hindlimb digits (p < 0.001). Similarly, as expected, we found that expression of Meis2 and Hoxd10 were significantly higher in Library Hand DII-V relative to the other two libraries (p < 0.001) and Krt17 was significantly lower in Library Hand DII-V relative to the other two libraries (p < 0.001).

Table 2

DEGseq analysis of marker genes as well as those highlighted in this study related to digit development

Differences in gene expression across the libraries

Additional files 2 and 3 list the results of DEGseq analysis for all the homologs of human (Additional file 2, Table S2) and mouse genes (Additional file 3, Table S3). For the comparison between elongated digits and the short digits, we found that 286 homologs of human genes and 89 homologs of mouse genes were up regulated in Library Hand DII-V compared with other two libraries (p < 0.001, Figure Figure4)4) and 97 homologs of human genes and 75 homologs of mouse genes were down regulated in Library Hand DII-V compared with other two libraries (p < 0.001, Figure Figure4).4). For the comparison between forelimb digits and the hindlimb digits, we found that 187 homologs of human genes and 163 homologs of mouse genes were down regulated in Library Foot compared with other two libraries (p < 0.001, Figure Figure4)4) and 198 homologs of human genes and 84 homologs of mouse genes were up regulated in Library Foot compared with other two libraries (p < 0.001, Figure Figure44).

Figure 4

The number of genes differentially expressed between distinct digit morphologies and digit limb association (p < 0.001). (A) Comparison between the elongated wing digits (Library Hand DII-V) and the short digits (Libraries Hand DI and Foot). (B) (more ...)

Among the differentially expressed genes, we found 14 genes likely related to digit elongation based on results of previous studies (see Discussion). The expression of Tbx3, Tbx15, Bmp3, Rgmb, Smad1, Smad4, Hoxd8, Hoxd9, Satb1, Hoxa1, Twist1, Tmeff2 and Enpp2 were significantly higher in Library Hand DII-V compared with other two libraries, and the expression of Nog was significantly lower in Library Hand DII-V (Table (Table2).2). Moreover, our results hightlighted 7 genes were likely related to the morphological and functional similarity between the thumb and the hindlimb digits. The levels of Tbx4, Pitx2, Acta1, Tnnc2, Atp2a1, Hrc and Myoz1 were significantly higher in Libraries Hand DI and Foot than in Library Hand DII-V (Table (Table22).

Real-time quantitative PCR (RT-qPCR)

To validate the results of DGE-tag profiling, we performed RT-qPCR of 14 genes which were differentially expressed between libraries (Table (Table2).2). All 14 genes exhibited the same overall expression pattern using the distinct methods (Figure (Figure5).5). Meis2, Hoxd10, Tbx3, Tbx15, Bmp3, Rgmb, Hoxa1, Twist1, Tmeff2 and Enpp2 showed significantly higher expression in the posterior forelimb digits (HDII-V) compared with the thumb (HDI) and all five hindlimb digits (Foot). Krt17, Tbx4, Acta1 and Tnnc2 showed significantly higher expression in both the thumb (HDI) and all five hindlimb digits (Foot) than in the posterior forelimb digits (HDII-V). The fold change in expression among the samples was often lower in RT-qPCR as compared to DGE-tag profiling; however, this is not surprising given the differences between the two methods in sensitivity of gene detection as well as the fact that the methods were performed in different species.

Figure 5

Patterns of gene expression in hand digit I (HDI), hand digits II-V (HDII-V), and all five foot digits (Foot). To validate DGE-tag profiling, expression of 14 genes, exhibiting differences in expression among the M. ricketti libraries, were examined in (more ...)

Discussion

Massively parallel sequencing is a powerful technology enabling extensive investigation of gene expression profiles. In this study, we obtained over 16.5 million sequences of expressed tags from three bat autopodial libraries, including the thumb (Hand DI), the posterior forelimb digits (Hand DII-V), and all five hindlimb digits (Foot). Using these data, we compared gene expression profiles among the libraries to identify candidate genes associated with elongation of the posterior forelimb digits (Hand DII-V). In addition, we examined gene expression patterns associated with the short digit morphology of the thumb and hindlimb digits as well as gene expression patterns associated with the fore- versus the hindlimb digits. The results of DGE-tag profiling were validated by examining gene expression in a subset of genes using RT-qPCR.

To assess the reliability of our gene expression profiles, we first examined expression of several genes with previously established expression patterns. We compared expression of Tbx5 and Pitx1, fore- and hindlimb marker genes, in the forelimb and hindlimb digits (Libraries Hand DI and Hand DII-V and Library Foot). Consistent with previous studies, we found that expression of Tbx5 was significantly higher in the forelimb digits and Pitx1 expression was significantly higher in the hindlimb digits although the embryonic stages used here are much later than those used in previous studies [19,20]. In addition, we compared expression of Meis2 and Hoxd10 in the fore- and hindlimb digits. Previously published microarray data revealed high levels of expression of these two genes in bat autopodia at Stages 16 and 17 [21]. We found that expression of both genes, Meis2 and Hoxd10, remained significantly higher in Library Hand DII-V in the Fetal Stage (Table (Table2).2). Finally, Krt17, a gene expressed in nails and important for the development of nails, exhibited significantly higher expression in the thumb and hindlimb digits [22]. Because only the thumb (Hand DI) and the hindlimb digits have claws, we would expect higher expression of Krt17 in the thumb and hindlimb digits relative to the posterior forelimb digits (Hand DII-V) that make up the wing.

In addition to comparing our gene expression profiles to genes with known expression patterns, we further validated our results by performing RT-qPCR using 14 genes which were identified as being differentially expressed between elongated digits and short digits in DGE-tag profiling. With the goal of identifying genes related to wing formation across bat taxa, we collected samples from a different species which belongs to the other suborder of Chiroptera (H. armiger). We selected 14 genes likely important for wing formation or for the similarity between the thumb and foot digits (see below). All the 14 genes exhibited the same overall expression pattern in the RT-qPCR as in the DGE-tag profiling (Figure (Figure5).5). This finding not only provides strong support for the DGE-tag profiling method used in this study, but also maintenance of gene expression patterns across these two phylogenetically distant species provides strong support for the role of these genes in digit elongation in the posterior forelimb digits or morphological and functional similarities between the thumb and hindlimb digits.

Overall, comparisons of gene expression profiles between digit morphologies and limbs identified hundreds of differentially expressed genes. Several interesting patterns have emerged from this data. Specifically, we highlight 21 genes likely related to wing formation or to morphological and functional similarities between thumb and hindlimb digits (Table (Table2).2). First, we found 14 genes that are likely associated with digit elongation in bats - two Tbx genes (Tbx3 and Tbx15), five genes from the BMP pathway (Bmp3, Rgmb, Smad1, Smad4 and Nog), four Homeobox genes (Hoxd8, Hoxd9, Satb1 and Hoxa1) and three other gene (Twist1, Tmeff2 and Enpp2) related to either digit malformation or cell proliferation. Next, we identified seven genes (Tbx4, Pitx2, Acta1, Tnnc2, Atp2a1, Hrc and Myoz1) that are likely associated with the morphological and functional similarities between the thumb and hindlimb digits. Expression patterns of these genes in the three distinct digit libraries as well as known function of these genes are described in detail below.

Genes associated with digit elongation

In general, Tbx3 and Tbx15 belong to T-box gene family and are involved in skeletogenesis [23,24]. In digit formation specifically, Tbx3 is thought to modulate posterior digit identity through Shh and BMP signaling, and Tbx15 insufficiency causes Cousin Syndrome which includes congenital dwarfism, moderate brachydactyly (or shortening of the digits) and leg shortening [25-27]. Mutants of Tbx15 have reduced bone size and changes in bone shape in the forelimb skeleton indicating an important role in the skeletal development of forelimbs [24].

Bmp3, Rgmb, Smad1, Smad4 and Nog are all components of the BMP pathway, a pathway critical for bone formation [28]. Bmp3 belongs to the transforming growth factor-beta (TGFβ) superfamily and stimulates cell proliferation and differentiation in a concentration dependent fashion [29]. Overexpression of Bmp3 in the chick wing bud has be shown to result in an increase in cell proliferation resulting in lengthening of skeletal elements [30]. Rgmb, also named as Dragon, enhances BMP signaling by directly binding to BMP proteins including Bmp2 [31]. Overexpression of Bmp2 and Bmp4 in the developing chick limb leads to a dramatic increase in the volume of cartilage elements [32]. In addition, Bmp2/4 has been shown to be up regulated in the bat forelimb embryonic digits relative to both mouse forelimb digits and bat hindlimb digits [11]. This up-regulation is correlated with increased rate of cartilage proliferation [11]. Beyond up-regulation of BMPs themselves, Smad1 mediates BMP signaling. Specifically, in its phosphorylated state Smad1 forms a complex with Smad4 and then accumulates in the nucleus. This SMAD complex regulates transcription of several target genes [28]. Embryonic manipulation of Smad1 has confirmed that its expression at the phalanx-forming region influences digit identity in the chick limb [33]. In addition, cartilage-specific knockout of Smad1 and Smad4 reduce chondrocyte proliferation and increase apoptosis resulting in defects of limb elements in the mouse [34,35]. While all the components of the BMP pathway mentioned above stimulate cartilage formation, Nog is an antagonist of BMPs and as such inhibits cartilage formation [36]. Mutations of Nog disturb the balance of BMP signaling leading to brachydactyly [37]. Thus, our results - exhibiting up regulation of Bmp3, Rgmb, Smad1, Smad4 and down regulation of Nog in the elongated bat forelimb digits - suggest a role for BMP signaling in the lengthening of the posterior bat forelimb digits.

Hoxd8, Hoxd9, Hoxa1 and Satb1 are homeobox genes which specify the anterior-posterior axis and have been implicated as digit identity determining genes [38]. Although Hoxd8 and Hoxd9 mainly contribute to proximal limb development and are only mildly expressed in the developing digits, we found that Hoxd8 and Hoxd9 were highly expressed in bat wing digits II-V relative to the thumb and hindlimb digits [39]. Our results also indicate that the homeobox gene Satb1 and Hoxa1 were highly expressed in bat wing digits II-V relative to the thumb and hindlimb digits. These two genes are highly expressed in cancer cells promoting tumor growth [40,41]. To our knowledge, there are no studies examining expression of Satb1 and Hoxa1 in limb development making them interesting genes for future investigation.

Similarly to Satb1 and Hoxa1, Tmeff2 and Enpp2 are up regulated in carcinomas [42,43]. Previous study showed Tmeff2 contributes to cell proliferation. Enpp2 which is expressed in precartilaginous condensations can be strongly activated by Hoxa13 and Hoxd13 which are crucial factors for digit identity and development [44]. Finally, knowledge of the role of Twist1 and the digit development comes from pathological studies [45]. Mutations Twist1 have been found in patients with Saethre-Chotzen syndrome which is associated with digit malformation [45]. Overexpression of Twist1 inhibits osteoblast differentiation and prevents premature fusion of skeletons [46].

Genes associated with similarities between the thumb and hindlimb digits

In addition to genes related to bat digit elongation, we found several genes that may contribute to the morphological and functional similarities between the bat thumb and hindlimb digits. First, Tbx4 and Pitx2 are primarily associated with hindlimb development. Tbx4 has also been shown to regulate bone size. Specifically, in the developing mouse forelimb, overexpression of these genes results in shortening of the forelimb skeletons [47,48]. Next, we found several genes that are likely associated with the distinct musculature and function of the short digits. Because the clawed thumb of bats is used to cling to the cave while roosting and the clawed hindlimb digits are used for hanging upside down these digits should be richer in skeletal muscles than the elongated digits. In fact, expression of several skeletal muscle associated genes - Acta1, Tnnc2, Atp2a1, Hrc and Myoz1 - was significantly higher in the thumb and hindlimb digits then in the elongated posterior forelimb digits [49].

Conclusions

Using genome-wide sequencing method we have identified many genes differentially expressed in the developing digits of the bat autopodia. Several of these genes are highly related to digit formation and thus likely play important roles in the development of bat wing. Our results provide robust candidate genes for future functional tests and molecular evolution comparison with other mammals to fully understand the mechanisms of wing formation and the evolution of bat flight.

Methods

Animal collection and breeding

All procedures were carried out in accordance with the Policy on the Care and Use of Animals, approved by the Ethical Committee, East China Normal University. Myotis ricketti, Rickett's big-footed bat, belongs to suborder Yangochiroptera and is widely distributed in China [50]. Their feet - which are used to forage on fishes - are larger than those of other bat species (Figure (Figure1A)1A) [51]. As in many species of bats, copulation in Myotis ricketti occurs in autumn and spermatozoa are stored in the female's uterus and oviducts through the winter hibernation. Ovulation and subsequent fertilization begin in the spring with parturition occurring in summer [52]. Because gravid M. ricketti could not be found during the summer breeding season, we capitalized on sperm storage by females and captured 10 hibernating female M. ricketti from a cave (39°42'N, 115°43'E) in Beijing on February 13, 2008. In addition to Myotis ricketti, we captured three gravid Hipposideros armiger from a cave (30°20'N, 117°50'E) in Anhui province of China on April 24, 2009 to use for quantitative real-time PCR confirmation sequencing results. H. armiger, the great leaf-nosed bat, belongs to the other suborder Yinpterochiroptera and is widely distributed in subtropical and tropical zones of Asia [53].

After capture, bats were kept in a large cage (1.2 × 1 × 1 m) covered with wire netting allowing them to hang. The cage was kept in a dark, temperature controlled room. The temperature was maintained between 18 and 24°C. Water was freely available and meal worms mixed with a powdered multivitamin and calcium tablets were provided daily.

Specimen processing

M. ricketti in the captive colony were euthanized by decapitation on May 28 and June 6, 2008 (Table (Table1).1). H. armiger were euthanized on May 20, 2009 (Table (Table1).1). Fetuses from the gravid females were put into ice cold PBS. The thumb, the remaining forelimb digits (metacarpi, phalanges and dactylopatagium) and all five hindlimb digits (metatarsi and phalages) were dissected from the fetuses (Figure (Figure1B1B and and1C)1C) and stored in liquid nitrogen until use. All other components of the fetus were fixed overnight in Bouin's fluid, washed with several changes of 70% ethyl alcohol and stored at room temperature. Crown rump length (CRL) and forearm length (FA) were measured for further estimation of the embryonic stages.

Digital gene expression tag profiling (DGE-tag profiling)

To reduce sex-specific gene expression patterns, we pooled the samples from a male and a female fetus of M. ricketti (Bat M1 and Bat M2 respectively, see Table Table1).1). Tag libraries were prepared using the DGE-Tag Profiling for NlaIII Sample Prep kit from Illumina according to the manufacturer's instructions. In brief, total RNA was extracted from forelimb digit I, forelimb digits II-V and all five hindlimb digits using TRIzol^®Reagent (Invitrogen, Cat. No. 15596026, US), treated with DNase I (Roche, Cat. No. 04716728001, Germany) and converted to double stranded cDNA using Oligo dT beads. Subsequently, the cDNA samples were digested using the restriction enzyme NlaIII, which recognizes and cuts the most 3' "CATG". Once digested into fragments, cDNA samples were ligated to Illumina specific adapter A, which contains a recognition site for enzyme MmeI. The Illumina specific adapter B was ligated after MmeI digestion (Illumina, Cat. No. FC-102-1005, US). Using this technique, three DGE libraries (Libraries Hand DI, Hand DII-V and Foot) from the three samples (Figure (Figure1B1B and and1C)1C) were constructed such that each molecule in the library is a 21 bp tag derived from a single transcript with Illumina specific adapters attached to both ends. Next, the 21 bp tags were enriched using PCR primers that anneal to the adaptors and clusters were created on a flow cell using a Solexa Cluster Station following the manufacturer's instructions (Illumina, Cat. No. FC-103-1004, US). Finally, massively parallel sequencing-by-synthesis was performed on the Illumina Genome Analyzer. Image analysis, base calling, extraction of 21 bp tags, and tag counting were performed using the Illumina pipeline. The number of times that a unique tag sequence is detected represents the quantitative expression level of the corresponding transcript in the tissue. The raw data, including tag sequences and counts, were deposited in NCBI Gene Expression Omnibus (GEO) database with the accession number [GEO: GSE20038]

In silico analysis

We analyzed saturation of each of the three libraries by examining the addition of new distinct tags associated with the increase of sequencing depth (frequency = 1, > 1 and all distinct tags). Because little gene annotation exists for the bat genome, all distinct tags were aligned against the well annotated NCBI human and mouse Refseq cDNA sequences using SOAP2 [54,55]. Taking into account the differences between species, one base pair mismatches were allowed during the mapping of tags to cDNA sequences. Tags that matched to more than one gene were not used in analyzing differential expression among the libraries. The identified cDNA sequences were mapped to NCBI Entrez GeneID. Differential gene expression analysis was performed using DEGseq software [56]. We compared libraries Hand DII-V with Hand DI and Foot to identify differentially expressed genes between elongated digits and short digits. We also compared libraries Hand DII-V and Hand DI with Foot to identify differentially expressed genes between the forelimb and hindlimb digits.

Real-time quantitative PCR (RT-qPCR)

Total RNA were extracted from forelimb digit I, forelimb digits II-V and all five hindlimb digits of two H. armiger (Bat H1 and H2, Table Table1),1), using TRIzol^®Reagent (Invitrogen, Cat. No. 15596026, US) and treated with DNase I (Roche, Cat. No. 04716728001, Germany). cDNA synthesis and qPCR was performed using the SYBR^®PrimeScript^®RT-PCR Kit (TaKaRa, Cat. No. DRR063A, Japan) on a Applied Biosystems 7300 Real-Time PCR System (ABI, US). RT-qPCR for each H. armiger (Bat H1 and H2, Table Table1)1) was performed in duplicate. β-actin was taken as a reference gene and it was amplified by a pair of previously published primers [57]. We designed 14 pairs of degenerate primers to amplify 14 target genes (Meis2, Hoxd10, Tbx3, Tbx15, Bmp3, Rgmb, Twist1, Hoxa1, Tmeff2, Enpp2, Krt17, Tbx4, Acta1 and Tnnc2). PCR products of each gene were ligated into a pMD^®19-T vector (TaKaRa, Cat. No. D102B, Japan), cloned, and sequenced using the Big Dye Terminator kit on an ABI 3730 DNA sequencer (Applied Biosystems, US). Using the obtained sequences, we designed gene specific primers for each target gene for qPCR. 2^-ΔΔ^CTmethod was used to quantify gene expression, data are expressed as mean ± SD [58]. Student's t-test was performed to examine the difference of gene expression between two different samples. Sequences of forward and reverse primers for normal PCR and RT-qPCR of each gene are shown in additional file 4, Table S4. The partial CDS of the genes were deposited at GenBank under accession number HM777024 - HM777038.

Abbreviations

CRL: Crown rump length; FA: forearm length; FRQ: frequency.

Authors' contributions

ZW participated in design of the study. ZW and DD carried out the data analysis and wrote the manuscript. ZW, BR and TG performed the experiments. RLY contributed to the final manuscript preparation. NH fed and examined the bats. SZ conceived of the study, participated in its design, and supervised the work. All authors read and approved the final manuscript.

Supplementary Material

Additional file 1

Table S1.xls. Tag frequency. The number of distinct and total tags counted at different frequencies in each library. FRQ indicates the frequency of distinct tags in the library.

Click here for file^{(16K, XLS)}

Additional file 2

Table S2.xls. Expression and analysis of human gene homologs in bat digits The expression of human gene homologs mapped by tags in each library and the results of DEGseq analysis for the study of gene differential expression [56]. Symbol indicates NCBI official gene symbol; GeneID indicates NCBI Entrez Gene ID. The raw count of each mapped gene is presented under its library name (Libraries Hand DI, Hand DII-V or Foot). All other parameters (log2(fold_change), log2(fold_change) nomalized, z-score, p-value, q-value and signature) were calculated using the DEGseq [56].

Click here for file^{(4.5M, XLS)}

Additional file 3

Table S3.xls. Expression and analysis of mouse gene homologs in bat digits The expression of mouse gene homologs mapped by tags in each library and the results of DEGseq analysis for the study of gene differential expression [56]. Symbol indicates NCBI official gene symbol; GeneID indicates NCBI Entrez Gene ID. The raw count of each mapped gene is presented under its library name (Libraries Hand DI, Hand DII-V or Foot). All other parameters (log2(fold_change), log2(fold_change) nomalized, z-score, p-value, q-value and signature) were calculated using the DEGseq [56].

Click here for file^{(3.3M, XLS)}

Additional file 4

Table S4.xls. Sequences of primers used in PCR and RT-qPCR Nucleotide sequences of the primers used in PCR and RT-qPCR of candidate genes underlying digit distinct morphologies (Figure 5). ^aindicates previous published primers [57]. The lengths of introns are inferred from the homologs of human genes.

Click here for file^{(24K, XLS)}

Acknowledgements

We thank Junpeng Zhang, Yinan Wang and Qian Yao for bat capture and care. This work was funded by a grant awarded to S. Zhang under the Key Construction Program of the National "985" Project and "211" Project.

References

Hunter P. The nature of flight - The molecules and mechanics of flight in animals. Embo Reports. 2007;8(9):811–813. doi: 10.1038/sj.embor.7401050. [PubMed] [Cross Ref]
Hedenstrom A, Johansson LC, Spedding GR. Bird or bat: comparing airframe design and flight performance. Bioinspir Biomim. 2009;4(1):015001. doi: 10.1088/1748-3182/4/1/015001. [PubMed] [Cross Ref]
Cretekos CJ, Weatherbee SD, Chen CH, Badwaik NK, Niswander L, Behringer RR, Rasweiler JJ. Embryonic staging system for the short-tailed fruit bat, Carollia perspicillata, a model organism for the mammalian order Chiroptera, based upon timed pregnancies in captive-bred animals. Dev Dyn. 2005;233(3):721–738. doi: 10.1002/dvdy.20400. [PubMed] [Cross Ref]
Hockman D, Mason MK, Jacobs DS, Illing N. The Role of Early Development in Mammalian Limb Diversification: A Descriptive Comparison of Early Limb Development Between the Natal Long-Fingered Bat (Miniopterus natalensis) and the Mouse (Mus musculus) Developmental Dynamics. 2009;238(4):965–979. doi: 10.1002/dvdy.21896. [PubMed] [Cross Ref]
Sears KE. Molecular determinants of bat wing development. Cells Tissues Organs. 2008;187(1):6–12. doi: 10.1159/000109959. [PubMed] [Cross Ref]
Cooper KL, Tabin CJ. Understanding of bat wing evolution takes flight. Gene Dev. 2008;22(2):121–124. doi: 10.1101/gad.1639108. [PubMed] [Cross Ref]
Behringer RR, Rasweiler JJt, Chen CH, Cretekos CJ. Genetic regulation of Mammalian diversity. Cold Spring Harb Symp Quant Biol. 2009;74:297–302. [PubMed]
Cretekos CJ, Wang Y, Green ED, Martin JF, Rasweiler JJ, Behringer RR, Progra NCS. Regulatory divergence modifies limb length between mammals. Genes & Development. 2008;22(2):141–151. [PubMed]
Hockman D, Cretekos CJ, Mason MK, Behringer RR, Jacobs DS, Illing N. A second wave of Sonic hedgehog expression during the development of the bat limb. Proc Natl Acad Sci USA. 2008;105(44):16982–16987. doi: 10.1073/pnas.0805308105. [PubMed] [Cross Ref]
Cretekos CJ, Deng JM, Green ED, Rasweiler JJ, Behringer RR. Isolation, genomic structure and developmental expression of Fgf8 in the short-tailed fruit bat, Carollia perspicillata. Int J Dev Biol. 2007;51(4):333–338. doi: 10.1387/ijdb.062257cc. [PubMed] [Cross Ref]
Sears KE, Behringer RR, Rasweiler JJ, Niswander LA. Development of bat flight: morphologic and molecular evolution of bat wing digits. Proc Natl Acad Sci USA. 2006;103(17):6581–6586. doi: 10.1073/pnas.0509716103. [PubMed] [Cross Ref]
Chen CH, Cretekos CJ, Rasweiler JJ, Behringer RR. Hoxd13 expression in the developing limbs of the short-tailed fruit bat, Carollia perspicillata. Evol Dev. 2005;7(2):130–141. doi: 10.1111/j.1525-142X.2005.05015.x. [PubMed] [Cross Ref]
Ray R, Capecchi M. An examination of the Chiropteran HoxD locus from an evolutionary perspective. Evol Dev. 2008;10(6):657–670. doi: 10.1111/j.1525-142X.2008.00279.x. [PubMed] [Cross Ref]
Velculescu VE, Zhang L, Vogelstein B, Kinzler KW. Serial analysis of gene expression. Science. 1995;270(5235):484–487. doi: 10.1126/science.270.5235.484. [PubMed] [Cross Ref]
't Hoen PA, Ariyurek Y, Thygesen HH, Vreugdenhil E, Vossen RH, de Menezes RX, Boer JM, van Ommen GJ, den Dunnen JT. Deep sequencing-based expression analysis shows major advances in robustness, resolution and inter-lab portability over five microarray platforms. Nucleic Acids Res. 2008;36(21):e141. doi: 10.1093/nar/gkn705. [PMC free article] [PubMed] [Cross Ref]
Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008;18(9):1509–1517. doi: 10.1101/gr.079558.108. [PubMed] [Cross Ref]
Feng L, Liu H, Liu Y, Lu Z, Guo G, Guo S, Zheng H, Gao Y, Cheng S, Wang J, Power of deep sequencing and agilent microarray for gene expression profiling study. Mol Biotechnol. pp. 101–110. [PubMed] [Cross Ref]
Wang Z, Han N, Racey PA, Ru B, He G. A comparative study of prenatal development in Miniopterus schreibersii fuliginosus, Hipposideros armiger and H. pratti. BMC Dev Biol. 2010;10(1):10. doi: 10.1186/1471-213X-10-10. [PMC free article] [PubMed] [Cross Ref]
Logan M, Tabin CJ. Role of Pitx1 upstream of Tbx4 in specification of hindlimb identity. Science. 1999;283(5408):1736–1739. doi: 10.1126/science.283.5408.1736. [PubMed] [Cross Ref]
Rodriguez-Esteban C, Tsukui T, Yonei S, Magallon J, Tamura K, Izpisua Belmonte JC. The T-box genes Tbx4 and Tbx5 regulate limb outgrowth and identity. Nature. 1999;398(6730):814–818. doi: 10.1038/19769. [PubMed] [Cross Ref]
Mason M, Hockman D, Jacobs D, Illing N. Differences in the wing and hindlimb transcriptomes of the natal long-fingered bat, Miniopterus natalensis, during embryonic development. Mech Dev. 2009;126:S17–03. doi: 10.1016/j.mod.2009.06.1015. [Cross Ref]
Mclean WHI, Rugg EL, Lunny DP, Morley SM, Lane EB, Swensson O, Doppinghepenstal PJC, Griffiths WAD, Eady RAJ, Higgins C. et al. Keratin-16 and Keratin-17 Mutations Cause Pachyonychia-Congenita. Nature Genetics. 1995;9(3):273–278. doi: 10.1038/ng0395-273. [PubMed] [Cross Ref]
King M, Arnold JS, Shanske A, Morrow BE. T-genes and limb bud development. Am J Med Genet A. 2006;140(13):1407–1413. [PubMed]
Singh MK, Petry M, Haenig B, Lescher B, Leitges M, Kispert A. The T-box transcription factor Tbx15 is required for skeletal development. Mech Dev. 2005;122(2):131–144. doi: 10.1016/j.mod.2004.10.011. [PubMed] [Cross Ref]
Nissim S, Allard P, Bandyopadhyay A, Harfe BD, Tabin CJ. Characterization of a novel ectodermal signaling center regulating Tbx2 and Shh in the vertebrate limb. Dev Biol. 2007;304(1):9–21. doi: 10.1016/j.ydbio.2006.12.010. [PMC free article] [PubMed] [Cross Ref]
Suzuki T, Takeuchi J, Koshiba-Takeuchi K, Ogura T. Tbx genes specify posterior digit identity through Shh and BMP signaling. Dev Cell. 2004;6(1):43–53. doi: 10.1016/S1534-5807(03)00401-5. [PubMed] [Cross Ref]
Lausch E, Hermanns P, Farin HF, Alanay Y, Unger S, Nikkel S, Steinwender C, Scherer G, Spranger J, Zabel B. et al. TBX15 mutations cause craniofacial dysmorphism, hypoplasia of scapula and pelvis, and short stature in Cousin syndrome. Am J Hum Genet. 2008;83(5):649–655. doi: 10.1016/j.ajhg.2008.10.011. [PubMed] [Cross Ref]
Anderson GJ, Darshan D. Small-molecule dissection of BMP signaling. Nat Chem Biol. 2008;4(1):15–16. doi: 10.1038/nchembio0108-15. [PubMed] [Cross Ref]
Carrington JL, Chen P, Yanagishita M, Reddi AH. Osteogenin (bone morphogenetic protein-3) stimulates cartilage formation by chick limb bud cells in vitro. Dev Biol. 1991;146(2):406–415. doi: 10.1016/0012-1606(91)90242-U. [PubMed] [Cross Ref]
Gamer LW, Ho V, Cox K, Rosen V. Expression and function of BMP3 during chick limb development. Dev Dyn. 2008;237(6):1691–1698. doi: 10.1002/dvdy.21561. [PMC free article] [PubMed] [Cross Ref]
Samad TA, Rebbapragada A, Bell E, Zhang Y, Sidis Y, Jeong SJ, Campagna JA, Perusini S, Fabrizio DA, Schneyer AL. et al. DRAGON, a bone morphogenetic protein co-receptor. J Biol Chem. 2005;280(14):14122–14129. doi: 10.1074/jbc.M410034200. [PubMed] [Cross Ref]
Duprez D, Bell EJ, Richardson MK, Archer CW, Wolpert L, Brickell PM, Francis-West PH. Overexpression of BMP-2 and BMP-4 alters the size and shape of developing skeletal elements in the chick limb. Mech Dev. 1996;57(2):145–157. doi: 10.1016/0925-4773(96)00540-0. [PubMed] [Cross Ref]
Suzuki T, Hasso SM, Fallon JF. Unique SMAD1/5/8 activity at the phalanx-forming region determines digit identity. Proc Natl Acad Sci USA. 2008;105(11):4185–4190. doi: 10.1073/pnas.0707899105. [PubMed] [Cross Ref]
Zhang J, Tan X, Li W, Wang Y, Wang J, Cheng X, Yang X. Smad4 is required for the normal organization of the cartilage growth plate. Dev Biol. 2005;284(2):311–322. doi: 10.1016/j.ydbio.2005.05.036. [PubMed] [Cross Ref]
Retting KN, Song B, Yoon BS, Lyons KM. BMP canonical Smad signaling through Smad1 and Smad5 is required for endochondral bone formation. Development. 2009;136(7):1093–1104. doi: 10.1242/dev.029926. [PMC free article] [PubMed] [Cross Ref]
Groppe J, Greenwald J, Wiater E, Rodriguez-Leon J, Economides AN, Kwiatkowski W, Affolter M, Vale WW, Belmonte JC, Choe S. Structural basis of BMP signalling inhibition by the cystine knot protein Noggin. Nature. 2002;420(6916):636–642. doi: 10.1038/nature01245. [PubMed] [Cross Ref]
Lehmann K, Seemann P, Silan F, Goecke TO, Irgang S, Kjaer KW, Kjaergaard S, Mahoney MJ, Morlot S, Reissner C. et al. A new subtype of brachydactyly type B caused by point mutations in the bone morphogenetic protein antagonist NOGGIN. Am J Hum Genet. 2007;81(2):388–396. doi: 10.1086/519697. [PubMed] [Cross Ref]
Montavon T, Le Garrec JF, Kerszberg M, Duboule D. Modeling Hox gene regulation in digits: reverse collinearity and the molecular origin of thumbness. Genes Dev. 2008;22(3):346–359. doi: 10.1101/gad.1631708. [PubMed] [Cross Ref]
Zakany J, Duboule D. The role of Hox genes during vertebrate limb development. Curr Opin Genet Dev. 2007;17(4):359–366. doi: 10.1016/j.gde.2007.05.011. [PubMed] [Cross Ref]
Zhang X, Emerald BS, Mukhina S, Mohankumar KM, Kraemer A, Yap AS, Gluckman PD, Lee KO, Lobie PE. HOXA1 is required for E-cadherin-dependent anchorage-independent survival of human mammary carcinoma cells. J Biol Chem. 2006;281(10):6471–6481. doi: 10.1074/jbc.M512666200. [PubMed] [Cross Ref]
Han HJ, Russo J, Kohwi Y, Kohwi-Shigematsu T. SATB1 reprogrammes gene expression to promote breast tumour growth and metastasis. Nature. 2008;452(7184):187–193. doi: 10.1038/nature06781. [PubMed] [Cross Ref]
Glynne-Jones E, Harper ME, Seery LT, James R, Anglin I, Morgan HE, Taylor KM, Gee JM, Nicholson RI. TENB2, a proteoglycan identified in prostate cancer that is associated with disease progression and androgen independence. Int J Cancer. 2001;94(2):178–184. doi: 10.1002/ijc.1450. [PubMed] [Cross Ref]
Nouh MA, Wu XX, Okazoe H, Tsunemori H, Haba R, Abou-Zeid AM, Saleem MD, Inui M, Sugimoto M, Aoki J. et al. Expression of autotaxin and acylglycerol kinase in prostate cancer: association with cancer development and progression. Cancer Sci. 2009;100(9):1631–1638. doi: 10.1111/j.1349-7006.2009.01234.x. [PubMed] [Cross Ref]
Williams TM, Williams ME, Kuick R, Misek D, McDonagh K, Hanash S, Innis JW. Candidate downstream regulated genes of HOX group 13 transcription factors with and without monomeric DNA binding capability. Dev Biol. 2005;279(2):462–480. doi: 10.1016/j.ydbio.2004.12.015. [PubMed] [Cross Ref]
Firulli BA, Krawchuk D, Centonze VE, Vargesson N, Virshup DM, Conway SJ, Cserjesi P, Laufer E, Firulli AB. Altered Twist1 and Hand2 dimerization is associated with Saethre-Chotzen syndrome and limb abnormalities. Nat Genet. 2005;37(4):373–381. doi: 10.1038/ng1525. [PMC free article] [PubMed] [Cross Ref]
Bialek P, Kern B, Yang X, Schrock M, Sosic D, Hong N, Wu H, Yu K, Ornitz DM, Olson EN. et al. A twist code determines the onset of osteoblast differentiation. Dev Cell. 2004;6(3):423–435. doi: 10.1016/S1534-5807(04)00058-9. [PubMed] [Cross Ref]
Marcil A, Dumontier E, Chamberland M, Camper SA, Drouin J. Pitx1 and Pitx2 are required for development of hindlimb buds. Development. 2003;130(1):45–55. doi: 10.1242/dev.00192. [PubMed] [Cross Ref]
Holmberg J, Ingner G, Johansson C, Leander P, Hjalt TA. PITX2 gain-of-function induced defects in mouse forelimb development. BMC Dev Biol. 2008;8:25. doi: 10.1186/1471-213X-8-25. [PMC free article] [PubMed] [Cross Ref]
Welle S, Brooks AI, Delehanty JM, Needler N, Bhatt K, Shah B, Thornton CA. Skeletal muscle gene expression profiles in 20-29 year old and 65-71 year old women. Exp Gerontol. 2004;39(3):369–377. doi: 10.1016/j.exger.2003.11.011. [PubMed] [Cross Ref]
Ma J, Dai Q, Zhang SY, Shen JX, Liang B. Distribution of Ricketti's big-footed bat (Myotis ricketti) Sichuan J Zool. 2003;22(3):155–156.
Ma J, Zhang J, Liang B, Zhang L, Zhang S, Metzner W. Dietary characteristics of Myotis ricketti in Beijing, North China. J Mammal. 2006;87:339–344. doi: 10.1644/05-MAMM-A-183R1.1. [Cross Ref]
Wang Z, Liang B, Racey PA, Wang Y, Zhang S. Sperm storage, delayed ovulation and menstruation of the female Rickett's big-footed bat (Myotis ricketti) Zool Stud. 2008;47:215–221.
Simmons NB. In: Mammal species of the world: a taxonomic and geographic reference. Wilson DE, Reeder DM, editor. Baltimore: Johns Hopkins University; 2005. Order Chiroptera; pp. 312–529.
Li RQ, Yu C, Li YR, Lam TW, Yiu SM, Kristiansen K, Wang J. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25(15):1966–1967. doi: 10.1093/bioinformatics/btp336. [PubMed] [Cross Ref]
Pruitt KD, Maglott DR. RefSeq and LocusLink: NCBI gene-centered resources. Nucleic Acids Res. 2001;29(1):137–140. doi: 10.1093/nar/29.1.137. [PMC free article] [PubMed] [Cross Ref]
Wang L, Feng Z, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–138. doi: 10.1093/bioinformatics/btp612. [PubMed] [Cross Ref]
Chen J, Sun M, Liang B, Xu A, Zhang S, Wu D. Cloning and expression of PDK4, FOXO1A and DYRK1A from the hibernating greater horseshoe bat (Rhinolophus ferrumequinum) Comp Biochem Physiol B Biochem Mol Biol. 2007;146(2):166–171. doi: 10.1016/j.cbpb.2006.10.095. [PubMed] [Cross Ref]
Schmittgen TD, Livak KJ. Analyzing real-time PCR data by the comparative C(T) method. Nat Protoc. 2008;3(6):1101–1108. doi: 10.1038/nprot.2008.73. [PubMed] [Cross Ref]

Articles from BMC Genomics are provided here courtesy of
BioMed Central