SERPINA1 gene identified in RNA-Seq showed strong association with milk protein concentration in Chinese Holstein cows

Animal ethics

All protocols for collection of the blood and frozen semen samples of experimental individuals were approved by the Institutional Animal Care and Use Committee (IACUC) at China Agricultural University (Permit Number: DK996). The procedures for sampling, laboratory analysis, and data processing were following our previously published article (Li et al., 2019).

Animal, phenotypes and traits

All samples and related data were collected from the cows from 17 farms of the Beijing Sanyuan Lvhe Dairy Farm Center (Beijing, China). A total of 1,027 Chinese Holstein cows from 17 sire families were involved in the present experiment with the family size ranging from 25 to 187 daughters with an average of 60 daughters per sire. All cows enrolled in the experiment were regularly tested using standard performance testing (Dairy Herd Improvement, DHI) and five milk production traits were recorded, including 305 d milk yield, 305 d protein yield, 305 d fat yield, average 305 d protein percentage, and average 305 d fat percentage.

The tubule frozen semen samples were collected for all 17 sires and prepared by Beijing Bull Station. Whole blood samples were collected from coccygeal vein of 1,027 Chinese Holstein cows using a 20-gauge BD Vacutainer needle (Beckton Dickinson, Franklin Lakes, NJ) and immediately stored at −20 °C prior to DNA extraction.

Genomic DNA extraction

Genomic DNA was isolated from semen samples using standard phenol-chloroform procedures as described in previous study (Li et al., 2016b), whereas genomic DNA was extracted from blood samples using a TIANamp Blood DNA kit (TIANGEN Biotech, Beijing, China) following the manufacturer’s instructions. The quality and quantity of genomic DNA were measured using the gel electrophoresis and a NanoDrop 2000 spectrophotometer (Thermo Scientific, Hudson, DE, USA).

Identification and genotyping of the SNPs

A total of 11 pairs of primers (Table S1) were designed using Primer3web Program (v.0.4.0) to amplify all exons, partial adjacent introns and 2,500 bp upstream of 5′ flanking sequences based on the genomic sequence of bovine SERPINA1 gene referring to the UMD3.1 assembly (NCBI reference sequence accession no. AC_000178.1). A DNA pool was generated with genomic DNA from the 17 sires’ semen with 50 ng/µL/bull and was amplified with Polymerase Chain Reaction (PCR). Each PCR reaction consisted 50 ng genomic DNA, 0.5 µL of forward primer, 0.5 µL of reverse primer, 2.5 µL 10× PCR buffer, 2.5 mM dNTP, and 1 U of Taq DNA polymerase (Takara Biotechnology Co., Ltd., Dalian, China) with a total reaction volume of 25 µL. Thermal cycling conditions were: 94 °C for 5 min, followed by 35 cycles at 94 °C for 30 s, 56 °C for 30 s, and 72 °C for 30 s, a final extension at 72 °C for 7 min for all primers. PCR products were confirmed by gel electrophoresis on 2% agarose gels and by photography under UV light, and were bi-directionally sequenced by ABI 3730XL DNA analyzer (Applied Biosystems, Foster, CA, USA). And the sequences were aligned to the bovine reference sequences (UMD3.1) using BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to identify potential SNPs.

The matrix-assisted laser desorption/ionization time of flight mass spectrometry assay (MALDI-TOF MS, Sequenom MassARRAY, Bioyong Technologies Inc. HK) was further applied for individually genotyping of the identified SNPs in 1,027 Chinese Holstein cows.

Linkage disequilibrium (LD) analysis and haplotype construction

Based on the criterion of D′, the extent of linkage disequilibrium (LD) was measured between each pair of SNPs that were genotyped within the SERPINA1 gene using the software Haploview (Barrett et al., 2005). Briefly, genotypes were imputed for each individual using the Beagle3.2 software program (Browning & Browning, 2007). After that, haplotype blocks with high LD of SNPs (D′ > 0.90) were defined based on confidence intervals methods (Gabriel et al., 2002). Haplotypes with frequencies >5% were considered distinguishable, whereas the haplotypes with relative frequencies <5% were pooled into a single group (Wang et al., 2013). Detailed procedures were described in our previously published research (Li et al., 2019).

Statistical analysis

The effects of single SNPs or haplotypes in SERPINA1 on five milk production traits were analyzed by the Mixed Procedure of SAS (SAS Institute Inc., Cary, NC) with the following model: ${\displaystyle y_{\text{ijklmn}}=\mu +F_i+Y{S}_j+P_k+b\times \mathrm{M}+G_l+{\alpha }_m+e_{\text{ijklmn}}}$

