Determination of genetic effects and functional SNPs of bovine HTR1B gene on milk fatty acid traits

Background Our previous genome-wide association study (GWAS) on milk fatty acid traits in Chinese Holstein cows revealed, the SNP, BTB-01556197, was significantly associated with C10:0 at genome-wide level (P = 0.0239). It was located in the down-stream of 5-hydroxytryptamine receptor 1B (HTR1B) gene that has been shown to play an important role in the regulation of fatty acid oxidation. Hence, we considered it as a promising candidate gene for milk fatty acids in dairy cattle. In this study, we aimed to investigate whether the HTR1B gene had significant genetic effects on milk fatty acid traits. Results We re-sequenced the entire coding region and 3000 bp of 5′ and 3′ flanking regions of HTR1B gene. A total of 13 SNPs was identified, containing one in 5′ flanking region, two in 5′ untranslated region (UTR), two in exon 1, five in 3′ UTR, and three in 3′ flanking region. By performing genotype-phenotype association analysis with SAS9.2 software, we observed that 13 SNPs were significantly associated with medium-chain saturated fatty acids such as C6:0, C8:0 and C10:0 (P < 0.0001 ~ 0.042). With Haploview 4.1 software, linkage disequilibrium (LD) analysis was performed. Two haplotype blocks formed by two and ten SNPs were observed. Haplotype-based association analysis indicated that both haplotype blocks were strongly associated with C6:0, C8:0 and C10:0 as well (P < 0.0001 ~ 0.0071). With regards to the missense mutation in exon 1 (g.17303383G > T) that reduced amino acid change from alanine to serine, we predicted that it altered the secondary structure of HTR1B protein with SOPMA. In addition, we predicted that three SNPs in promoter region, g.17307103A > T, g.17305206 T > G and g.17303761C > T, altered the binding sites of transcription factors (TFs) HMX2, PAX2, FOXP1ES, MIZ1, CUX2, DREAM, and PPAR-RXR by Genomatix. Of them, luciferase assay experiment further confirmed that the allele T of g.17307103A > T significantly increased the transcriptional activity of HTR1B gene than allele A (P = 0.0007). Conclusions In conclusion, our findings provided first evidence that the HTR1B gene had significant genetic effects on milk fatty acids in dairy cattle. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-021-07893-8.


