All of the animal procedures were performed in strict accordance with the guidelines proposed by the China Council on Animal Care and the Ministry of Agriculture of the People’s Republic of China.

After slaughter, gender, age, CW, CL (from the basis of tail to the basis of the neck), TVN, LVN, and thoracolumbar vertebral number (TLVN) were collected in the abattoir in Yining City of Xinjiang, China, from November to December in 2014.The age which was estimated by inspection of the teeth at slaughter ranged from 1 to 6 years.

DNA was extracted from ear tissue using the standard phenol-chloroform method [14]. The quality and quantity of the DNA were measured with a NanoDrop 2000 (Thermo Fisher Scientific Inc., Wilmington, DE, USA). DNA was at the concentrations >50 ng/μL, total volume >50 μL, and the ratio of light absorption (A260/280) between 1.8 and 2.0.

According to sheep VRTN gene sequence in GenBank (Accession number: NC_019464.1), two sets of polymerase chain reaction (PCR) primers (Table 1) were designed using the Oligo 6.0 software and were synthesized by Sangon Biotech (Shanghai) Co., Ltd. Shanghai, China.

Two PCR fragments which were produced by use of primers above harbored the complete VRTN gene and PCR amplifications were performed with EasyTaq PCR SuperMix (−dye) (Transgen Biotech Co. Ltd. Beijing, China) in accordance to manufacturer’s instructions. The conditions for the PCR were as follows: 94°C for 5 min; 30 cycles of 94°C for 45 s, 59°C for 45 s, and 72°C for 90 s; followed by a further 5 min extension at 72°C. All PCR products were bi-directionally sequenced with original PCR primers on the ABI-PRISM3730 sequencer at Sangon Biotech (Shanghai) Co., Ltd. The sequence alignment and detection of mutations were conducted using the MEGA 6.0 software and Chromas 2.5 software.

All statistical analyses were performed by SPSS 19.0 for Windows (SPSS Inc., Chicago, IL, USA). Univariate analysis of variance (general linear model, GLM) was applied to estimate the effect of vertebral number variation on CW and length. All data were presented as estimated marginal means (EMMEANS)±standard error [15]. The model was:

