Abstract The molecular events leading to the development of the bat wing remain largely unknown, and are thought to be caused, in part, by changes in gene expression during limb development. These expression changes could be instigated by variations in gene regulatory enhancers. Here, we used a comparative genomics approach to identify regions that evolved rapidly in the bat ancestor, but are highly conserved in other vertebrates. We discovered 166 bat accelerated regions (BARs) that overlap H3K27ac and p300 ChIP-seq peaks in developing mouse limbs. Using a mouse enhancer assay, we show that five Myotis lucifugus BARs drive gene expression in the developing mouse limb, with the majority showing differential enhancer activity compared to the mouse orthologous BAR sequences. These include BAR116, which is located telomeric to the HoxD cluster and had robust forelimb expression for the M. lucifugus sequence and no activity for the mouse sequence at embryonic day 12.5. Developing limb expression analysis of Hoxd10-Hoxd13 in Miniopterus natalensis bats showed a high-forelimb weak-hindlimb expression for Hoxd10-Hoxd11, similar to the expression trend observed for M. lucifugus BAR116 in mice, suggesting that it could be involved in the regulation of the bat HoxD complex. Combined, our results highlight novel regulatory regions that could be instrumental for the morphological differences leading to the development of the bat wing.

Author Summary The limb is a classic example of vertebrate homology and is represented by a large range of morphological structures such as fins, legs and wings. The evolution of these structures could be driven by alterations in gene regulatory elements that have critical roles during development. To identify elements that may contribute to bat wing development, we characterized sequences that are conserved between vertebrates, but changed significantly in the bat lineage. We then overlapped these sequences with predicted developing limb enhancers as determined by ChIP-seq, finding 166 bat accelerated sequences (BARs). Five BARs that were tested for enhancer activity in mice all drove expression in the limb. Testing the mouse orthologous sequence showed that three had differences in their limb enhancer activity as compared to the bat sequence. Of these, BAR116 was of particular interest as it is located near the HoxD locus, an essential gene complex required for proper spatiotemporal patterning of the developing limb. The bat BAR116 sequence drove robust forelimb expression but the mouse BAR116 sequence did not show enhancer activity. These experiments correspond to analyses of HoxD gene expressions in developing bat limbs, which had strong forelimb versus weak hindlimb expression for Hoxd10-11. Combined, our studies highlight specific genomic regions that could be important in shaping the morphological differences that led to the development of the bat wing.