Background
Milk fat, a vital nutritional ingredient of milk, is considered as one of the economic traits of milk production in dairy cattle [1]. Triglyceride synthesized by fatty acids and α-glycerophosphate in mammary epithelial cells is the main component of milk fat [2]. Fatty acids contain saturated fatty acids (SFAs) and unsaturated fatty acids (UFAs). In SFAs, C12:0, C14:0, C16:0 increases lowdensity lipoprotein cholesterol and risk of cardiovascular diseases [3], while C15:0, C17:0 are inversely associated with cardiometabolic risk [4]. UFAs are beneficial for reducing the risk of heart and other diseases [5,6]. Many previous studies have shown that the phenotypic variation of milk fatty acid compositions were genetically controlled and the heritability estimates were around 0.14~0. 33 for SFAs and 0.08~0. 29 for UFAs in Holstein cattle [7][8][9][10][11].
One of genome-wide significant SNPs, BTB-01556197 associated with C10:0 (P = 0.0239) identified in previous GWAS [12], was located in the down-stream of 5hydroxytryptamine receptor 1B (HTR1B) gene. The bovine HTR1B gene has merely one exon spanning 4305 bp and was involved in the c-AMP signaling pathway that was related to PI3K-Akt pathway, a well-known pathway for fat synthesis and metabolism [13]. Hence, we considered the HTR1B gene as a promising candidate for milk fatty acid traits in dairy cattle. In the present study, we aimed to further preform association analysis in a different Chinese Holstein population to confirm the genetic effects of HTR1B on milk fatty acids and to identify potential functional genetic variations.
Changes of the HTR1B protein secondary structure and function caused by the missense mutation g.17303383G > T Using SOPMA software, we predicted that the missense mutation in exon 1, g.17303383G > T, changed the HTR1B protein secondary structure, including α-helix (46.90 to 41.44%), extended strand (13.40 to 16.38%), βturn (2.48 to 1.74%), and random coil (37.22 to 40.45%) with the allele from G to T. While, the HTR1B protein function was not altered by the missense mutation with the scores 0.25 (SIFT) and − 1.42 (PROVEAN).
Changes of transcriptional activity caused by g.17307103A > T, g.17305206 T > G and g.17303761C > T By searching the TFBSs of the two SNPs in 5′ UTR and one SNP in 5′ flanking region with Genomatix, we discovered that the allele T of g.17307103A > T created two TFBSs for HMX2 (Hmx2/Nkx5-2 homeodomain transcription factor) and PAX2 (Zebrafish PAX2 paired domain protein), and the allele G of g.17305206 T > G created a TFBS for FOXP1ES (Alternative splicing variant of FOXP1, activated in ESCs). For g.17303761C > T, the allele C created two TFBSs for MIZ1 (Myc-interacting Zn finger protein 1, zinc finger and BTB domain containing 17) and CUX2 (Cut-like homeobox 2, dimeric binding site), and the allele T created two TFBSs for DREAM (Downstream regulatory element-antagonist modulator, Ca2 + −binding protein of the neuronal calcium sensors family that binds DRE sites as a tetramer) and PPAR-RXR (PPAR/RXR heterodimers, DR1 sites). The detailed results were shown in Table 2.
Further, we utilized the luciferase assay (Fig. 2) to confirm the above prediction results. The luciferase activities of the construct T of g.17307103A > T was observed significantly higher than those of the blank control (P < 0.0001), the empty vector PGL4.14 (P = 0.0022), and the construct A (P = 0.0007), indicating that the allele T of g.17307103A > T increased the transcriptional activity of HTR1B compared with allele A. However, the luciferase activities of four constructs (T and G of g.17305206 T > G, and C and T of g.17303761C > T) were not strongly changed than those of the blank and empty vector, implying that g.17305206 T > G and g.17303761C > T did not significantly alter the transcriptional activity of HTR1B gene.

