<title>Introduction</title> <p>After dispersal from African, modern human migrated to the rest of the world and adapted rapidly to a variety of environmental challenges, such as climate, food supply, etc. Many phenotypes may undergo apparent adaptation ##REF##16778047##[1]##. Skeletal system has been observed to evolve rapidly during the past 10,000 years among human populations during the time of rapid population growth ##UREF##0##[2]##. However, unlike other beneficial traits, e.g. speech, cognitive ability, diverse skin and others ##REF##16778047##[1]##, virtually no single gene has been identified under a major selection event accounting for the rapid evolution of skeletal system of human.</p> <p>Polydactyly, characterized by the addition or part of a digit, is one of the most common congenital limb malformations which are relatively common human abnormalities occurring at an incidence of one per ∼500–1000 live births ##REF##17300748##[3]##. Polydactyly occurs in many phenotypes, like preaxial polydactyly II (PPD II), triphalangeal thumb-polysyndactyly syndrome (TPT-PS) and isolated triphalangeal thumb (OMIM 174500) ##REF##17300748##[3]##–##UREF##1##[7]##. <italic>LMBR1</italic> gene, which contains 17 exons spanning approximately 200 kb of genomic DNA, encodes a 490-amino acid protein containing 9 putative transmembrane and one coiled-coil domains ##REF##10945466##[8]##. Level of <italic>LMBR1</italic> activity had been associated with the number of digits across vertebrates ##REF##11606546##[4]##. The crucial functional element of <italic>LMBR1</italic> is located within the intron 5, which serves as the long-range regulatory element of the adjacent <italic>SHH</italic> gene, a key development gene in the nervous system, skeletal system and others. Disruption of this intron, leading to dysregulation of <italic>SHH</italic>, can cause all kinds of above mentioned polydactyly ##REF##17300748##[3]##–##UREF##1##[7]##. Phylogenetic analysis also indicated conservation of the intron 5 region in teleost fishes and throughout the tetrapod lineage ##REF##14727139##[9]##.</p> <p>Considering the profound role of <italic>LMBR1</italic> gene, particularly the intron 5, in the development of limb and skeletal system, we checked the evolutionary pattern by sequencing a ∼10 kbp region in the intron 5 in 41 East Asian individuals. Tajima's D value is significantly higher than neutrality as after considering human population growth. Additionally, the advent of large-scale genome and polymorphism data in human population supports specific selection effect during human evolution. Extensive linkage disequilibrium and lower genetic differentiation were found in this region in East Asian and European populations. We concluded that balancing selection at the region occurred during the evolution of modern human.</p>
<title>Materials and Methods</title> <title>Samples and sequencing</title> <p>41 unrelated East Asian individuals (19 Han Chinese, 15 minorities of China, five Thais, 1 Filipino and 1 Lao) were chosen for sequencing. Ethical approval for this study was provided by the Ethics Committee of Kunming Institute of Zoology, Chinese Academy of Sciences, and all participants provided written informed consent. The first ten kbps of <italic>LMBR1</italic> intron 5 was amplified by LA-PCR method with two pair primers (<named-content content-type="gene">5′-AGAAAGGAGGTCATTGTAG-3′</named-content> as first sense primer and <named-content content-type="gene">5′-AGATTGAGGTCCAGGTAT-3′</named-content> as first antisense primer; <named-content content-type="gene">5′-CGTATGGGAACTCAGAAA-3′</named-content> as second sense primer and <named-content content-type="gene">5′-ACGCAAGCCAAATAAGAC-3′</named-content> as second antisense primer), and sequenced by ABI PRISM 3730xl DNA analyser (Applied Biosystems) with ABI BigDye Terminator Cycle Sequencing Kit, Version 3.1 (Applied Biosystems, Foster City, California, USA). The thermal cycling condition of two pair primers both are: 95°C, 4 min; 94°C, 1 min, 57°C, 5 min, 72°C, 5 min, 30 cycles; 72°C, 10 min. The resulting sequences were analyzed by the DNASTAR software (DNASTAR). Detailed information on sequencing primer sequences is available on request.</p> <title>Statistical analysis on the sequenced intron 5 region</title> <p>The haplotypes were constructed by PHASE program ##REF##11254454##[10]##, ##REF##15700229##[11]##, and the network was constructed by using median joining algorithm ##REF##10331250##[24]## implemented in Network. The average number of pairwise difference (π), Watterson's estimator (θw) ##REF##1145509##[25]## and haplotype diversity were calculated. Tajima's D (1989) ##REF##2513255##[26]## was used to test the evolutionary pattern by Arlequin program ##UREF##3##[27]##. Allowing for human population growth, we applied the Tajima's D test on the background of different magnitudes of population growth and the growth beginning at different times simulated by the algorithm of Rogers (1995) ##REF##28565146##[15]##. Pairwise mismatch analysis was performed by DnaSP ##REF##10089204##[28]## under population growth assumptions with initial theta as 3.442, final theta as 1000, and final tau as 2.267. LD extent was analyzed by R<sup>2</sup> of all pairwise comparisons between SNPs, and the significances were identified with χ<sup>2</sup> tests by using a Bonferroni correction for multiple testing. The nucleotide diversity (π) values of previous reported balancing selection genes were obtained from ##REF##15996168##[29]## (<italic>LDLR</italic>), ##REF##15889130##[30]## (<italic>HAVCR1</italic>), ##REF##15361935##[20]## (<italic>ABO</italic>, <italic>IL10R</italic>B, <italic>IL1A</italic>, and <italic>ACE2</italic>), ##REF##12149450##[21]## (5′ <italic>CCR5</italic>), ##REF##15306844##[31]## (<italic>MBL2</italic>), ##REF##16252244##[32]## (<italic>CD209L</italic>), ##REF##15845028##[33]## (<italic>C6</italic>), ##REF##14997422##[14]## (<italic>PTC</italic>), ##REF##17227474##[34]## (<italic>FSHB</italic>), ##UREF##4##[35]## (<italic>FMO3</italic>), ##REF##12378426##[36]## (<italic>G6PD</italic>), ##REF##14569463##[37]## (<italic>FUT2</italic>).</p> <title>LD analysis and Fst comparison based on the HapMap data</title> <p>SNPs of <italic>LMBR1</italic> location in chr7:155920–156290 kbp (NCBI35) were chosen from HapMap with the criteria: minor allele frequency ≥10% and consistent with Hardy-Weinberg equilibrium in 0.01 level. LD measures between pairs of SNPs were quantified using statistic D' ##REF##14569463##[37]##, which were calculated by Haploview program ##REF##15297300##[39]##. The results of pairwise D' were visualized by the GOLD program ##REF##10842743##[40]##. SNPs in the <italic>LMBR1</italic> gene region were used to calculate Fst values between Caucasians, Africans, and East Asians. Total 113 SNPs in the <italic>LMBR1</italic> region fulfilled the criteria were used to calculate the Fst values among human populations as described in ##REF##12466284##[17]##, ##REF##28563791##[41]##.</p>
<title>Results</title> <title>Genetic variation of the sequenced region in the intron 5 in the East Asians</title> <p>We sequenced one 9256 bp region in the 5th intron of <italic>LMBR1</italic> gene in 41 East Asian individuals, and identified 21 SNPs, 8 of which had the minor allele frequencies &gt;0.4. The sequence for each individual was submitted to GenBank under accession numbers EU880543-EU880583. One SNP was difficult to be confirmed in about one fourth individuals for sequencing technology difficulties, the SNP and 300 bp sequence around it were excluded from analysis. The other total 20 SNPs were used to construct the haplotypes by the PHASE program ##REF##11254454##[10]##, ##REF##15700229##[11]##, and 13 haplotypes were obtained (##FIG##0##Figure 1A##). All SNPs were not deviated from Hardy-Weinberg equilibrium.</p> <p>The nucleotide diversity π is 6.4×10<sup>−4</sup>, lower than the human genome average (∼7.5×10<sup>−4</sup>) ##REF##11237013##[12]##, and the Watterson's estimator θw is 4.49×10<sup>−4</sup>. The haplotype diversity is 0.828. There are two major haplotypes found at high frequencies, and the network analysis indicates that they are separated by relative long branch length (##FIG##0##Figure 1B##). The “twin-peak” phenomenon is observed in the mismatch distribution of pairwise differences between the haplotypes (##FIG##0##Figure 1C##). Tajima's D value is 1.25 (p&gt;0.1), but the value is significantly higher at the empirical 95% level (by one sided p-value) in an empirical distribution available from a study of 313 genes ##REF##11452081##[13]##. These data are suggestive of the direction of balancing selection.</p> <p>Evidences have suggested that human populations have grown dramatically, which have strong effects on the genetic diversity. The assumption of constant population size is highly conservative in the detection of balancing selection, and will cause type II statistical error (i.e. a failure to reject the null hypothesis of neutrality when it is false) ##REF##14997422##[14]##. Considering the population size change in human populations, we performed 250 different tests of Tajima's D using the algorithm of Rogers ##REF##28565146##[15]## across different magnitudes of population growth (from 1-fold growth [i.e., no growth] to 250-fold growth), with the growth beginning at different times (0 years ago to 250,000 years ago). The hypothesis of neutrality is rejected by Tajima's D under those models assuming magnitudes of growth greater than 1.3, from an ancient effective population size of 10,000, beginning more than about 65,000 years ago. The sequence is large enough that recombination may have occurred, and it would make Tajima's D test conservative.</p> <title>Age estimation</title> <p>The time of the most recent common ancestor (T<sub>MRCA</sub>) is calculated as T<sub>HC</sub>*D<sub>H</sub>/D<sub>HC</sub>. T<sub>HC</sub> is the time of divergence between human and chimpanzee, D<sub>HC</sub> is the divergence between human and chimpanzee sequences (from NCBI) and D<sub>H</sub> is the average difference of human haplotypes. The values are 6×10<sup>6</sup> years, 95.00, and 5.71 respectively, and the T<sub>MRCA</sub> is 3.61×10<sup>5</sup> years, much longer than 65,000 years.</p> <title>Extensive linkage disequilibrium of the intron 5 in the East Asian and European populations</title> <p>Extended linkage disequilibrium is usually observed in the selected region because recombination does not have enough time to break it down ##REF##12397357##[16]##. Large scale polymorphism data are facilitating the studies of evolutionary patterns in human genome. We examined the characteristics of <italic>LMBR1</italic> by the HapMap data. All pairwise D' measures among these HapMap SNPs were estimated, and the graphical representation of LD level is illustrated in ##FIG##1##Figure 2##, which demonstrates strong LD of the intron 5 region in the East Asian and European populations (##FIG##1##Figure 2##). We also estimate the LD extent by analyzing the R<sup>2</sup> of all pairwise comparisons between the 20 SNPs, and found 56 significant pairs at 5% level (##FIG##0##Figure 1D##).</p> <title>Lower genetic differentiation among human populations</title> <p>In general, positive directional selection promotes the divergence among populations, but balancing selection will decrease the differentiation of selected loci compared with neutrality. The average Fst values of total 113 SNPs in the <italic>LMBR1</italic> gene were 0.036 (East Asian vs European), 0.081 (East Asian vs African), and 0.064 (European vs African), significantly lower than the average values in ##REF##12466284##[17]##, 0.098, 0.128 and 0.102 respectively by t test (p = 9.36E-44, 1.57E-07, 2.08E-06 respectively with degree of freedom 112; p = 0.00098, 0.016, p = 0.02 respectively with degree of freedom 2; and p = 0.014, 0.058, 0.065 respectively with degree of freedom 1). Such weak genetic differentiation indicates balancing selection might have occurred in this region rather than population subdivision.</p>
<title>Discussion</title> <p>Evidences have indicated that the size of human population increased in the Upper Pleistocene ##REF##28565146##[15]##. Populations that have grown are expected to have an excess of low-frequency alleles and thus low pairwise difference between sequences, which will lead to the reduction of common statistics used to detect from neutrality, e.g. Tajima's D, ##REF##12196407##[18]##. Therefore, it is inappropriate to detect natural selection, e.g. conservative to detect balancing selection, under the model of constant population size ##REF##14997422##[14]##. In this study, we identified significant deviation of Tajima's D from neutrality under models incorporating different human population growth parameters. Another two pieces of evidences, strong linkage disequilibrium and lower genetic differentiation among human ethnical populations also support the existence of a balancing selection, because population subdivision, another competing hypothesis, could also lead to significantly high Tajima's D for divergent haplotypes existing in different geographical regions ##REF##14997422##[14]##, ##REF##1552836##[19]##. However, the nucleotide diversity is low in the region, which is not usually observed in the genes under balancing selection (##FIG##2##Figure 3##). Perhaps, the intron 5 region is highly conservative during evolution for its essential function and does not allow accumulation of new mutations. For example, the nucleotide diversity of <italic>ACE2</italic> gene, subject to balancing selection, is even lower than that found in this intron ##REF##15361935##[20]##.</p> <p>Cis-regulating modules have been known contributing profoundly to the genome evolution as they are key regulator of gene expression level. Up to now, cis-regulatory regions of many genes showed evidence of positive and/or balancing selection. For example, the cis-regulator of <italic>CCR5</italic>
##REF##12149450##[21]## was involved in a strong balancing selection, whereas, the encoding sequence of the gene was subject to positive selection. The intron 5 and the downstream <italic>SHH</italic> gene appear evolving in a similar manner as <italic>CCR5</italic>, in that evidence of positive selection has been found in <italic>SHH</italic> gene in primate, particularly lineage leading to human, and human population ##REF##16687440##[22]##.</p> <p>It is appealing that testing for selection should be based on a strong biological working hypothesis ##UREF##2##[23]##. The lower genetic differentiation of the loci among human populations indicated that the balancing selection had occurred before the divergence between the main modern human subpopulations. Accordingly, we proposed that polydactyly has maintained for a long time during the evolution of human. However, considering the strangeness appearance of the phenotype, it may be a slightly side effect and disadvantage during the adaptive evolution of limb and skeletal system, like limb size, hand bone strength, finger flexibility or others. Presumably, the variants maintained by balancing selection may be associated additional intermediate phenotypes other than polydactyly. In addition, the balancing selection on the region probably also plays a role in maintaining the diversity of skeletal system, e.g. sizes, among different populations and different individuals. Although we showed the evidence of balancing selection in this region, the mechanism under which balancing selection occurred at the intron 5 of <italic>LMBR1</italic> gene is unclear. It may become more apparent when more insight into function of the intron 5 is available by the future functional studies.</p>
<p>Conceived and designed the experiments: FH DDW QPK YPZ. Performed the experiments: FH. Analyzed the data: FH DDW. Wrote the paper: FH DDW YPZ.</p> <title>Background</title> <p>The intron 5 of gene <italic>LMBR1</italic> is the cis-acting regulatory module for the sonic hedgehog (<italic>SHH</italic>) gene. Mutation in this non-coding region is associated with preaxial polydactyly, and may play crucial roles in the evolution of limb and skeletal system.</p> <title>Methodology/Principal Findings</title> <p>We sequenced a region of the <italic>LMBR1</italic> gene intron 5 in East Asian human population, and found a significant deviation of Tajima's D statistics from neutrality taking human population growth into account. Data from HapMap also demonstrated extended linkage disequilibrium in the region in East Asian and European population, and significantly low degree of genetic differentiation among human populations.</p> <title>Conclusion/Significance</title> <p>We proposed that the intron 5 of <italic>LMBR1</italic> was presumably subject to balancing selection during the evolution of modern human.</p>
<p>We are grateful to Dr. Stephen Wooding from the Department of Human Genetics, University of Utah School of Medicine for helping stimulating and analyzing the data, Dr. Yong-Gang Yao for helpful suggestions and comments, and Dr. Nelson Leung-sang Tang for reading and polishing the manuscript. We also thank Li Hui and Xie Liang for technical helping; and Gou S.-K., Wu S.-F., and Zhu C.-L. for technical assistance.</p>
<fig id="pone-0002948-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002948.g001</object-id><label>Figure 1</label><caption><title>Genetic variation analyses in the <italic>LMBR1</italic> intron 5 among 41 East Asian individuals.</title><p>A: The 13 haplotypes constructed by PHASE program, and the right-most column shows the number of each haplotypes among 41 subjects. B: Median joining network of haplotypes. Each circle represents a haplotype indicated in ##FIG##0##Figure 1A##, and the size of the circle is the relative frequency. Beside the branches are labels of the SNPs in the haplotypes counted from left to right. C: Graph of pairwise differences between the haplotypes. The dash line represent the observed sequence pairwise difference, and the real line represent the expected distribution of pairwise difference simulated by DnaSP under population growth with initial theta as 3.442, final theta as 1000, and final tau as 2.267. The “twin-peak” of observed mismatch distribution is suggestive of balancing selection. D: LD extent analyzed by R<sup>2</sup> of all pairwise comparisons between the 20 SNPs. The shadows indicate significant pairwise comparison identified with χ<sup>2</sup> tests by using a Bonferroni correction for multiple testing.</p></caption></fig> <fig id="pone-0002948-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002948.g002</object-id><label>Figure 2</label><caption><title>Linkage disequilibrium pattern of chr7: 155920–156290 kbp (NCBI35) region in the East Asian, European and African populations based on the HapMap Data.</title><p>
<italic>LMBR1</italic> gene and the intron 5 are showed.</p></caption></fig> <fig id="pone-0002948-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pone.0002948.g003</object-id><label>Figure 3</label><caption><title>Nucleotide diversity (π) of previous reported balancing selection genes and the <italic>LMBR1</italic> intron 5 studied here.</title><p>It shows that <italic>LMBR1</italic> intron 5 had low π among these genes with documented evidence of balancing selection. The data are from ##REF##15996168##[29]## (<italic>LDLR</italic>), ##REF##15889130##[30]## (<italic>HAVCR1</italic>), ##REF##12196407##[18]## (<italic>ABO</italic>, <italic>IL10R</italic>B, <italic>IL1A</italic>, and <italic>ACE2</italic>), ##REF##1552836##[19]## (5′ <italic>CCR5</italic>), ##REF##15306844##[31]## (<italic>MBL2</italic>), ##REF##16252244##[32]## (<italic>CD209L</italic>), ##REF##15845028##[33]## (<italic>C6</italic>), ##REF##11452081##[13]## (<italic>PTC</italic>), ##REF##17227474##[34]## (<italic>FSHB</italic>), ##UREF##4##[35]## (<italic>FMO3</italic>), ##REF##12378426##[36]## (<italic>G6PD</italic>), ##REF##14569463##[37]## (<italic>FUT2</italic>).</p></caption></fig>
<fn-group><fn fn-type="COI-statement"><p><bold>Competing Interests: </bold>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p><bold>Funding: </bold>This work was supported by the Key State Research Program of China (2006CB943900), Science and Technology Committee of Yunnan Province and the National Natural Science Foundation of China.</p></fn></fn-group>
<graphic xlink:href="pone.0002948.g001"/> <graphic xlink:href="pone.0002948.g002"/> <graphic xlink:href="pone.0002948.g003"/>
[{"label": ["2"], "element-citation": ["\n"], "surname": ["Armelagos", "Harper"], "given-names": ["GJ", "KN"], "year": ["2005"], "article-title": ["Genomics at the origins of agriculture, part one."], "source": ["Evol Anthropol"], "volume": ["14"], "fpage": ["68"], "lpage": ["77"]}, {"label": ["7"], "element-citation": ["\n"], "surname": ["Gurnett", "Bowcock", "Dietz", "Morcuende", "Murray"], "given-names": ["CA", "AM", "FR", "JA", "JC"], "year": ["2007"], "article-title": ["Two novel point mutations in the long-range SHH enhancer in three families with triphalangeal thumb and preaxial polydactyly."], "source": ["Am J Med Genet A"], "volume": ["143"], "fpage": ["27"], "lpage": ["32"]}, {"label": ["23"], "element-citation": ["\n"], "surname": ["Kreitman", "Di Rienzo"], "given-names": ["M", "A"], "year": ["2004"], "article-title": ["Balancing claims for balancing selection."], "source": ["Trend Genet"], "volume": ["20"], "fpage": ["300"], "lpage": ["304"]}, {"label": ["27"], "element-citation": ["\n"], "surname": ["Excoffier", "Laval", "Schneider"], "given-names": ["L", "G", "S"], "year": ["2005"], "article-title": ["Arlequin version 3.0: An integrated software package for population genetics data analysis."], "source": ["Evolutionary Bioinformatics Online"], "volume": ["1"], "fpage": ["47"], "lpage": ["50"]}, {"label": ["35"], "element-citation": ["\n"], "surname": ["Allerston", "Shimizu", "Fujieda", "Shephard", "Yamazaki"], "given-names": ["CK", "M", "M", "EA", "H"], "year": ["2007"], "article-title": ["Molecular evolution and balancing selection in the flavin-containing monooxygenase 3 gene (FMO3)."], "source": ["Pharmacogenet Genom"], "volume": ["17"], "fpage": ["827"], "lpage": ["839"]}]
PLoS One. 2008 Aug 13; 3(8):e2948