Citation: Booker BM, Friedrich T, Mason MK, VanderMeer JE, Zhao J, Eckalbar WL, et al. (2016) Bat Accelerated Regions Identify a Bat Forelimb Specific Enhancer in the HoxD Locus. PLoS Genet 12(3): e1005738. https://doi.org/10.1371/journal.pgen.1005738 Editor: Claude Desplan, New York University, UNITED STATES Received: July 10, 2015; Accepted: November 23, 2015; Published: March 28, 2016 Copyright: © 2016 Booker et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: ChIP-seq data is available as part of BioProject PRJNA252737 in the sequence read archive (http://www.ncbi.nlm.nih.gov/sra) as experiment ID SRX793524. Funding: This work was supported by NICHD award R01HD059862, an NSF predoctoral fellowship (to TF), and funding from Gladstone Institutes. NA is also supported in part by NIDDK award number 1R01DK090382, NINDS award number 1R01NS079231 and NCI award 1R01CA197139. NI received funding from the University of Cape Town, and the National Research Foundation (South Africa). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Vertebrate limbs show a wide range of morphological variety ranging from fins to limbs. The developing tetrapod limb is made up of three skeletal elements: the stylopod (humerus/femur), zeugopod (ulna/tibia, radius/fibula), and autopod (carpals/tarsals; metacarpals/metatarsals; phalanges) [1,2]. Autopods are highly specialized, composed of different numbers and lengths of digits, and exhibit varying degrees of interdigital soft tissue (webbing). Autopods are a hallmark of tetrapod diversity and are essential for adaptation to life on land, in the sea and in the air. Bats are an extreme example of this. To form a wing, bat forelimbs have gone through three major changes: elongation of digits II-V, retention of membranous tissue forming the inter-digital patagia (chiropatagium) and a relative reduction in the diameter of the ulna [3–5]. These morphological innovations are clearly apparent in bat fossils from 52.5 million years ago [6,7]. The genetic changes that led to the development of these specialized limb structures and mammalian flight are likely to have occurred prior to the radiation of the Chiroptera, one of the most diverse mammalian orders. Enhancers can regulate spatiotemporal gene expression during vertebrate development [8] and nucleotide changes within them can lead to phenotypic differences, such as limb malformations [9]. For example, regulatory regions in the 5’Hoxd locus have been implicated in digit specification during mammalian autopod development and loss of interactions with these regions can result in limb phenotypes, similar to Hoxd10-Hoxd13 deletions [10]. Nucleotide changes in enhancers have also been linked to morphological differences between species [11]. One such example is the Prx1 limb enhancer. The replacement of the mouse sequence of this enhancer with the homologous bat Prx1 sequence resulted in mice with longer forelimbs [12]. The recent availability of several bat genomes (Myotis lucifugus, Myotis davidii, Pteropus vampyrus, and Pteropus alecto) [13–16] now make it possible to identify specific nucleotide changes in the bat lineage, as compared to other mammals, that could have a role in the development of the unique limb morphology of the bat. Various computational approaches have been used to identify regulatory elements that could be involved in species-specific morphological changes [17–21]. These include human accelerated regions (HARs) and human accelerated conserved noncoding sequences (HACNSs), which are highly conserved sequences that have acquired a disproportionate number of nucleotide substitutions since humans diverged from our common ancestor with chimpanzees [20,22,23]. Based on epigenetic marks, Capra and colleagues predicted that at least 30% of these noncoding HARs are developmental enhancers [24]. So far, 62 out of 92 tested HARs have shown enhancer activity in mouse transgenic assays, and 7 out of 26 HARs, where the activity of the human and chimp sequences were compared, showed differential enhancer activity [25]. These include the limb enhancer sequences HAR2/HACNS1, which showed no limb specific activity for the non-human homologous sequence [22], and 2xHAR.114, which displayed restricted limb activity for the human sequence compared to the chimpanzee sequence [24]. These findings indicate that the identification of accelerated regions could serve to detect sequences that function as gene regulatory elements and could possibly give rise to characteristic phenotypes among species. Here, we set out to identify enhancers whose alteration could have contributed to bat wing development. We utilized mouse limb-specific ChIP-seq datasets and available bat genomes to identify bat accelerated regions (BARs). We identified 166 BARs with potential enhancer activity in developing limbs and functionally analyzed five of them in mouse, finding all five M. lucifugus cloned BARs (BAR2, BAR4, BAR61, BAR97, BAR116) to be functional limb enhancers. Comparison of the enhancer activity of mouse and M. lucifugus orthologous BAR sequences revealed expression differences for three of the four tested sequences (BAR4, BAR97, and BAR116), suggesting that these sequences could be accelerated in bats due to functional differences. Amongst them, M. lucifugus BAR116, which resides in a gene desert on the telomeric side of the HoxD locus, showed robust forelimb and weak hindlimb expression, a trend similar to bat Hoxd10 and Hoxd11 gene expression as we determined using whole-mount in situ hybridization on bat and mouse embryos.