Discussion
According to our previous GWAS results and biological functions, the HTR1B gene has been identified as one of novel promising candidates for milk fatty acids in dairy cattle [12]. In the present study, we firstly confirmed that the HTR1B gene showed significant genetic effects on medium-chain saturated fatty acids in dairy cattle, providing a basis for further verification. Previous studies reported that protein secondary structure could be used to build safe starting cores to generate the complete protein fold [14], and to set structural constraints for protein threading [15,16]. Studies reported that the missense mutations caused by the sequence variations were related to the protein function to account for the phenotype variations [17][18][19][20]. Here, we identified a missense mutation (g.17303383G > T), and it changed the protein secondary structure by prediction with the SOPMA. While we used SIFT and PROVEAN softwares to detect that the missense mutation did not alter the HTR1B protein function. Hence, the significant associations of g.17303383G > T with milk fatty acids might be due to the LD with the true causal mutations.
Regulatory region, including promoter, enhancer, silencer, and insulator etc. [21], are important for the gene regulation and expression [22]. Transcription factors (TFs) are the sequence-specific DNA-binding proteins that regulate gene expression in all organisms [23,24], and approximately 10% are implicated in diverse diseases in human [25]. In eukaryotes, multiple TFs cooperatively bind regulatory DNA to temporally and spatially control gene expression [26]. Generally speaking, 5' UTR plays regulatory roles in the gene expression through binding transcription factors or unknown regulatory mechanisms. Hence, we wish to investigate if the SNPs in the 5' UTRs changed the expression of HTR1B gene. It has been well-known that   [27,28]. Previous studies found that both 5' UTR and 5' flanking region can combine with TF, and some SNPs within these regions may change the TF-binding thereby leading to changes in gene transcription activity and eventually affecting related traits [29][30][31]. In this study, we predicted that three SNPs (g.17307103A > T, g.17305206 T > G and g.17303761C > T) changed the TFBSs. While only g.17307103A > T actually changed the transcriptional activity of HTR1B gene. The allele T of g.17307103A > T was predicted to create two TFBSs for HMX2 and PAX2. The hmx homeoboxcontaining TF gene family, containing contains hmx2 [32], is highly conserved across species [33][34][35][36]. HMX2 is involved in a feedback loop of EGF signaling and located in upstream of the PAX5 of the utricular maculae to affect the inner ear development in zebrafish [37,38]. The pax gene family encodes DNA binding TFs that control vital steps in embryonic development and differentiation of specific cell lineages in human [39]. PAX2 as a TF can promote the expression of ADAM10 to negatively regulate the epithelia to mesenchyme transition in human [40]. According to the significant associations and transcriptional activity caused by g.17307103A > T, we suggested that g.17307103A > T might be a potential causal mutation regulating the HTR1B gene expression by altering the binding sits for the TFs HMX2 and PAX2 to impact the milk fatty acid traits in dairy cattle. Further investigation is needed to validate the regulatory of specific transcription factors. Regarding the expression of the HTR1B gene in multiple tissues, based on the RNA-seq database, Cattle Gene Atlas (http://cattlegeneatlas.roslin.ed.ac.uk/) [41], we observed that the HTR1B gene was expressed in 82 tissues/ cell types, including mammary gland, while HTR1B gene was moderately expressed in mammary gland.

Conclusions
In this study, we confirmed the significant genetic effects of the HTR1B gene on milk fatty acids using post-GWAS strategy, and identified a potential functional mutation in 5′-flanking region, g.17307103A > T, that altered the transcriptional activity of HTR1B. Our findings provided valuable molecular information for genetic improvement programs in dairy cattle.

Gas chromatograph
Total milk fat were extracted from approximately 2 ml of each milk sample. Fatty acids have high boiling points, so they are unstable and easy to crack at high temperatures. The high temperature in gas chromatographic analysis will cause the loss of fatty acids, so pre-treatment is required. Therefore, when analyzing fatty acids and fats, especially fatty acid components, to reduce the boiling point and improve stability, we react the fatty acids or fats with methanol to prepare fatty acid methyl esters and then perform gas chromatography analysis.
The measurement conditions of the gas chromatograph are as follows: the injector temperature is 260°C; carrier gas (helium) flow rate is 45 mL/min; split ratio is 100:1; chromatographic column conditions: keep at 100°C for 10 min, heat up at 6°C /min to 160°C and hold for 10 min, heat up at 5°C /min to 200°C and hold for 20 min, heat up at 4°C /min to 240°C and hold for 12 min; the detector temperature is 260°C.

SNP identification and genotyping
We extracted genomic DNA from the blood samples of 1065 cows and the semen samples of 44 sires using TIA-Namp Blood DNA Kit (Tiangen, Beijing, China) and salt-out procedure, respectively. DNA quantity and quality were measured by NanoDrop™ ND-2000 Spectrophotometer (Thermo Scientific, Hudson, DE, USA) and 2.0% agarose gel electrophoresis.
A total of 15 pairs of primers (Additional file 4: Table S4) were designed for PCR amplification using the Primer 3.0 (http://primer3.wi.mit.edu/) based on the sequences of all the exons, and 3000 bp of 5′ and 3′ flanking regions of the bovine HTR1B gene (Gene ID: 317707), and were synthesized in the Beijing Genomics Institute (Beijing, China). By using the DNA samples of the abovementioned 44 sires with the same concentration of 50 ng/μl, two DNA pools were constructed and 22 sires were included in each pool. Then, PCR amplification was performed with abovementioned 15 pairs of primers and PCR procedure was as follows: initial denaturation at 94°C for 5 min; annealing at 94°C for 30 s, 60°C for 30 s and 72°C for 30 s, for 35 cycles and final extension at 72°C for 7 min. By sequencing PCR products, we identified potential polymorphic sites.
Then, 1065 cows were individually genotyped by using matrix-assisted laser desorption/ionization time of flight mass spectrometry (MALDI-TOF MS, Sequenom Mas-sARRAY, Bioyong Technologies Inc. HK). As for each identified SNP, PCR amplification was first performed with sequence-specific extension primers, then 1 base was extended targeting two alleles of the identified SNP. According to different mass-to-charge ratios of two alleles, different mass spectrum peaks could be observed to detect the genotype of each SNP.

Statistical analysis
First, we used the Haploview 4.1 software (Broad Institute of MIT and Harvard, Cambridge, MA, USA) to identify the LD extent among the identified SNPs of the HTR1B gene.
Subsequently, we performed single SNP-based and haplotype-based association analysis. We traced the pedigrees of the 1065 Chinese Holstein cows back to three generations, as a result, a total of 3335 individuals were included for association analysis, which kinship matrix (A-matrix) were constructed with SAS9.2 (SAS Institute, Cary, NC, USA). Then, associations between the identified SNPs and haplotype blocks with 23 milk fatty acid traits were performed by SAS9.2 on the basis of the following mixed animal model: Here, Y ijklm is the phenotypic value of each milk fatty acid trait; μ is the overall mean; G i is the fixed effect corresponding to the genotype or haplotype combination; h j is the fixed effect of farm (j = 1~23); l k is the fixed effect of stage of lactation (k = 1~4); a l is the random polygenic effect; M m is the fixed effect of age at calving (m = 1~293); b is the regression coefficient of covariate M; and e ijklm is the random residual. Bonferroni correction was performed according to the number of multiple tests, in which the adjusted significance levels of P < 0.05 for the single SNP and haplotypebased analysis were 0.0002 and 0.0011, respectively.
Further, we calculated the additive effect (a), dominant effect (d), and substitution effect (α) of SNP on the milk fatty acid traits according to the formulas [43]: ; and α ¼ ɑ þ dðq−pÞ , in which, p and q were the frequencies of A and B, respectively; and AA, AB and BB were the least square means of fatty acids corresponding to the genotypes.
Prediction of the secondary structure and function changes of the HTR1B protein We used the NPSA SOPMA SERVER program (https:// npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_ sopma.html) to predict whether the identified missense mutation in coding region changed HTR1B protein secondary structure, and set the parameters with similarity threshold (8), and number of states (4-Helix, Sheet, Turn, Coil). Also, we used the SIFT (http://sift.bii.a-star. edu.sg/) and PROVEAN (http://sift.jcvi.org/index.php) to investigate whether the missense mutation altered the protein function. The score thresholds of the SIFT and PROVEAN were 0.05 [44] and − 2.5 [45], respectively. When the score is below the threshold, the protein function is considered changed.
Prediction of the changes of transcription factor binding sites (TFBSs) We predicted whether the SNPs in 5′ UTR and flanking region of the HTR1B gene impacted on TFBSs by using the Genomatix suite v3.9 (http://www.genomatix.de/cgibin/welcome/welcome.pl?s=d1b5c9a9015b02bb3b1a806 f9c03293f).

Construction of recombinant plasmid, cell culture and luciferase assay
We constructed six luciferase reporter gene fragments with Kpn1 and Nhel restriction sites at the 5′ to 3′ termini (Figs. 3 and 4), which contained alleles A and T of g.17307103A > T, T and G of g.17305206 T > G, and C and T of g.17303761C > T. The six fragments were synthesized in Genewiz (Suzhou, China), and were cloned into the pGL4.14 Luciferase Assay Vector (Promega, Madison, USA). Subsequently, the plasmids were purified by Endofree Plasmid DNA Mini Kit II (OMEGA, omega bio-tek, Norcross, Georgia, USA), and then were sequenced to confirm the integrity of each construct's insertions. Human Embryonic Kidney (HEK)-293 T cells were cultured in Dulbecco's modified Eagle's medium (DMEM) including 10% heat-inactivated fetal bovine serum (FBS; Gibco, Life Technologies) at 5% CO 2 and 37°C. We seeded approximately 2 × 10 5 cells per well into 24-well plates, and transfected the cells using Lipofectamine 2000 (Invitrogen, CA, USA). For each well, we transfected 500 ng of the constructed plasmid DNA along with 10 ng of pRL-TK renilla luciferase reporter vector (Promega), and conducted three replicates for each construct. The cells were cultured for about 36-48 h after transfection and then were measured the activities of firefly and renilla luciferases using a Dual-Luciferase Reporter Assay System (Promega) with a Modulus microplate multimode reader (Turner Biosystems, CA, USA). Finally, average statistics of three replicates were calculated as the normalized luciferase data (firefly/renilla).