Where Y _ijk was the phenotypic value of each trait (CW or CL), with μ as general mean; A_i = fixed effect of the ith age at traits recorded (i = 1 to 6 ); G_j = fixed effect of jth gender (j = 1 or 2, representing male or female, respectively); V_k = fixed effect of vertebral number, was replaced by T_k (k = 1, 2, or 3, representing 12, 13, or 14 thoracic vertebrae, respectively), L_K (k = 1 or 2, representing 5 or 6 lumbar vertebrae, respectively ), P_k (k = 1, 2, or 3, representing 18, 19, or 20 thoracolumbar vertebrae, respectively), or F_k (k = 1, 2, 3, or 4, representing T12L6 (12 thoracic and 6 lumbar vertebrae), T13L5 (13 thoracic and 5 lumbar vertebrae), T13L6 (13 thoracic and 6 lumbar vertebrae) or T14L6 (14 thoracic and 6 lumbar vertebrae, respectively), when the effect of each kind of vertebra trait (TVN, LVN, TLVN or vertebral formula) was assessed (one at a time). E_ijk was the random residual effect.

The association between the mutations in VRTN gene and vertebral number was performed by one-way analysis of variance (ANOVA) analyses. The additive effect of allele was calculated with: (AA–BB)/2 (AA: mean phenotypic value of homozygous individuals, BB: mean phenotypic value of the other homozygous individuals).

In present investigation, a total of 624 Chinese Kazakh sheep were sampled for phenotype data including gender, age, CW, CL, and the number of thoracic and lumbar vertebrae. The ratio of TLVN variation is presented in Figure 1. The variation range observed for TVN was 12 to 14, LVN was 5 to 6, and TLVN was 18 to 20. The range of TLVN variation (18 to 20) in China Kazakh sheep was smaller than that in Texel and Scottish Blackface sheep (17 to 21) based on previous reports [3]. The further analysis showed that the range of TVN variation in Kazakh sheep was the same as that in Texel sheep, but the range of LVN (5 to 6) is much smaller than that in Texel sheep (4 to 8). Remarkably, individuals with 7 lumbar vertebrae accounted for relatively larger proportion in Texel and Scottish Blackface sheep (35.4% and 60.8%, respectively) [3,16], but those have not been found in Kazakh sheep.

Frequency distributions revealed that the percentage of in dividuals with 19 thoracolumbar vertebrae was up to 79.2%, obviously higher than that of the total samples with 18 or 20 thoracolumbar vertebrae (15.1% and 5.8%, respectively) (Figure 1a). The vertebral formulae consisted by T12L6, T13L5, T13L6, and T14L6, whereby T13L6 was the most common thoracolumbar vertebral composition (Figure 1b). In mammals, the number of thoracolumbar vertebrae was exclusively fixed at 19 (T13L6), so T13L6 was thought to be the primitive vertebrae form [1]. However, the changes in TLVN (from 21 to 23) were common in western modern pig breeds, such as Duroc, Landrace, and Large White, because these breeds have long been selectively bred for enlargement of body size and have presumably increased the number of vertebrae [6,17]. The percentages of animals possessing the extra vertebral counts were very low (Figure 1) in our investigation, which was reasonable, because Kazakh sheep are an old indigenous breed that was seldom subjected to artificial selection (cross-breeding). On the other hand, an investigation reported decades ago indicated that animals with more thoracolumbar vertebrae suffered from movement problems [18,19]. We speculated that the locomotion problem caused by extra vertebrae maybe was too slight to be observed when animals were grazing on flat pasture or in the breeding house. But Kazakh sheep originally lived in diverse terrain and animals with movement problem caused by more thoracolumbar vertebrae were eliminated by nature selection. This was the main second reason that the frequency of Kazakh sheep with more thoracolumbar was obviously lower than that of Texel and Scottish Blackface sheep breeds.

Statistical analysis by GLM indicated that the TLVN affected CW and CL apparently (Table 2). The effect of TLVN on CL was more significant than on CW, and the effect on CL was highly significant (p<0.01).

The TVN and LVN independently were significantly corre lated with CL (Supplementary Table S1, S2). Moreover, the effect of TVN on CL was larger than that of LVN. An extra thoracic vertebra increased the CL by around 17 mm, and an extra lumbar vertebrae by 13 mm respectively. Notablly, the only variation in vertebral formulae with the same number of thoracolumbar vertebrae had little affect on CW and CL (Supplementary Table S3), however, the extra thoracic vertebrae not only apparently increased CL and longissimus muscle length which potentially positively associated with longissimus muscle weight, but also increased a pair of spareribs, hence it was a very important economical trait.

Furthermore, both extra thoracic and extra lumber vertebrae tended to increase CW, while it didn’t reach to significant level statistically (Supplementary Table S1–S2), which was in line with the results in pigs [20]. We infered that CW was a complex quantity trait and was affected by multiple factors, especially by the fat weight [21]. Kazakh sheep has the character of typical fat-rumped tail. The rump weight is about 3.5 kg in average at maturity and is account for 16% of CW, which may have more influence on CW rather than on CL. So the effect of TLVN on CL was more significant than on carcass weight statistically.

In addition, previous studies showed that vertebra numbers were associated positively with leaner carcass, loin lengths and growth rate, and associated negatively with intramuscular fat content and fat deposition [2,11,20,22]. However, the association in sheep has not been addressed comprehensively.

A total of 9 single nucleotide polymorphisms (SNPs) were detected in VRTN CDS, including 2 missense and 7 synonymous SNPs (Supplementary Table S4) by sequencing of 624 lambs. Of the 9 mutations, the SNP1259 (rs426367238) had significant effect on the number of thoracic vertebrae by ANOVA analyses (Table 3). Sheep with the C/C genotype had significantly more TVN (13.67±0.21) than those with the C/G genotype (13.22±0.06, p = 0.022), and highly significantly more TVN than those with the G/G genotype (12.98±0.04, p = 0.004).

Subsequent analyses revealed that the SNP1259 (rs426367238) had an additive effect on the thoracic vertebrae number. The allele C with lower frequency (0.23, calculated from Table 2) in our investigation was most presumably a causal genotype for increasing of vertebrae number. The effect of SNP1259 was an increase of approximate 0.4 thoracic vertebrae per allele, which was less than the additive effect of VRTN mutation in pigs (0.5 to 0.6 additional thoracic vertebrae per allele) [11]. Therefore, we speculated that there could be a number of other loci that affect TVN. On the other hand, none of the 9 SNPs had significant effect on the number of lumbar vertebrae. It inferred that the VRTN gene had strong association with thoracic vertebrae number variation, rather than with the number of lumbar vertebrae, which was consistent with previous result in pig [10,12].

Our findings suggested that VRTN was the responsible gene for the variation of thoracic vertebrae number in sheep. VRTN was therefore likely to be an essential regulator for somitogenesis in the development of embryo in a wide range of vertebrate [9,10]. Nevertheless, the detailed underlying mechanism of VRTN in the development of vertebral formula remains largely unknown at present and needs further investigation.