A single-nucleotide mutation within the TBX3 enhancer increased body size in Chinese horses

Summary Chinese ponies are endemic to the mountainous areas of southwestern China and were first reported in the archaeological record at the Royal Tomb of Zhongshan King, Mancheng, dated to approximately ∼2,100 YBP.1 Previous work has started uncovering the genetic basis of size variation in western ponies and horses, revealing a limited number of loci, including HMGA2,2LCORL/NCAPG,3ZFAT, and LASP1.4,5 Whether the same genetic pathways also drive the small body size of Chinese ponies, which show striking anatomical differences to Shetland ponies,6 remains unclear.2,7 To test this, we combined whole-genome sequences of 187 horses across China. Statistical analyses revealed top association between genetic variation at the T-box transcription factor 3 (TBX3) and the body size. Fine-scale analysis across an extended population of 189 ponies and 574 horses narrowed down the association to one A/G SNP at an enhancer region upstream of the TBX3 (ECA8:20,644,555, p = 2.34e−39). Luciferase assays confirmed the single-nucleotide G mutation upregulating TBX3 expression, and enhancer-knockout mice exhibited shorter limbs than wild-type littermates (p < 0.01). Re-analysis of ancient DNA data showed that the G allele, which is most frequent in modern horses, first occurred some ∼2,300 years ago and rose in frequency since. This supports selection for larger size in Asia from approximately the beginning of the Chinese Empire. Overall, this study characterized the causal regulatory mutation underlying small body size in Chinese ponies and revealed size as one of the main selection targets of past Chinese breeders.

In brief Liu et al. discover a regulatory mutation in TBX3 driving size variation in Chinese horses. Reporter assays and mice models confirm that it can increase TBX3 transcription and the limb length. They find this mutation first occurred 2,300 years ago and rose in frequency since, revealing size as one of the main selection targets of past Chinese breeders.

Genome sequencing and phylogenetic analysis
We sequenced the genomes of 48 ponies and nine horses at an average depth of coverage of $11.13 and combined them with 130 genomes previously characterized 7,8 (Data S1A; STAR Methods). Our dataset includes 52 Chinese ponies from four breeds, 100 Chinese native horses from 12 breeds, 23 imported Falabella ponies, seven Thoroughbred horses, and five Przewalski's horses. Stringent quality filtering identified a total of 20,428,318 SNPs, 65.1% of which were within introns and 20.2% within intergenic regions (Data S1B; STAR Methods). Principal component analysis (PCA) and phylogenetic tree reconstruction revealed that all Chinese ponies clustered together in a single, monophyletic group (hereafter referred to as Southern ponies, SC; Figures S1A and S1B). Two other main clusters consisting of Northern horses (NC) and Qinghai-Tibetan horses (QT) were identified, in line with previous analyses. 8 Pony breeds showed genome-wide heterozygosity levels generally higher than those found in horse breeds (Data S1C; STAR Methods). Their genomes were also characterized by shorter runs of homozygosity fragments (ROHs) ( Figure S1D), possibly reflecting relaxed selective pressures in Chinese ponies relative to horses.