where, y_ijklmnwas the phenotypic value of each trait of cows (n = 1,027 for each trait);  μ was the overall mean; F_iwas the fixed effect of farm; YS_jwas the fixed effect of year-season; P_kwas the fixed effect of parity; M was the covariate effect of calving month; b was the regression coefficient of M; G_l was the fixed effect corresponding to the genotype of polymorphisms or haplotype; α_mwas the random polygenic effect, distributed as N (0, A ${\sigma }_a^2$), with the additive genetic relationship matrix A and the additive genetic variance ${\sigma }_a^2$; and e_ijklmn was the random residual, distributed as N (0, I ${\sigma }_e^2$), with identity matrix I and residual error variance ${\sigma }_e^2$.

The Bonferroni adjustment was used for single SNP and haplotype analyses according to the number of SNP loci or haplotype blocks. Associations were considered as significant if a raw P value < 0.05/N, where N is the number of SNP loci or haplotype blocks tested in analyses. Hardy-Weinberg equilibrium (HWE) tests were performed on each identified SNP. The SNP allele frequencies were calculated and the expected genotype numbers were estimated using the expected genotype frequencies under HWE. Chi-square analysis was used to compare the number of expected genotypes and observed genotypes, using 0.05 as the significance threshold value. The additive (a), dominance (d), and allele substitution (α) effects were estimated according to the equation proposed by Falconer & Mackay (1996), i.e.,  a = (AA − BB)∕2, d = AB − (AA + BB)∕2 and $\alpha =a+d\left(q-p\right)$, where AA and BB represent the two homozygous genotypes, AB is heterozygous genotype, and p and q are the allele frequencies of corresponding loci.

Phenotypic variance

The effect of each individual SNP on a specific trait was measured as the proportion of phenotypic variance of the trait explained by the SNP. The proportion of variance explained by a SNP was calculated as follows, ${\displaystyle \text{Phenotypic variance ratio}=2p\left(1-p\right){\alpha }^2∕{\sigma }_p^2.}$

Where, p is the allele frequency of the analyzed SNP, α is the average effect of gene substitution calculated using the linear mixed model, and ${\sigma }_p^2$ isthe estimate of the phenotypic variance using the complete DHI data of the Chinese dairy cattle population.

Phenotype data

Descriptive statistics of five milk production traits from 1,027 Chinese Holstein cows were presented in Table 1. All phenotypic values of five milk production traits followed approximately normal distributions and were able to be used for the following association studies.

SNPs identification

Through resequencing the entire coding sequences, partial adjacent introns and 2,500 bp upstream of 5′ flanking sequences, a total of four SNPs were identified for the SERPINA1 gene. One (rs208607693) was located in the 5′ -UTR, one (rs207601878) was intronic and the other two (rs210222822 and rs41257068) were located in exonic regions (Table 2). Both exonic SNPs were synonymous substitutions. All four SNPs were in Hardy-Weinberg equilibrium (P > 0.05), and the locations and allele frequencies of the four SNPs were summarized in Table 3.

Associations between SERPINA1 gene and milk production traits

Associations between the four SNPs of SERPINA1 and five milk production traits are presented in Table 4. All four SNPs (rs208607693, rs207601878, rs210222822 and rs41257068) had significant (P = 0.0458 −  < 0.0001) associations with five milk production traits, with the exception that no significant associations were observed between two SNPs (rs208607693 and rs210222822) and milk fat percentage. Phenotypic variances (>1%) explained by the four SNPs in SERPINA1 gene were observed in 11 significant pairs of SNP-trait, among which rs208607693 and rs210222822 presented phenotypic variances of approximately 7% in milk fat yield and protein yield.

The additive, dominant and substitution effects of four SNPs of SERPINA1 on five milk-production traits are shown in Table 5. There were significant (P < 0.05) additive and substitution effects of rs210222822 on milk protein percentage. Significant additive and substitution effects were also observed (P < 0.05) in rs41257068 on protein and fat percentages and observed (P < 0.05) in rs207601878 on fat yield, fat percentage and protein percentage. Significant dominant effects of rs20860769 on protein percentage, rs41257068 and rs207601878 on milk yield and protein yield were also observed (P < 0.05). In addition, significant additive, dominant and substitution effects of rs208607693 and rs210222822 on milk yield, fat yield and protein yield, rs41257068 on fat yield were observed (P < 0.05).

LD and haplotype analysis

A markedly strong linkage (D′ = 0.99–1.00) was observed between the four SNPs in SERPINA1, as shown in Fig. 1. A 9 Kb haplotype block composed of four SNPs was inferred (Fig. 1) and three main haplotypes were formed. The common haplotypes GTGT, CCCC and CCGT occurred at the frequencies of 37.73%, 30.62% and 29.21% respectively (Table 6). Haplotype association analysis showed that haplotypes were highly associated with all five milk production traits (P = 0.0003 −  < 0.0001, Table 7).