Discussion Using comparative genomics, developing limb ChIP-seq datasets and mouse enhancer assays, we characterized genomic regions in bat genomes that could play a role in mediating gene expression changes underlying the unique morphological development of the bat wing. We identified 166 BARs which showed a global enrichment for Nr2c2, Zfp281, Zfp740, Zic2/3 and Egr1, and depletion for Osr1, Osr2, Tgif1 and Meis1 TFBS when comparing their mouse sequences to the inferred ancestral bat sequences (S2 Table). Analysis of five M. lucifugus BARs using a mouse transgenic assay showed all of them to have enhancer activity in the developing limb at E12.5. Examination of the mouse orthologous sequences for four of these BARs showed three to be differentially expressed compared to the M. lucifugus sequence, including BAR116 that showed strong FL versus HL expression, matching the differential expression pattern observed for Hoxd10 and Hoxd11 in developing bat limbs. The M. lucifugus BAR4 sequence showed enhancer expression in the proximal FL [68], while the mouse orthologous sequence BAR4 was negative for enhancer activity at this time point (S1 and S2 Figs). BAR4 is near Spry1, a gene involved in skeletal and muscle development in mice that when overexpressed leads to chondrodysplasia, a skeletal disorder leading to arrested development [57,69]. Interestingly, 5 of our 166 identified BARs reside in the Spry1 locus (Fig 2, S1 Table), suggesting that altered regulation of this gene could have been important for bat wing development. It is also worth noting that all five Spry1 BARs are located in the same topological associating domain [TAD; [70]](Fig 2). The presence of multiple BARs in this region could implicate a co-operative and tightly coordinated regulation of Spry1 during bat limb myogenesis and digit elongation via Fgf signaling [29,57]. The transcription factor Plag1 could be differentially binding these enhancers and affecting the regulation of Spry1 or other nearby genes. Interestingly BAR61, the well-characterized Shh limb enhancer (ZRS), did not show differential enhancer activity between M. lucifugus and mice at E12.5. Both BAR61 sequences had similar enhancer activity in the ZPA of the autopod (S1 and S2 Figs), however, both also seemed to have an expanded domain compared to the previously characterized E11.5 expression pattern [59]. These analyses do not exclude differences in expression pattern that may occur at later stages of development, or quantitative differences, which cannot be picked up through these transient transgenic assays. BAR116 is of extreme interest, due to its telomeric location relative to the HoxD locus (Fig 2). The HoxD locus is conserved across vertebrates and has a critical spatiotemporal role in skeletal development [67]. Hoxd10 and Hoxd13 in particular are known to directly interact with Shh during limb outgrowth [71,72] and mutations in these genes lead to various limb malformations including synpolydactyly, split hand and foot, and distally located skeletal elements [67,73]. Several cis regulatory elements for the HoxD cluster have been shown to drive limb development and are conserved among vertebrates [10,74]. The M. lucifugus BAR116 portrayed enhancer activity in the FL mesenchyme but had reduced activity in HL (5/5 embryos). Its orthologous mouse sequence was negative for enhancer activity, despite showing evolutionary conservation in vertebrates. It is worth noting that a 700 bp partially overlapping region to mouse BAR116, termed CNS9 found in the HoxD telomeric region (Fig 2), was previously examined for enhancer activity in a lentiviral-mediated mouse transgenesis system and was also negative [75]. M. lucifugus BAR116 showed strong enhancer activity throughout the developing FL and weak expression in the proximal portion of the HL. Using whole-mount in situ hybridizations, we analyzed the expression patterns of Hoxd10-13 at CS15-CS17 for M. natalensis and mice at matching developmental time points. Our results showed an overall similarity to previously characterized Hoxd10-13 bat and mouse limb expression patterns [15,76]. The robust FL expression and weaker HL expression of BAR116 showed a similar trend to that of Hoxd10 and Hoxd11 (Fig 6). Interestingly, Hoxd9 also showed a reduction in HL expression at CS16 [15], similar to the one we observed for Hoxd10 and Hoxd11, suggesting that BAR116 could possibly be regulating these and other HoxD genes that are located 5’ to these genes. Additionally, by using 4C combined with HoxD telomeric deletion assays in mice, Hoxd9-Hoxd11 were suggested to interact, during early phase of HoxD expression, with telomeric enhancers promoting forearm/arm development [75]. BAR116 (CNS9) lies within this telomeric region and is 850kb from the telomeric boundary of this topological domain suggesting that it could regulate HoxD genes. There are also two functional limb enhancers on both sides of it, CNS39 (60kb centromeric to CNS9; Fig 2) and CNS65 (229kb telomeric to CNS9; Fig 2), that are thought to interact with HoxD genes [75]. The differential enhancer activity we observed for M. lucifugus BAR116 compared to the negative enhancer activity of its mouse sequence and bat-mouse BAR116 composite sequence at E12.5, could imply that bats have acquired a novel enhancer function or a temporal specific chromatin conformation essential for forelimb morphology during autopod development in this locus. Our data suggests that accelerated regions could be used to identify species-specific developmental enhancers that serve as critical determinants during morphoevolution as has been demonstrated previously [20,22,24]. In bats, an examination of the evolution of echolocation identified Foxp2 as a major constituent of vocal and orofacial development pathways, finding conserved noncoding elements in close proximity of Foxp2 that changed significantly in echolocating bats when compared to non-echolocating species [77]. Similarly, our study utilizes different species to test orthologous sequences, but focuses specifically on regions that are predicted to be developmental limb enhancers through ChIP-seq. There have been previous reports that analyzed bat-specific limb enhancers [12,76]. However, our study is the first to examine, in a genome-wide manner, putative limb enhancers in bats. Interestingly, the Prx1 known bat limb enhancer whose replacement in the mouse led to longer forelimbs [12], was not identified in our analyses as a BAR element. It is worth noting that our study also had many caveats. The mouse transgenic enhancer assay that we used is not quantitative and can be inconsistent due to differences in integration sites and transgene copy number. We also could not test sequences in bat embryos and so were limited to observing expression changes only in mice. For our in situs, due the scarcity of these embryos, we were only able to use M. natalensis embryos whereas in our mouse transgenic assays we used sequences from M. lucifugus, which could also lead to differences in expression patterns. Moreover, we only examined 5 of 166 BARs, which does not represent the majority of the BARs found in our pipeline. For our global TFBS analysis, we analyzed what we determined to be the ancestral bat sequence that could lead to us missing several TFBS changes. In addition, the TFBS matrixes we used were mainly human and mouse based and could differ in bats. In M. lucifugus BAR sequences we also noticed repetitive regions that were not present in the mouse BARs, which could explain the enrichment for specific TFBSs during our motifDiverge analysis. Despite these caveats, our study shows that the use of tissue specific ChIP-seq datasets combined with sequence acceleration can be an efficient means to identify sequences that are important in determining morphological changes between species.