Selection signatures among horses and ponies
We scanned the genomes of 152 Chinese native ponies and horses for signatures of positive selection that may underpin their different statures. To achieve this, we calculated three complementary statistics along the horse reference genome using 100-kb-long sliding windows (step size = 15 kb). The first statistic was the pairwise genetic differentiation F ST index for identifying those genomic regions maximizing allelic frequency differences between ponies and horses. The second statistic measured relative differences in nucleotide diversity between the two groups (qp(pony/horse)), with high or (low) qp values indicative of positive selection in Chinese horses or ponies, respectively. Selection candidates were also identified through the calculation of crosspopulation extended haplotype homozygosity (XP-EHH). In order to limit false positive identifications, we considered the top 1% windows from all three scans, which provided a total of 153 windows encompassing 4.66 Mb and 69 protein-coding genes (Figure 1A; Data S2A-S2D; STAR Methods). These selection candidates showed limited genetic diversity in horses, but not in ponies, suggesting positive selection for horse-related phenotypes ( Figure 1A). Most of the protein-coding genes were significantly enriched for functional Gene Ontology categories related to the anterior/posterior pattern specification (p = 7.77eÀ7, following g:SCS ''Set Counts and Sizes'' correction 9 ), the embryonic skeletal system development (p = 3.07eÀ6), embryonic skeletal system morphogenesis (p = 1.27eÀ5), and pattern specification process (adjusted p = 3.59eÀ5) (Data S2E). Importantly, all the four outlier genomic regions providing the strongest signatures of selection contained TBX3 (ECA8:20,640,001-20,740,  000 bp, ECA8:20,655,001-20,755,000 bp, ECA8:20,625,001-20,725,000 bp, and ECA8:20,610,001-20,710,000 bp) ( Figure 1A; Data S2D) and were further refined into an 80-kb region (ECA8: 20,640,000-20,720,000 bp) repeating calculations in 20-kb sliding windows ( Figure 1B; Data S2F-S2H).
TBX3 is a T-box gene involved in the development of the anterior/posterior axis, 10 which was previously reported to undergo selection in Debao ponies. 2 This indicated that TBX3 could play an important role in driving body size differences between Chinese ponies and horses. Interestingly, the strongest selection signature identified in Falabella ponies, which represent the smallest miniature horse breed on the planet, 11 was located around another gene, HMGA2 ( Figure S2A), which was previously reported to be associated with short stature in Shetland ponies. 12 Including other western pony breeds to the selection analysis confirmed previously reported selection candidates, such as FGF12, DLX3, and PDK2, 13,14 but not TBX3 (Data S2I-S2K). This suggested a potentially different genetic basis for size-related phenotypes in Chinese and western ponies (Data S2J and S2K).
In addition to showing selection signatures in Debao ponies, 2 TBX3 was previously described to cause asymmetric deposition of skin pigments characteristic of Dun phenotypes when present in the ancestral, non-deleted allelic state. 15 In order to check whether the selection signature identified above could in fact target Dun phenotypes, we screened the presence of not deleted (Dun and non-dun1) and deleted (non-dun2) alleles in all 187 genomes present in our panel ( Figure S2B; Data S2L). We found Dun alleles in five Przewalski's horses (PrZ) and one Chinese Hequ horse, but not in any of the remaining 99 Chinese horses, 52 Chinese ponies, 23 Falabella ponies, and 7 Thoroughbred horses. This ruled out Dun phenotypes driving the selection signatures identified.   Report Genome-wide association study We next carried out a genome-wide association analysis for the wither height of all 152 Chinese samples (STAR Methods). This analysis accounted for the underlying population structure by adjustment of above mentioned PCA covariates and genetic relatedness among all pairs of individuals (STAR Methods). It revealed two statistically significant regions (p < 6eÀ8), encompassing 82 SNPs ( Figures 1C and 1D; Data S2M). The first region spanned the TBX3/TBX5 locus (ECA8:20,566,781-20,865,109), which also provided the top-selection candidate. The second region comprised the HMGA2 locus ECA6:82,576,739-82, 701,867. This result indicated that TBX3, TBX5, and/or HMGA2 may provide the genetic basis for the short stature of Chinese ponies.
In order to further identify the main underlying genetic driver, we investigated patterns of differential placental gene expression in four Chinese YiLi horses and three NiQi ponies (Figure S2C). We focused on the placenta as this tissue is essential for fetal growth, 16 and the placental surface structure is known to be associated with the neonatal body size. 17 Overall, we identified a total of 620 differentially expressed genes in the transcriptomes of YiLi horses and NiQi ponies (Data S2N). Interestingly, TBX3 was highly expressed in horse placenta. In contrast, none of the other two candidate genes (TBX5 and HMGA2) showed substantial expression levels (fragments per kilobase of exon per million fragments mapped, FPKM < 0.3). TBX3 placental expression was also 4.3-fold higher in YiLi horses than NiQi ponies (p = 0.04; Figure 1E), indicating positive association of TBX3 expression levels and adult wither height. This finding was in line with the absence of missense mutations in the TBX3 region examined, and the lower nucleotide diversity and negative Tajima's D statistics found in this region in horses ( Figure 1B). Altogether, this suggested that the genetic variation present within the TBX3 region could regulate gene expression and potentially drive size differences in Chinese horses and ponies. We thus decided to focus the following experimental validation efforts on the TBX3 locus.
Identifying TBX3 causative variant As the TBX3 region showing selection signatures and associated with wither height encompassed approximately $80 kb, we carried out a number of analyses to narrow down the list of possible causative variants (STAR Methods). We first explored the patterns of linkage disequilibrium (LD) in the region using a haplotype block analysis (Data S2O), which revealed nine blocks (H1-H9, haplotype p < 1eÀ5 and #SNPs > 1) containing 16 candidate SNPs ( Figure 1F). Of these, only the H2 block, which contained two SNPs (TBX3-EN1, ECA8:20,644,525 and TBX3-EN2, ECA8:20,644,555), overlapped with a 243-bp region showing extremely elevated conservation GERP scores in 90 mammals (score = 504.5, p = 1.51eÀ186). These two SNPs also overlapped an enhancer region previously identified in large-scale mouse embryonic studies 18 ( Figure 1F). Together, this suggested that the TBX3-EN1 and TBX3-EN2 SNPs may play a strong role in the regulation of TBX3 expression ( Figure 1F; Data S2P). We next measured the allelic frequency at the TBX3-EN1 and TBX3-EN2 SNPs in a comprehensive panel of 189 Chinese ponies and 574 Chinese horses using kompetitive allele specific PCR (KASP) assays (Data S3A-S3C; STAR Methods). Chisquare tests confirmed strong statistical association with the horse wither height, which was stronger for TBX3-EN2 (p = 2.34eÀ39) than TBX3-EN1 (p = 8.18eÀ29) ( Figure 2). Our association model suggested that the variation at TBX3-EN2 (TBX3-EN1) could explain as much as $20.3% ($15.1%) of the height variation ($10 cm) in this panel ( Figure 2A).
Interestingly, the A allele at TBX3-EN2 was most common in seven pony breeds (frequency = 65.6%, N = 189) and the Przewalski's horses (frequency = 100%, N = 13). The G allele, however, reached a higher frequency in 18 horse breeds (frequency = 81.2%, N = 561). The frequency of the G allele was most limited (21.1%) in Debao ponies (wither height < 100 cm) but the highest in YiLi horses (85.7%, wither height > 150 cm) and Northern horse  . This confirmed a model in which the A variant could represent the ancestral allele at TBX3-EN2 and the derived G allele increased in frequency among Chinese horses in response to selection for higher wither height. The distribution of the genetic variation at TBX3-EN1 in our panel suggested a possibly similar scenario for the T allele.
We further reconstructed the frequency trajectories for the TBX3-EN1 and TBX3-EN2 derived alleles during the last $6,000 years, leveraging an extensive ancient genome time series previously released for horses. 19 We found that most ancient samples carried the A allele at TBX3-EN2 until $2,300 YBP, a time after which the G allele increased in frequency steadily until the present time ( Figures 3A and 3B; STAR Methods). Importantly, post-mortem DNA decay can lead to the accumulation of G/A substitutions in ancient DNA data. 20 The allelic trajectory reconstructed here was, however, not correlated to postmortem DNA damage levels, as measured through mapDam-age2 deamination parameters 21 (Figures 3B and S3). Therefore, our analyses supported ongoing selection for the G allele at TBX3-EN2 in relation with increased wither height from the late first millennium BCE. The T allele present at the TBX3-EN1 also strongly associated with horse wither height and was randomly distributed at the beginning but followed a similar trajectory with the increase of the TBX3-EN2 G frequency, in line with the selection acting on linked genetic variants ( Figure 3B). Moreover, the genotyping profiles of 245 ancient horses 19 within the 10-kb region flanking the TBX3 enhancer containing TBX3-EN1 and TBX3-EN2 were closer to those found in Chinese ponies (N = 52) than in Chinese horses (N = 100) ( Figure 3C). Modern Chinese horses were characterized by limited sequence variation, in line with the positive selection acting in the region.

Functional validation of the TBX3 causative variant
The analyses based on ancient and modern genomic data supported a model in which the selection for the T mutation at TBX3-EN1 and G mutation at TBX3-EN2 may have driven increased sizes in Chinese horses over the last 2,300 years. Reporter assays in HEK293T cells, however, indicated no significant difference in expression levels between the vector carrying both derived alleles at TBX3-EN1 and TBX3-EN2 (T and G), and that carrying the derived allele at TBX3-EN2 only (C and G; Figure 4A). This contrasted with the significant drop detected for vectors carrying the ancestral A allele at TBX3-EN2, regardless of the allele associated at TBX3-EN1 (T and A or C and A). Together, this indicated that the G mutation at TBX3-EN2 alone was sufficient to increase TBX3 expression levels and, thus, likely represented the causative allele. We then speculated that the G allele at TBX3-EN2 could affect the TBX3 enhancer region and TBX3 expression and may have driven larger wither height in Chinese horses during the last $2,300 YBP. Interestingly, the horses examined for placental transcriptome expression supported such a model, as all four YiLi horses were homozygous G at TBX3-EN2 while all three NiQiang ponies were homozygous A at TBX3-EN2 ( Figure 1E).
In order to further validate the functional consequences of TBX3 mutations, we constructed knockout (KO) mice for the TBX3 enhancer using the CRISPR/Cas9 technology (STAR Methods), in which a core region showing high sequence similarity among vertebrates and including the TBX3-EN2 SNP was deleted (Data S4A; Figure 4B). Such TBX3 enhancer À/À mice showed a 6-fold reduction in TBX3 forelimb transcription relative to the WT (p = 0.004). Similar reduction in TBX3 transcription levels was also detected in embryos and 2-week leg bones (p < 0.05) (Figures 4C and S4A). No differences in TBX3 transcription were detected in muscle, skin, or heart tissues. Our results thus supported that the 35-bp region acts as a limb-specific enhancer for TBX3 expression, in line with dual-luciferase reporter experiments in the HEK293T cell lines ( Figure 4A).
Finally, we used micro-CT and X-ray imaging techniques on 11-and 2-week mice to measure the length of all bone structures involved in limb formation (i.e., the humerus, radius, metacarpal, and fingers for forelimbs and the femur, tibia, metatarsal, and toes for hindlimbs) ( Figure S4B). The 11-week mice included the littermate of six males (WT:KO = 3:3) and six females (WT:KO = 3:3), while 2-week mice include the littermate of ten males (WT:KO = 3:7) and eight females (WT:KO = 5:3). KO and WT mice showed statistically significant differences in the length of all measured limb bones, with size reduction in KO mice ranging from 1.06% to 21.42% (p < 0.05). This finding was consistent across sexes and the two developmental stages investigated (Data S4B-S4D; Figures 4D, 4E, and S4C). Interestingly, the size reduction was more pronounced in forelimbs than hindlimbs, especially in 2-week female mice, reaching almost 20% for the third and fourth metacarpal digits (Figures 4E and  S4C). This was in agreement with previous studies reporting forelimb malformation in TBX3 mice mutants. 22,23 In this study, we used whole-genome sequence data to identify candidate regions underlying the striking differences in wither height found among Chinese horse and pony breeds. One such region encompassed the TBX3 gene and showed strong signatures of positive selection, indicating possible functional relevance. We further restricted and functionally validated genotype-phenotype association for one single A/G nucleotide variant located within the TBX3 enhancer (ECA8:20,644,555), with the G allele predominantly found in horses with increased TBX3 transcription in cell reporter assays, the placenta, and during limb development. Mice models deleted for the 35-bp region carrying the mutation within the TBX3 enhancer also showed reduced limb sizes at 2-and 11-week developmental stages. Altogether, our work demonstrated that TBX3 expression played a crucial role in driving body size variation in Chinese horses and ponies. Interestingly, the A allele appeared ubiquitous among (E) Metacarpal length of TBX3 enhancer in KO À/À and WT mice at the age of 2 weeks; the x axis 1-5 means the five digits (refer to D) and sex separately.

ll
OPEN ACCESS ancient horses predating the Iron Age in the most extensive ancient genome time series available for horses. 19 The G allele was first detected among Berel' Pazyryk horses some $2,300 YBP in the Altay region 24 and steadily increased to presentday frequencies thereafter. This supported a scenario in which past breeders started to select horses of larger sizes concurrently with the beginning of the Chinese Empire during the Qin dynasty. 25,26 Importantly, all Przewalski's horses investigated carried the A ancestral allele, which was also predominant among Chinese ponies. Przewalski's horses and Chinese pony breeds were, thus, not subject to selection for size during their recent evolutionary history. As the G allele, most frequent in horses, is also dominant in western ponies and miniature horses, including Falabellas, our work unveiled different mechanisms underlying size shifts in the history of horse breeding. It also offered a new genetic marker to the horse breeding industry for size improvement and/or pony conservation.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors have a patent related to this work. Report RESOURCE AVAILABILITY

Lead contact
Further information and requests for resources, material and reagents should be addressed and will be fulfilled by the lead contact, Lin Jiang (jianglin@caas.cn).

Materials availability
All plasmids and mouse lines generated in this study are available upon request.

Data and code availability
The NCBI accession number for the genomic data reported in this paper is PRJNA780752. All other previously published genomic data used in this study are available in the key resources table. This paper does not report original code. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

EXPERIMENTAL MODEL AND SUBJECT DETAILS
Candidate SNP genotyping in a larger population We used KASP (Kompetitive Allele-Specific PCR) 43 as a cost-effective approach to genotype a panel of 763 horses and ponies originating from 26 Chinese populations at the for two TBX3 enhancer positions (TBX3-EN1 and TBX3-EN2; Data S3A and S3B). The allele-specific primers were designed by carrying the standard FAM and HEX tails and with the targeted SNP at the 3 0 end (Data S3C). PCR cycling was performed using the following protocol: 1) hot start, 95 C, 15 min; 2) ten touchdown cycles, 95 C, 20 s; 3) touchdown at 65 C initially and decreasing by À1 C per cycle for 25 s; and 4) annealing, 95 C, 10 s, 57 C, 60 s, 30 cycles. The genotyping of KASP markers was conducted at the China Golden Marker (Beijing) Biotech.

