|Home | About | Journals | Submit | Contact Us | Français|
In the development of the Drosophila embryo, gene expression is directed by the sequence-specific interactions of a large network of protein transcription factors (TFs) and DNA cis-regulatory binding sites. Once the identity of the typically 8–10 bp binding sites for any given TF has been determined by one of several experimental procedures, the sequences can be represented in a position weight matrix (PWM) and used to predict the location of additional TF binding sites elsewhere in the genome. Often, alignments of large (>200 bp) genomic fragments that have been experimentally determined to bind the TF of interest in Chromatin Immunoprecipitation (ChIP) studies are trimmed under the assumption that the majority of the binding sites are located near the center of all the aligned fragments. In this study, ChIP/chip datasets are analyzed using the corresponding PWMs for the well-studied TFs, CAUDAL, HUNCHBACK, KNIRPS and KRUPPEL, to determine the distribution of predicted binding sites. All four TFs are critical regulators of gene expression along the anterio-posterior axis in early Drosophila development. For all four TFs, the ChIP peaks contain multiple binding sites that are broadly distributed across the genomic region represented by the peak, regardless of the prediction stringency criteria used. This result suggests that ChIP peak trimming may result in the exclusion of functional binding sites from subsequent analyses.
The specification of the embryonic body plan in Drosophila melanogaster, as in many other organisms, is accomplished by the regulation of gene expression during development through the activity of cis-regulatory sequences in the genome. Protein transcription factors (TFs) bind to such sequences, regulating expression of target genes . The molecular components of this spatio-temporal network have been extensively characterized (for a review, see ). At the top of the cascade that controls anterio-posterior patterning, translation of spatially localized maternal mRNAs, deposited in the unfertilized egg cell during oogenesis, establishes the early TF gradients in the embryo [3, 4]. In turn, these maternal TFs bind at target embryonic regulatory sequences for gap genes, directing gap TF expression patterns in the developing embryo [5, 6]. Gap TFs further regulate downstream target genes, such as those for pair-rule and homeotic TFs [7, 8]. At each step in the cascade, gene expression patterns are tightly controlled by the binding of TFs to specific clusters of activator and repressor binding sites within embryonic cis-regulatory modules (CRMs). The transcriptional output is mediated by the specific molecular properties of individual CRMs. Whether a given TF acts as an activator or repressor when it binds to a CRM can be context-dependent [9, 10] and the DNA sequences within CRMs may bind TFs with varying affinities [6, 11]. Therefore, the accurate identification of potential TF binding sites in known and putative CRMs is critically important to further our understanding of the molecular control of gene regulation in the developing embryo.
The binding affinity of a TF’s DNA binding domain for any potential TF binding site is a function of the sequence of base pairs that make up that binding site [12, 13]. There is extensive variation in the type and number of DNA binding domains utilized by TFs that can modulate both the sequence specificity and relative affinity for different binding sites [14–17]. As such, the organization of binding sites at genomic regions known to recruit a specific TF has been the subject of many investigations aimed at determining the influence of DNA sequences on TF binding affinity and the computational prediction of binding site locations.
A multitude of bioinformatic tools have been designed to aid in the process of analyzing and predicting genomic fragments that contain TFBSs . The use of position weight matrices (PWMs) has been shown to be very effective in such analyses [18–23]. The datasets required to construct PWM-based models are obtained from protein binding microarrays, yeast and bacterial one-hybrid assays, DNA footprinting assays, and in vitro SELEX experiments [24, 25]. The known binding regions obtained through these experimental procedures are aligned and trimmed to a minimal length of highest conservation (L) and a 4 X L matrix is created containing the frequencies of each nucleotide at each position in the binding site . A PWM is then constructed from this frequency matrix, with the primary underlying assumption that the probabilities for nucleotides at each position are independent of one another. PWMs often also take into consideration the background frequency of the nucleotides in a given genome and are constructed using log likelihoods. The resulting PWM-model can reduce the need to conduct complex in vitro binding assays by predicting the location of unconfirmed binding sites in the genome in silico [19, 22].
In the last decade Chromatin Immunoprecipitation-sequencing (ChIP-seq) and ChIP/chip datasets have significantly increased the ability to identify regions of DNA that bind TFs in vivo at the whole genome level [26–29]. However, a limitation of these studies is the resolution of the binding regions identified. Experimentally identified ChIP peaks are often in the range of 500–1000 bp , while the TF binding sites themselves are usually only 8–10 bp in size [30, 31]. As a result, it is not uncommon for ChIP peaks to be trimmed to the middle 100 bp before using bioinformatic tools to make de novo predictions of TF binding site locations [18, 32]. This approach assumes that the most relevant binding regions for the corresponding TF in vivo are located within this middle 100 bp region of any ChIP peak. This assumption follows from the fact that ChIP peaks are called based on alignments of multiple overlapping immunoprecipitated fragments of randomly sheared chromatin fragments of similar size . The most highly represented regions among the fragments, in theory, therefore should be those that bind the actual TF most often. Thus, assuming a normal distribution of fragment lengths, the most highly represented binding region should be captured in the middle 100 bp. (see Fig. 1 and description in Methods).
In this study we attempt to validate this assumption using PWMs for four well-studied TFs (CAUDAL, HUNCHBACK, KNIRPS and KRUPPEL) involved in the early patterning of the anterio-posterior axis in the Drosophila embryo to predict binding site locations in ChIP peaks from four corresponding ChIP/chip datasets. The results of this analysis indicate that predicted binding sites are in fact broadly distributed throughout the ChIP peaks regardless of the stringency of the predictive thresholds applied, suggesting that ChIP peaks should be left untrimmed to ensure inclusion of all possible binding sites in the dataset.
At any genomic region represented by a ChIP peak, it is possible that the region contains either a single binding site or multiple binding sites for the TF of interest. While the TF binding sites themselves are usually only 8–10 bp in size, for the four Drosophila TFs investigated here: CAUDAL (CAD), HUNCHBACK (HB), KNIRPS (KNI) and KRUPPEL (KR), the mean length of a ChIP peak is ~800 ± 280 bp (Table 1). Trimming ChIP peaks to only consider the center 100 bp, based on the assumptions that a) the genomic fragments recovered for a given peak in a ChIP experiment have a normal distribution around a biologically important binding site(s), and b) that the binding site(s) is likely to be centered around the midpoint of the peak and therefore contained in the central region (Fig. 1), is potentially an efficient approach to help narrow the search for DNA sequences that represent the in vivo binding sites. To examine the distribution of predicted TF binding sites in these ChIP peaks we utilized three different measurements (see Methods for detailed description): directional midpoint distance, position ratio and absolute midpoint distance (Fig. 1).
If we consider the directional midpoint distance for only the top ranked (top 1) binding site in each ChIP peak for each of the four TFs, the average location is very close to the center of the corresponding ChIP peaks, with the mean distance from the midpoint in the narrow range of +1 to +14 bp (Table 2). However, the distribution for the binding sites in each case is surprisingly large (~250 bp) (Fig. 2). Extending this analysis to include the top 2 or 3 scoring sites for each ChIP peak only expands the range of the mean distance from the midpoint to −3 to +14 bp, while the distribution remains broad (Table 2 and Fig. 2).
Previous studies examining the distribution of binding sites for three mammalian TFs in ChIP-seq datasets have identified more restricted patterns of binding site distribution . In each case, if only the single highest scoring predicted site (equivalent to the top 1 binding site in our current study) is considered across a trimmed peak region of 500 bp (−250 bp to +250 bp) then 90% of all binding sites fall in the middle 100 bp (see Figure 7 in ). Comparable analysis of the top 1, 2 or 3 binding sites in each ChIP peak for the four Drosophila TFs in our study reveals a similar general trend to the Wilbanks et al. study , with a mean binding site distance from the middle of the peak of 0 bp (Fig. 3). However, the distribution for the Drosophila TF binding sites is larger, the 25th and 75th percentile values average around −100 bp and +100 bp (Fig. 3), than the less than 50 bp most programs predict for the mammalian TF binding sites . It should be noted that the data in the mammalian studies is generated from ChIP-seq , while the Drosophila studies utilize ChIP-chip . The use of different experimental approaches will likely introduce some systematic differences in the overall size of the identified chromatin fragments, which may impact the ability to accurately compare the results.
It is possible that the “Top n” analyses may have limited biological significance because they do not require the top predicted binding sites to meet any threshold relative to the PWM sequence, which should be indicative of their ability to bind the TF. To address this issue, for each TF, thresholds of varying stringency across a broad range were used in separate analyses (see Methods for details). The use of a range of thresholds is paramount to a robust and unbiased analysis because it has been shown that the performance of a PWM can vary significantly across a range of thresholds . Intriguingly, in each case the results largely remain consistent, with a narrow range for the mean distance to midpoint and a wide distribution (−7 to +49 ± ~270 bp) (Table 2 and Fig. 2). It is worth noting that the pattern is maintained despite the fact that at different thresholds there are vastly different numbers of predicted binding sites (Table 3).
As a control, we also examined the distribution of binding sites for the DORSAL (DL) TF on each of the four ChIP peak datasets. DL is a TF active in the dorsal-ventral patterning network in Drosophila development. Although there is some evidence that DL can bind to a number of genomic regions known to also bind TFs involved in the anterio-posterior patterning of the embryo , binding sites for DL may not be expected to be significantly enriched at ChIP peaks for anterio-posterior TFs when compared to the overall background level in the genome. In almost all cases the total number of DL binding sites predicted is considerably lower than the corresponding experimental TF, particularly at weaker thresholds (Table S2). The exception is KNI, which has an order of magnitude fewer ChIP peaks than any of the other three experimental TFs and a corresponding decrease in the number of predicted KNI binding sites (Table 1 and and3).3). The distribution of the DL binding sites mirrors the pattern observed for the four experimental TFs, but with a larger range at stronger thresholds as a result of the reduced number of sites included in the analysis (Fig. S2 and Table S3).
Analyzing the location of CAD, HB, KNI and KR binding sites using a different measure in the position ratio approach, which allows for normalization of the binding site location relative to the ChIP peak size in each case, reveals a very similar distribution pattern to the directional midpoint approach (Table S1 and Fig. S1). However, an additional feature of the binding site distribution is revealed by measuring the binding site locations using a third approach. Measuring the absolute distance from the midpoint of a ChIP peak reveals a mean of ~220 ± 140 bp (Table 4). Once more this distribution pattern is largely consistent irrespective of which TF or which threshold is analyzed (Fig. 4) and suggests that on average the binding site(s) in any particular ChIP peak are located over 200 bp away from the center of the peak. In combination with the results from the directional and position ratio analyses this indicates that multiple TF binding sites are centered around the midpoint of a given ChIP peak with a broad distribution that may include sites both upstream and downstream of the midpoint.
To investigate this distribution more carefully we analyzed the number of binding sites contained in incremental 100 bp intervals centered on the middle of the ChIP peaks. For each of the four TFs at both weak and strong “Footprint” thresholds the distribution is very consistent. Only ~12.5% of all binding sites are included in the middle 100 bp and only ~25% are in the middle 200 bp (Table 5). Each additional 100 bp interval expanding out from the middle adds ~10% of the total binding sites. As a result, including 800 bp of sequence around the center of the ChIP peaks (an absolute distance of 400 bp in either direction) accounts for ~85% of all the binding sites (Table 5) and the resulting distribution plots are very broad (Fig. 5). The distribution of binding site locations using the strong “Footprint” thresholds are significantly different from a normal distribution (D values for CAD, HB, and KR are 0.04, 0.05, and 0.1 respectively, all with p-values <0.001, see Methods for details), except in the case of KNI (D value of 0.1), which only has 41 predicted binding sites at this threshold (Table 3).
To address the potential biological significance of the observed genome-wide TF binding site distributions we also analyzed the distribution patterns specifically at two well-characterized CRMs from the Drosophila bithorax complex. The 1 kb IAB5 and 1.7 kb IAB8 CRMs are both embryonic enhancers responsible for driving precise spatio-temporal expression of the Abdominal-B gene in early development . Examining the TF binding profile using Berkeley Drosophila Transcription Network Project ChIP/chip tracks  reveals that both CRMs overlap with binding peaks for CAD, HB, KNI and KR, and a number of additional TFs in the anterio-posterior regulatory network (Figure S3). In almost all cases the ChIP peaks are broader than the CRM and many do not align precisely with each other or the center of the defined CRM (Figure S3). In addition, it is clear that genomic regions from within the bithorax complex that do not function as embryonic CRMs  can also recruit these TFs, as evidenced by the presence of separate ChIP peaks outside of the annotated CRMs (Figure S3). This profile supports the observation in earlier studies that many of the key TFs in early Drosophila development bind thousands of active and inactive genomic regions .
To quantitatively assess the distribution pattern for the binding sites in the ChIP peaks associated with the two IAB CRMs we compared the profile to that in the neighboring upstream and downstream genomic flanking regions of equal size. In each case the neighboring flanking regions on either side of the ChIP peak demonstrate no in vivo TF binding (Figure S3) . The analysis is conducted at the all and weak “Footprint” thresholds, as the paucity of predicted binding sites at more robust thresholds prevents meaningful comparison. Despite the fact that there is no obvious pattern to the total number of predicted binding sites in the ChIP peak when compared to the flanking regions at these thresholds, the broad distribution for these sites remains consistent and in the majority of cases (9 out 14) does not show any statistically significant difference (p-values >0.05) (Table S5). In only 1 out of the 14 cases investigated is there a statistically significant narrower distribution in the ChIP peak when compared to the flanking genomic regions (p-value <0.05). In 7 out of the 14 cases, including 3 statistically significant cases (p-value <0.05), the distribution pattern is in fact broadest in the ChIP peak (Table S5). Together these results indicate that the binding site distribution in the ChIP peaks associated with the IAB5 and IAB8 CRMs is in many cases broader than neighboring genomic regions (Table S5) and the genome-wide profile observed if all ChIP peaks are considered (Table 4).
The results from our analysis on four Drosophila developmentally important TFs indicate that the average ChIP peak contains multiple computationally predicted binding sites for a given TF (mean number range = 1.3 – 6.1) if we consider all sites above the defined weak thresholds (Table 3). This raises the question of whether the ChIP peaks may actually represent clusters of binding sites for a TF. If we only consider binding sites above the strong predictive thresholds, the mean number per peak remains greater than one (range = 1.14 – 2.27), indicating the presence of multiple sites in many ChIP peaks and the potential importance of clustering of binding sites. This clustering has been previously recognized as a key feature of regulatory regions in the genomes of a number of different species [36–38], including at cis-regulatory modules in the even-skipped locus that bind the KNI and KR repressors during early Drosophila development [39–41], and may be responsible for mediating homotypic cooperativity of a TF [42, 43] or heterotypic interactions between proteins in ChIP peaks . Indeed, homotypic clustering of binding sites has been shown to be a common feature of cis-regulatory modules active in embryonic development . The conservation of TF binding site organization at a number of cis-regulatory modules in Drosophila and sepsid species [30, 34, 46–48], despite underlying differences in overall genome size , suggests that specific spacing of multiple binding sites in a single ChIP region may in fact be an important evolutionary feature. ChIP datasets are currently not very extensive in insect species other than Drosophila melanogaster, preventing a comprehensive genome-wide investigation of this issue. However, the very precise arrangement of binding sites at the few well-studied cis-regulatory modules, including the IAB5 and IAB8 embryonic enhancers [30, 34], indicates that the local architecture of binding site clusters may be critical to their function.
The results also reveal that the binding sites for these TFs within genomic sequences previously identified in ChIP studies as in vivo binding regions are on average centered in the middle of ChIP peak regions (Fig. 2 and and3),3), but are widely distributed across the peaks (Fig. 5). Comparison to the previously assumed normal distribution for binding sites in ChIP peaks [18, 32] demonstrates that in fact the distribution pattern is significantly broader, with only ~25% of all binding sites contained in the middle 200 bp of an average peak (Table 5). The four Drosophila TFs analyzed in this study therefore demonstrate generally similar binding site profiles to previously studied mammalian TFs in ChIP-seq peaks [33, 49] in their characteristic centering on the middle of the peak. However, given the broad distribution of binding sites it may be wise to exhibit caution when restricting any analysis only to the middle 100 bp, as has been previously applied [18, 32]. Rather, the distribution of binding sites observed suggests it may be prudent to consider an extended sequence in a ChIP peak in future studies if the aim is to retain as much information on TF binding sites as possible. In future studies, the application of next-generation ChIP technologies such as ChIP-exo , where an exonuclease trims ChIP DNA to a precise distance from the crosslinking site, will be of significant value as this enables the mapping of TF binding sites at single base pair resolution.
ChIP/chip peaks were obtained from the Berkeley Drosophila Transcription Network Project (BDTNP) for the TFs; HUNCHBACK (HB), CAUDAL (CAD), KNIRPS (KNI), and KRUPPEL (KR). The datasets  were downloaded from the University of California at Santa Cruz Genome Browser (http://hgdownload.soe.ucsc.edu/goldenPath/dm3/database/). The total number of peaks for each TF is displayed in Table 1. Only those peaks that were equal to or greater than 100 bp in length were included in the analysis because shorter peaks would not be subject to potential information loss when predicting binding sites in the center 100 bp of a peak. The PWMs used for HB, CAD, KNI and KR are as previously published . A PWM for DORSAL (DL)  was utilized as a control.
The spatial arrangement of computationally predicted binding sites was determined using PATSER analysis  with the PWM for each TF on the corresponding set of ChIP peaks. In all cases, PATSER was run with a background sequence A/T content of 0.56, C/G content of 0.44 [51, 53] and, with the exception of the threshold criteria, the default settings.
As each individual ChIP peak should, in theory, contain at least one real binding site for the corresponding TF, the initial analysis was performed such that only the top 1, 2 or 3 scoring PATSER-predicted binding sites for each peak were included (“Top n” analysis). The same “Top n” analyses were also performed using the control DL PWM on each ChIP peak set.
In addition to the “Top n” analyses, binding site score thresholds of different stringencies were used in separate PATSER runs. The ln(p-value) thresholds were obtained by running each PWM on the set of DNA footprinting assay-derived sequence fragments that had originally been used to create the PWM. This produces a range of ln(p-value) scores across all fragments. The 50th percentile of these scores was used as the “weak” threshold in subsequent analyses and the 75th percentile was used as the “strong” threshold, as previously published . In addition to the “weak” and “strong” thresholds, the 0th percentile of the scores was used as the “all” threshold, as it includes all sequences that score at or above the lowest-scoring footprinted binding site and the 100th percentile of the scores is used as the “strongest” threshold (“Footprint” threshold analysis). A second set of analogous thresholds were derived using the same approach, but percentiles were defined from the scores of the binding sites PATSER output for each “top 1” analysis on the ChIP peaks (“ChIP” threshold analysis).
Functional thresholds for the control DL PWM were obtained by running an analysis on the 330 bp rhomboid (rho) neuroectodermal enhancer (NEE). The rho NEE has been shown to contain four binding sites for DL and their locations have been accurately determined via DNA footprinting . The lowest-scoring of these binding sites was used as the “weak Footprint” threshold, while the highest scoring was used as the “strong Footprint” threshold. Using percentile cutoffs for thresholds in this case would have been superfluous as there are only four binding sites for DL; however, these thresholds can be compared to those derived from the experimental PWMs  because both sets of thresholds were derived from footprint data. A second set of “ChIP” thresholds was derived for DL from the “top 1” analysis using the same percentile cutoff approach employed for the experimental PWMs.
For each TF, three different measures of computationally predicted binding site location in each ChIP peak were used to assess the distribution of binding sites (see Fig. 1 for a graphical representation). A “directional distance” (in bp) for each binding site position from the midpoint of each peak was calculated for each of the analyses described above. An “absolute distance” was calculated in the same way but did not differentiate between which side of the midpoint the binding site was located. A “position ratio” value was calculated by dividing the binding site position in the ChIP peak by the total length of the peak. In this case, the ratio values 0 and 1 correspond to the ends of the peak and a value of 0.5 represents the center of the peak. Boxplots of the data were generated in MATLAB with default settings for each TF.
The average number of binding sites per ChIP peak was calculated for each “Footprint” threshold and “ChIP” threshold analysis. In addition, the total number and percentage of weak and strong “Footprint” threshold binding sites in different regions centered on the middle of the ChIP peaks was also calculated. The regions range from the middle 100 bp to the middle 800 bp, in increments of 100 bp.
To compare the distribution of predicted binding site locations from the midpoint at the strong “Footprint” threshold to the expected normal distributions, we performed a Kolmogorov-Smirnov test. For each TF, the position of predicted binding sites was tested against the null hypothesis that the positions follow a normal distribution with a mean of 0 and a standard deviation that results in the same percentage of binding sites falling into the middle 800 bp as is predicted using the strong “Footprint” threshold. The statistical analysis was done at a 5% significance level.
To compare the distribution of predicted binding site locations in the IAB5 and IAB8 CRMs to the distributions in the neighboring upstream or downstream genomic flanking regions of equal size, we performed two-sample Kolmogorov-Smirnov tests. For each TF, the position of predicted binding sites within the CRM-associated ChIP peak region was tested first against the null hypothesis that the position of predicted binding sites within the neighboring upstream region came from the same distribution and secondly against the null hypothesis that the position of predicted binding sites within the neighboring downstream region came from the same distribution. The statistical analysis was done at a 5% significance level.
For each of the three different measures of computationally predicted binding site locations, a box-and-whisker plot was generated in MATLAB. In each plot, the central mark in the box represents the median value, while the bottom and top of the box correspond to the 25th and 75th percentiles of the values. The whiskers extend to include all measurement values not considered outliers. Values are defined as outliers if they are greater than y + 1.5 (y – x) or less than x − 1.5 (y – x), where x and y are the 25th and 75th percentile values, respectively.
Boxplots for each TF, with the y-axis representing the position ratio of binding sites in the ChIP peak, for “Top n”, “Footprint”, and “ChIP” analyses. The red line indicates the median position ratio and the middle 50% of binding sites are contained within the upper and lower bounds of each box. Whiskers extend to include all positions not marked as outliers.
The mean and standard deviation of the position ratio for binding sites, and the ln(p-value) threshold used in the “Top n”, “Footprint” threshold and “ChIP” threshold analyses are listed. A ‘−’ corresponds to entries that are not applicable.
The total number of predicted DORSAL (DL) binding sites and the mean number of sites per ChIP peak are listed for each “Footprint” and “ChIP” threshold analysis.
The mean and standard deviation of the directional midpoint distances for DORSAL (DL) binding sites, and the ln(p-value) threshold used in the “Top n”, “Footprint” threshold and “ChIP” threshold analyses are listed. A ‘−’ corresponds to entries that are not applicable.
The number and percent of total DORSAL (DL) binding sites obtained in the weak and strong “Footprint” analyses in ChIP peak regions ranging from the middle 100 bp to the middle 800 bp in increments of 100 bp is shown for each of the four corresponding anterio-posterior TF peak regions. The number of binding sites over the entire length of the ChIP peaks (all) is also listed.
Boxplots for DORSAL (DL) analysis on each set of ChIP peaks, with the y-axis representing the directional distance of binding sites from the ChIP peak midpoint, for “Top n”, “Footprint”, and “ChIP” analyses. The red line indicates the median directional midpoint distance and the middle 50% of binding sites are contained within the upper and lower bounds of each box. Whiskers extend to include all positions not marked as outliers.
The location of the characterized CRMs in the bithorax complex are shown as a custom track in the UCSC Genome Browser (http://genome.ucsc.edu/; ). The Berkeley Drosophila Transcription Network Project ChIP/chip track  shows the location of verified binding for selected anterio-posterior gap/terminal (green) and pair-rule (yellow) transcription factors in stage 4–5 embryos (1% false discovery rate). The center 100 bp of each CRM is indicated (red rectangle).
The total number of predicted binding sites in the corresponding ChIP peak region mapping to each CRM and the ln(p-value) threshold used in the “Footprint” threshold analyses are listed. Color coding for the number of predicted sites indicates that the ChIP peak contains more sites than either of the neighboring upstream or downstream genomic flanking regions of equal size (green), less sites than either flanking region (red) or an intermediate number (orange). The mean and standard deviation of the directional midpoint distance and normalized absolute midpoint distances for binding sites are also shown. A ‘−’ corresponds to entries that are not applicable. Color coding for the normalized absolute position ratio indicates that the value is higher (representing a broader distribution) in the ChIP peak than either of the neighboring genomic flanking regions of equal size (green), lower than either flanking region (red) or an intermediate value (orange). Statistical significance of the differences in the binding site distribution patterns between the ChIP peak and the neighboring genomic flanking regions of equal size was assessed using two-sample K-S tests using the directional midpoint distance values. The resulting normalized absolute position ratios are indicated as significantly different to both (**) or one (*) of the flanking regions as appropriate.
The research in this paper was funded by National Institutes of Health (GM090167) and National Science Foundation (IOS-0845103) grants to RAD, and a National Institutes of Health (GM110571) grant to RAD and JMD. K.P.P. was funded as a Fellow from a National Science Foundation Interdisciplinary Training Award to Amherst College (DBI1129152).
The authors report no competing interests.
Authors’ contributionsK.P.P. participated in the design of the study, performed the analysis, and drafted the manuscript. J.M.D. jointly conceived the study, participated in the design of the study and coordination, and drafted the manuscript. R.A.D. jointly conceived the study, participated in its design and coordination, and drafted the manuscript. All authors read and approved the final manuscript.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.