Materials and Methods Ethics statement Mouse work was approved by the UCSF Institutional Animal Care and Use Committee (protocol number AN100466) and was conducted in accordance with AALAC and NIH guidelines and also by the University of Cape Town Faculty of Health Sciences animal ethics committee application number FHS AEC 012/052. Ethical approval to collect bats was given by the University of Cape Town, Faculty of Science Animal Experimentation Committee (2006/V4/DJ, 2008/V16/DJ and 2012/V39/NI) with permission to sample granted by the Western Cape Nature Conservation Board (AAA004-00030-0035 (2006), AAA007-00041-0056 (2012–2014). Computational analyses To identify BARs, we employed a statistical phylogenetic test for accelerated nucleotide evolution in the common ancestor of all extant bats. This is an extension of a previously proposed likelihood ratio test for acceleration in a single species or clade [28]. This new ancestral lineage version of the likelihood ratio test is implemented in the PhyloP function (option–-branch) in the open source software package PHAST [78]. The input to PhyloP is a multiple sequence alignment for each genomic region to be tested for acceleration, plus a phylogenetic tree of the species in the alignment that is estimated from genome-wide data (in this case, four-fold degenerate sites). To apply this statistical test to bat limb development, we first identified a collection of candidate enhancers for limb development genes by intersecting evolutionarily conserved elements with enhancer-associated histone modifications and transcription factor binding events measured in the developing mouse limbs (Fig 1). Specifically, we took the union of all peaks from two previously published ChIP-seq experiments targeting H3K27ac or p300 [18,27] and an H3K27ac dataset generated for this project. Next, we generated a set of vertebrate conserved elements that were agnostic to the rate of nucleotide substitutions in bats. We started with 60-way vertebrate multiple sequence alignments with mouse as the reference species (UCSC Genome Browser, mm10 assembly). We dropped the two bat genomes (M. lucifigus and P. vampyrus) from the alignments to ensure that high rates of nucleotide differences between the bats and other vertebrates would not prevent us from identifying conservation in other species. Finally, we ran the PhastCons program with default settings [26] on the resulting genome-wide alignments. This analysis identified 4,384,943 conserved elements, many of which were less than 100 bp long and, thus, too short for statistical tests for acceleration [28]. However, we observed that many short elements frequently clustered together on the chromosome and that known functional elements (e.g., coding exons) were often tiled with multiple conserved elements separated by short gaps. Hence, we iteratively merged adjacent elements until the ratio of the distance between the elements merged over the total length of the region was less than or equal to 0.1. This merging algorithm was the result of empirical experiments aimed at producing one or a small number of merged elements per exon. We also experimented with adjusting the parameters of PhastCons to produce longer elements, but found that post-processing, by merging, recapitulated exons more effectively. Next, we intersected all merged regions greater than 100 bp with the ChIP-seq peaks and unmasked the M. lucifigus and P. vampyrus sequences from the multiple alignments. Regions with more than 50% missing sequence from either bat or more than 25% of nucleotides overlapping a coding exon were dropped to produce a collection of 20,057 candidate limb enhancers. Prior to PhyloP analysis, we integrated sequences from two additional bat genomes into the candidate enhancer alignments. We obtained assembled contigs for two bats, M. davidii and P. alecto, that were sequenced to high coverage (100x)[13]. We used the BLAST algorithm to identify alignments of the mouse sequence from each candidate enhancer to contigs from M. davidii and P. alecto [79]. The single best hit with an e-value less than or equal to 0.01 was then blasted back to the mouse genome. If this produced a reciprocal best hit (i.e., the top scoring alignment to the mouse genome overlapped the original candidate enhancer sequence), we added the M. davidii or P. alecto sequence to the 60-way multiple alignment for that candidate enhancer. This produced alignments with between two and four bats present per enhancer. The two additional bat species were added to the phylogenetic tree corresponding to the 60-way alignments (UCSC Genome Browser) and their branch lengths were adjusted using their relationship to M. lucifigus and P. vampyrus. We then restricted our analysis to regions containing at least one bat. Finally, we used PhyloP to test each candidate enhancer for accelerated nucleotide substitutions along the ancestral bat lineage. The resulting p-values were adjusted for multiple testing using a false discovery rate (FDR) controlling procedure [80,81]. We call all candidate enhancers with FDR < 5% Bat Accelerated Regions (BARs) (S1 Table). Their genomic distribution and sequence composition were analyzed using custom Python scripts. Significant associations with functions and phenotypes of nearby genes were identified using GREAT after lifting BARs over to mm9 coordinates [82]. We curated a list of limb-associated genes by exhaustively looking through the literature for evidence found in mouse or human and used resampling tests to assess associations between BARs and these genes compared to random sets of PhastCons elements. TFBS analyses To look for TFBS differences, we manually curated a list of limb-associated TFs (S2 Table). BARs were analyzed for loss and gain of binding sites for each TF using motifDiverge [35]. We first compared the ancestral bat sequence to mouse. We used prequel to computationally infer the sequence of the common ancestor of extant bats using our multiple alignments [78]. We created the corresponding aligned mouse sequence from these alignments. We then called a TFBS a hit if its FDR exceeded a threshold of 0.01. We then used motifDiverge [35] to test if the total number of TFBS in the bat ancestor was significantly different than the number of TFBS in mouse for each TF in each individual BAR. We repeated these tests collectively over all BARs. ChIP-seq ChIP was performed using the LowCell# ChIP kit (Diagenode) as previously described [83]. About 70,000 cells were pooled per IP and sonicated with a Covaris sonicator (S220 Focused-ultrasonicator, Covaris). Of the sheared chromatin, 30ul was used for each ChIP experiment with the antibody anti-acetyl histone H3 (Lys27) clone CMA309 (Milipore 05–1334). Following the manufacturer’s directions, each library was constructed using the Rubicon ThruPLEX library construction kit. Each library included 10 ul of ChIP material for a total of 14 cycles of amplification. Sequencing was carried out using an Illumina HiSeq and FASTQ files were aligned to the Mus Musculus genome (mm9) using Bowtie 0.12.8 [84]. A single base pair mismatch was permitted and reads with multiple alignments were discarded. The ChIP-seq library was sequenced to a depth of 168M total reads with 137M aligning uniquely. The input sample was sequenced to a depth of 111M reads total and 81M aligning uniquely. In each case, approximately 18% of sequences failed to align. We sorted and indexed the alignments using SAMtools 0.1.18 [85] and then converted to BED files with the bam2bed utility, a part of bedtools 2.17.0 [86]. To identify enriched H3K27ac islands in the limb samples, the peak-finding tool SICER 1.1 [87] was used. Mouse transgenic enhancer assays PCR was carried out either on M. lucifugus or M. musculus DNA using primers that were designed to amplify candidate enhancer peak sequences with additional 100–500 bp outside of predicted regions (S1 Table). The bat-mouse BAR116 composite sequence was synthesized (Biomatik) and sequence validated. PCR products and the synthetic sequence were cloned into the Hsp68-LacZ vector [63] and sequence verified. All transgenic mice were generated by Cyagen Biosciences using standard procedures [88], and harvested and stained for LacZ expression at E12.5 as previously described [89]. Pictures were obtained using an M165FC stereo microscope and a DFC500 12-megapixel camera (Leica). To be designated as an enhancer, we required consistent spatial expression patterns present in at least two embryos. Bat and mouse whole mount in situ hybridization. M. natalensis embryos were collected from the De Hoop Nature Reserve, Western Cape Province, South Africa in 2008 and 2012 and were harvested, staged and stored as described previously [64,66,90]. Whole-mount in situ hybridization was performed on wild-type Parkes wild-type strain (PKS) mouse embryos. Bat specific in situ probes were generated using primers that were designed to the M. lucifugus genomic sequence (Myoluc 2.0) and synthesized by Source Bioscience. PCR was performed using cDNA (S3 Table) and amplicons were gel extracted using the Wizard SV Gel PCR Clean-up System (Promega). Purified products were cloned using the pGEM-T Easy Vector System (Promega) in XL-Blue cells. Purified plasmids (Promega Pure Yield™ Plasmid Miniprep System) from positively selected colonies were sequenced by Source Bioscience to confirm insert identity. NCBI blast analyses confirmed sequence similarity to the expected bat transcripts (S3 Table). Mouse probes (mHoxd10, mHoxd11, mHoxd12, mHoxd13) were generated from previously published probe templates [91]. Plasmids were linearized by digestion; bHoxd10: NcoI, bHoxd11-13: SphI, mHoxd10: EcoRI, mHoxd11: SalI, mHoxd12: Bam, mHoxd13: pVUII (Thermo Fisher Scientific) and purified to form an antisense probe template. DIG-labeled probes were generated using In Vitro Transcription (IVT). SP6 polymerase (Roche) was used for all bat probes and mHoxd10 and T7 polymerase (Roche) was used for mHoxd11-13 as per the manufacturer’s instructions. Reactions were DNase treated using DNA-free (Ambion, Life Technologies) and purified using a SigmaPrep Spin Column (Sigma-Aldrich). Embryos were sectioned along the sagittal plane to conserve samples and allow direct comparisons of different probes. Whole mount in situ hybridization was performed as described by [92] including Proteinase K digestion [93] and an overnight hybridization step at 70°C followed by detection using NBT-BCIP. Accession numbers ChIP-seq data has been made publically available through NCBI (ChIP-seq BioProject ID: PRJNA252737 as experiment ID SRX793524).

Supporting Information S1 Fig. M. lucifugus BAR LacZ PCR positive embryos. Insets showing higher magnification images of all embryos that had limb LacZ staining are shown next to the whole embryo. https://doi.org/10.1371/journal.pgen.1005738.s001 (PDF) S2 Fig. Mouse BAR and bat-mouse BAR116 composite LacZ PCR positive embryos. Insets showing higher magnification images of all embryos that had limb LacZ staining are shown next to the whole embryo. https://doi.org/10.1371/journal.pgen.1005738.s002 (PDF) S1 Table. BARs identified through our computational pipeline. BAR ID, ChIP-seq datasets where they originated from, ChIP-seq peak coordinates (mm10), PhastCons coordinates (mm10), PhyloP scores, p-Value (see computational methods), False Discovery Rates (FDR), and all genes within 1Mb upstream and downstream of BAR element are shown. The 38 BARs near limb-specific genes are highlighted in red font. The five mouse BARs that were found within gene deserts (gene desert defined by distance to a TSS that is greater than 500kb in either direction) are depicted with green highlighting. The M. lucifugus and mouse BARs chosen for enhancer assays are in bold and their primer sequences used for cloning are in the subsequent worksheet. To identify the nearest gene, we used LiftOver of the M. lucifugus sequence to mm9 coordinates. https://doi.org/10.1371/journal.pgen.1005738.s003 (XLSX) S2 Table. Limb-associated transcription factor binding site analysis. The cloned sequences for M. lucifugus and mouse BARs were analysed for TFBS enrichment and depletion of known limb-associated TFs using motifDiverge [35,95]. We defined a TFBS if it exceeded an FDR of 0.01. TFBSs were tested for significant enrichment or depletion (Adjusted p-value < 0.05) in each cloned M. lucifugus BAR compared to its corresponding cloned mouse BAR. We also inferred the sequence of the common ancestor of the four extant bats for each of the 166 BARs using the prequel program in the PHAST package and compared these to their corresponding mouse sequences. We asked if TFBSs within these sequences were significantly enriched or depleted. All motifDiverge tests were corrected for multiple testing (to control the false discovery rate with the Benjamini-Hochberg method). Lists of all enriched/depleted TFs, limb-associated genes, and transcription factors are also provided. https://doi.org/10.1371/journal.pgen.1005738.s004 (XLSX) S3 Table. Bat in situ probe primers and GenBank numbers for the HoxD gene expression analysis. Bat specific in situ hybridization probes were generated using primers given. Here, we provide the amplicon lengths and GenBank identification numbers for Bat Hoxd10-13 WISH probe sequences. https://doi.org/10.1371/journal.pgen.1005738.s005 (DOCX)

Acknowledgments We would like to thank Alisha Holloway and Dennis Kostka for lending us scripts to calculate significance of BAR and HAR clustering. We would also like to thank Dr. David Ray (Texas Tech University) for the M. lucifugus DNA used in this work.

Author Contributions Conceived and designed the experiments: BMB TF KSP NA. Performed the experiments: BMB TF MKM JEV JZ. Analyzed the data: BMB TF MKM JEV JZ WLE KSP NA. Contributed reagents/materials/analysis tools: ML NI KSP NA. Wrote the paper: BMB TF MKM NI KSP NA.