Dual-luciferase expression assays
Dual-luciferase assays were carried out in HEK293T (Human embryonic kidney 293 cells). HEK293T cell line was purchased from Beijing university medical college hospital. HEK293T cells were cultivated in Dulbecco's Modified Eagle Medium DMEM basal medium (Glibco, USA), supplemented with 10% heat-inactivated fetal bovine serum and penicillin (0.2 U/mL)/streptomycin (0.2 mg/mL)/ L-glutamine (0.2 mg/mL) (GIBCO). The 63-bp TBX3 enhancer region carrying both TBX3-EN1 and TBX3-EN2 mutations was isolated by PCR with the primers listed in Data S3D, following purification and digested with KpnI and NheI restriction enzymes before being inserted into PGL4.23-Basic vectors (EasyGeno, QIAGEN). A total of four vectors displaying different combinations of mutations of TBX3-EN1 and TBX3-EN2 were constructed and tested with reporter assays for regulatory expression effects. These included: the horse derived mutation (T and G), the TBX3-EN2 single mutation (C and G), TBX3-EN1 single mutation (T and A), and the double TBX3-EN1 and TBX3-EN2 ancestral allele (C and A), where the first letter provides the TBX3-EN1 allele, and the second letter the TBX3-EN2 allele. The TBX3 enhancer region was constructed into dual luciferase reporter vector pGL4.74 and transfected into HEK293T cell by Lipofectamine3000 (ThermoFisher Scientific, USA), following the procedure provided by the manufacturer.