Novel Insight Into The Potential Role of AGPAT Gene Family in Triacylglycerols Synthesis in Buffalo

Xiaoya Ma Guangxi Buffalo research institute Anqin Duan Guangxi buffalo research institute Xingrong Lu Guangxi buffalo research institute Md Mahmodul Hasan Sohel Erciyes University Hamdy Abdel-Shafy Cairo University Ahmed Amin Cairo University Shasha Liang Guangxi buffalo research insititute Tingxian Deng (  dtx282000@163.com ) Guangxi Buffalo Research Institute https://orcid.org/0000-0002-7442-6739


Introduction
Water buffalo (Bubalus bubalis) is important dairy livestock producing more than 15% of the world's total milk production [1]. Buffalo milk is a highly nutritious food containing 89% more fat, 28% more protein, and 100% more iron as well as 1% less lactose, and 77% less cholesterol content compared to those in dairy cattle [2]. Milk fat is composed of 98% triacylglycerols (TAG), which are the most abundant lipids in humans and animals. It is well-reported that the TAG determines the functional and physical properties of several dairy products such as the spreadability of butter [3]. The TAG is formed by a combination between glycerol and three fatty acids. The availability of fatty acids presents in milk epithelial cells along with the enzymes involved in the metabolism is largely affecting the characteristics of TAG, which is partially controlled by several genetic factors [4,5]. The synthesis of TAG in eukaryotes is controlled by two major pathways: the Kennedy pathway (de novo pathway) and monoacylglycerol pathway. Most TAG synthesis is performed via the de novo pathway in most cells. Evidence showed that four rate-limiting steps play a vital role in determining the triglyceride content [6], which are controlled by the glycerophosphate acyltransferase (GPAT), acylglycerophosphate acyltransferase (AGPAT), lipid phosphate phosphohydrolase (LPIN), and diacylglycerol acyltransferase (DGAT) gene families, respectively [7].
AGPAT family acts on the second acylation step in the de novo pathway to catalyze the conversion of lysophosphatidic acid (LPA) to phosphatidic acid (PA). To date, genome-wide identi cation, evolution, gene structure, and expression pro le analyses have identi ed eleven gene members (AGPAT1-11) of the AGPAT gene family in different species [8][9][10]. Evidence revealed that the AGAPT family was ancient and experienced different origins, with many eukaryotic species having multiple genes originated by duplication events [8]. For example, the duplication events resulted in the origination of two isoforms AGPAT3 and AGPAT4 in animal species. Gene structure revealed that the number of introns per gene varied from one to thirteen [11]. In addition, four acyltransferase motifs (motifs I~IV) were also found in the animal species. Among them, motifs I and IV contain a conserved NHxxxxD and proline (P) sequence, respectively, involving in the acyl-CoA binding and catalysis, whereas motifs II and III are involved in LPA binding, with a conserved arginine (R) and EGTR sequence, respectively [7,10,11]. Moreover, numerous studies have characterized the expression pro les of the AGPAT gene family in different animals based on different datasets, such as the microarray expression data of human and mouse [8], RNA-seq data of chicken [9], and real-time quantitative PCR (qRT-PCR) data for cattle and pigs [12,13], implying that different members had spatial and temporal speci c expression pattern in different tissues. For example, AGPAT6 appears to be the most abundant AGPAT isoform in the bovine mammary gland, followed by AGPAT1 [13]. Numerous studies also found that AGPAT6 polymorphisms were highly associated with milk production traits in cattle including milk fat percentage [14], fat compositions [15], and fatty acids [16]. Notably, Na kov et al. [17] also found that haplotypes of the AGPAT1 gene were associated with higher PUFA and linoleic acid concentrations in cattle. These ndings suggested that both AGPAT1 and AGPTA6 play a vital role in the TAG synthesis. However, information on how many members of the AGPAT family involved in the TAG synthesis of mammary epithelial cells is still extremely limited.
Since milk production of dairy cows is signi cantly higher than that of buffaloes and the functional importance of AGPATs, it is of considerable importance to characterize the function of buffalo AGPAT members on the TAG synthesis. In addition, the identi cation of excellent functional genes will be of great signi cance to the molecular breeding practices of this species. In this regard, taking dairy cattle as a reference, the objectives of the present study were rst to characterize the differential expression pattern of AGPAT gene family in both cattle and buffalo using comparative analysis; subsequently, to investigate the potential functional role of AGPAT genes on cell growth and TAG synthesis in buffalo mammary epithelial cells (BuMEC) and bovine mammary epithelial cells (BoMEC) using RNA interference (RNAi) technology.

Materials And Methods
Identi cation and sequence analysis of buffalo AGPAT genes To perform the genome-wide identi cation of AGPAT genes in buffalo, we downloaded public genome datasets of seven representative mammals including the human (GRCh38.p12), cattle (ARS-UCD1.2), river buffalo (UOA_WB_1), goat (ARS1), sheep (Oar_rambouillet_v1.0), and horse (EquCab3.0) from the National Center for Biotechnology Information (NCBI) Genome database, as well as the swamp buffalo (accession number: GWHAAJZ00000000) deposited in the Genome Warehouse in the Beijing Institute of Genomics (BIG) Data Center. The known sequences of some AGPAT proteins were downloaded from the UniPort [18] and used to build the hidden Markov models (HMM) pro le with the HMMER v3.3.1 software [19]. All AGPAT proteins were examined in the studied species by HMMER v3.3.1 software [19] with default parameters. Subsequently, these protein sequences were further used for multiple alignments and phylogeny tree construction implemented in MEGA-X [20] software. Identi ed AGPAT protein sequences were subjected to the ExPASy proteomics server for the molecular weight (MW) and isoelectric points (pI).
The motif pattern and element annotation of AGPATs was analyzed using the MEME platform and Pfam programs, respectively. The conserved domain of buffalo AGPATs was predicted using CDD tools of NCBI. The motif pattern and conserved domains of buffalo AGPATs were visualized by TBtools v1.051 software [21].
Chromosomal mapping and collinearity analysis for AGPAT family Chromosome locations of the AGPAT gene family in buffalo and cattle were obtained from their genome resources. Identi cation of orthologous AGPAT genes between buffalo and cattle was analyzed and visualized by TBtools v1.051 software [21] with one-step MCScanX command. Synonymous (Ks) / nonsynonymous (Ka) substitution (Ka/Ks) rate for orthologous AGPAT pairs was calculated by the TBtools v1.051 software [21]. The divergence time for each orthologous gene pairs was evaluated by T = Ks/2λ ×10 -6 million years ago (Mya) [22], where T is the absolute time of divergence, Ks is the synonymous substitution rate, and λ is the clock-like rates in buffalo of 1.26 × 10 −8 [23].

Comparative transcriptomic analyses for orthologous AGPAT genes
To explore the differential expression of AGPAT gene family between buffalo and cattle, two published RNA-seq data (BioProject: PRJNA419906 and PRJNA453843) from milk samples were utilized for conducting the comparative transcriptome analyses. Overall, the quality control of raw data was performed by the Trim galore ver0.6.6 software [24]. Mapping of the cleaned data from buffalo and cattle was conducted by the HISAT ver.2.2.1 software [25], corresponding to the UOA_WB_1 and ARS-UCD1.2 version as reference genomes, respectively. The count matrix of gene or transcript was constructed by the StringTie ver.2.1.4 software [26], and Transcripts per million (TPM) values for each gene were obtained using the DESeq2 R-package [27]. Finally, the differential expression analysis of orthologous AGPATs between buffalo and cattle milk samples was performed. Clustering and generation of a heat map of TPM values for the selected genes were performed using the TBtools v1.051 software [21].

Cell culture and transfection
Both BuMECs and BoMECs were preserved in the Buffalo Research Institute laboratory. BuMECs were cultured in DMEM/F-12 medium (Thermo Fisher Scienti c, USA) supplemented with 10% fetal bovine serum, 5 μg/mL insulin, 1 μg/mL progesterone, 1ug/mL Hydrocortisone, 5 μg/mL prolactin (Sigma-Aldrich, USA), and 10 ng/mL epidermal growth factor 1 (Thermo Fisher Scienti c, USA) and incubated at 37°C in a humidi ed atmosphere of 5% CO 2 level. At the same time, BoMECs were also cultured under the same cultural conditions. The siRNA fragments for AGPAT1 and AGPTA6 were designed and synthesized by GenePharma (Shanghai, China). These siRNA fragment sequences are listed in the Additional le 1: Table S1. Both BuMECs and BoMECs were either transfected with AGPATs siRNAs (siAGPAT) or negative control (NC) siRNAs using the RNAiMAX following the manufacturer's protocol (Thermo Fisher Scienti c, USA), respectively.

Cell viability assay
Cell viability was detected after the knockdown of AGPAT1 or AGPAT6 using CellTiter-Glo® Reagent (Promega, USA) according to the manufacturer's instructions. Brie y, cells were seeded in a 96-well culture plate. At 60% con uency, cells were transfected. After treatments with three time-points (24 h, 48 h, and 72 h), cells were then equilibrated at room temperature for 30 min. Next, 100 μL compound reagents were added to 100 μL of medium containing cells, mixed for 2 min on an orbital shaker to induce cell lysis, incubated at room temperature for 10 min to stabilize the luminescent signal. Control wells containing medium without cells were also prepared to obtain a value for background luminescence, this was followed by luminescence measurement.
Cell apoptosis assay Cell apoptosis assay was performed using Annexin V-EGFP Apoptosis Detection kit (Abcam, Cambridge, USA). According to the manufacturer's protocol, the cells were collected by centrifugation, resuspended in 500 μL of binding buffer, added to 5 μL of Annexin V-EGFP and 5 μL of propidium iodide, and incubated at room temperature for 5 min in the dark. After that, the cell suspension was placed in a glass slide, and results were detected under the uorescence microscope, where the cells which bound Annexin V-EGFP were stained with green in the plasma membrane, while cells that have lost membrane integrity were shown in red staining through the nucleus and a halo of green staining (EGFP) on the cell surface.
Triglyceride assay Cellular TAG content in BuMECs and BoMECs was evaluated by using a TAG assay kit (Applygen, Beijing, China). Brie y, cells were seeded in 24 well plates and cultured until 60-70% con uency. Subsequently, the cultured cells were transfected with siRNA fragments. After 48 h, the culture medium was removed and the cells were washed three times with PBS. The TAG level was calibrated with protein concentration determined with a BCA protein assay kit (Applygen, Beijing, China) and expressed as total TAG per cellar protein. Each experiment was performed in triplicate and repeated at least three times.
Quantitative Real-Time PCR Total RNA was isolated and puri ed with PureLinkTM RNA Mini Kit (Thermo Fisher Scienti c, USA). After that, the First Strand cDNA Synthesis Kit (Thermo Fisher Scienti c, USA) was used to synthesize the rststrand of cDNA. A LightCycler 480 II sequence detection system instrument (Roche, Vienna, Austria) was used to quantify the transcript abundance of the selected genes. The qRT-PCR reactions were set up in 20 μL containing 10 μL Power-up SYBR Green Master Mix (Thermo Fisher Scienti c, USA), 2 μL rst-strand cDNA template, 0.3 μM of forward and reverse gene-speci c primers, and 7.4 μL deionized H 2 O. The expression analysis was performed using a comparative CT(2 -∆∆Ct ) method [28]. The GAPDH gene was used as an endogenous control. The statistical signi cance between different groups was analyzed using a t-test and signi cance was declared at P<0.05. All the primers used for qRT-PCR are list in the Additional le 2: Table S2.

Statistical analysis
Statistical analyses were conducted using the GraphPad Prism 9.0 software [29]. The analyzed data was presented as the mean and standard error of the mean (SEM). Each test was used in triplicate. Statistical signi cance between the two groups was determined by the Tukey post hoc test. The signi cance level was declared at P<0.05.

Results
Identi cation and sequence analysis of buffalo AGPAT genes A total of 32 and 14 AGPAT isoform protein sequences encoded by 13 AGPAT genes were predicted from the river and swamp buffalo genome using the HMMER and BLAST software, respectively (Additional le 3: Table S3). The open reading frames (ORFs) of buffalo AGPAT protein isoforms ranged from 762 to 2,136 bp in length encoding the protein of 253 to 711 residues, with the predicted MW from 28.93 to 78.25 kDa. The pI values of these isoforms ranged from 6.20 to 11.26.
Phylogenetic analysis revealed that all buffalo AGPAT genes could be divided into four clusters ( Figure  1). Cluster IV was the top one with the larger numbers of AGPATs (n = 5), followed by Cluster I with 4 AGPAT genes, while Cluster II and III were the smaller ones (n = 2). The constructed dendrogram further showed that the buffalo AGPAT gene family was usually the most closely evolutionary relationship with the other ve representative mammals (Cattle, Goat, Sheep, Horse, and Human). Moreover, the motif analysis showed that a total of 10 conserved motifs were detected in the identi ed buffalo AGPAT genes (Additional le 6: Figure S1). Here, 4 motifs (MEME-1, MEME-3, MEME-5, and MEME-7) were annotated as the collagen domain after the Pfam search (Additional le 4: Table S4). Interestingly, we observed that the AGPATs in Cluster I had 3 acyltransferase domains (MEME-1, MEME-3, and MEME-7) and one EF-hand_1 domain. The AGPATs in the remaining Cluster II and III had two acyltransferase domains with the MEME-1 and MEME-3 motifs, but the Cluster II AGPATs had an EF-hand_1 domain. The Cluster IV AGPATs had the MEME-3 and MEME-5 acyltransferase domains. Moreover, conserved protein domain analysis revealed that a total of 8 domains were found in the analyzed AGPAT protein sequences (Additional le 6: Figure  S1). AGPATs in Cluster I had the LPLAT_LPCAT1-like conserved domain. The AGPATs in Cluster II and III had the acyltransferase C-terminus and LPLAT_AGPAT-like domains, respectively. The AGPATs in Cluster IV had the LPLAT_LPCAT1-like and AGPAT conserved domains.

Comparative transcriptomic analyses of orthologous AGPATs between buffalo and cattle
Prior to the comparative transcriptomic analysis of the AGPAT gene family between buffalo and cattle, we rst performed collinearity analysis between the two species. The chromosomal mapping revealed that a total of 13 and 12 AGPAT genes were found to be randomly distributed on 10 chromosomes, which are mainly located on the proximate or the distal ends of the chromosomes in the buffalo and cattle, respectively (Additional le 7: Figure S2). Collinearity analysis showed that 12 AGPAT gene pairs were orthologous between the two subspecies ( Figure 2A). Divergence times of all orthologous gene pairs between buffalo and cattle ranged from 0.465 to 2.937 Mya. All the orthologous gene pairs had Ka/Ks ratios less than 0.5 (Additional le 5: Table S5).
Using the RNA-seq data from the milk samples, we observed that three AGPAT genes (AGPAT1, AGPAT6, and LPCAT1) were highly expressed in the three lactation stages (early-, mid-, and late-lactation) in buffalo milk ( Figure 2C). By contrast, two AGPAT genes (AGPAT1 and AGPAT6) were found to be highly expressed in the same three lactation stages in cattle milk ( Figure 2B). Interestingly, the expression pattern of AGPAT1 and AGPAT6 genes existed a complementation relationship in buffalo during lactation. In cattle, AGPAT6 gene was always highly expressed during lactation, followed by AGPAT1 gene. Moreover, gene expression analysis by using qRT-PCR test also showed that both AGPAT1 and AGPAT6 were highly expressed in mammary gland tissue ( Figure 2D). It can be inferred that both AGPAT1 and AGPAT6 genes might have an important in uence on milk production performance for both buffalo and cattle.

Knockdown of AGPAT1 and AGPAT6 decreased TAG concentration in BuMEC and BoMEC
For further exploration of the potential effect of AGPAT1 and AGPAT6 on milk performance in both buffalo and cattle, we rst investigated the TGA content of BuMEC and BoMEC after the two genes silencing at 48 h post-transfection. As shown in Figure 3A, the interference e ciency of AGPAT1 and AGPAT6 in the BuMECs was 95% and 76% (P <0.0001), while their interference e ciency in the BoMECs was 89% and 82% (P <0.0001), respectively. The results suggested that these siRNA fragments could be used for further analysis. Subsequently, we observed that knockdown of AGPAT1 or AGPAT6 genes signi cantly decreased the TGA concentration in BuMECs and BoMECs ( Figure 3B and 3C; P<0.05). We found that the AGPAT1 knockdown in both BuMECs and BoMECs decreased FASN, SCD, GPAM, DGAT1, PLIN2, ACSL1, ACSS1, and LPIN1 lipogenic pathways related genes at the mRNA expression levels, while increased the expression levels of FADS2 and ACACA ( Figure 3D; P<0.05). As for the knockdown of AGPAT6 gene, we observed that ve lipid metabolism related genes (SCAP, FASN, FADS2, ACSL1, and ACACA) had higher expression levels compared to the control group, while another six genes (SCD, GPAM, DGAT1, PLIN2, LPIN2, and ACSS1) had lower mRNA expression than that of control ( Figure 3E; P<0.05).

Knockdown of AGPAT1 and AGPAT6 suppress BuMECs proliferation
The effect of AGPAT1 and AGPAT6 gene knockdown on BuMECs proliferation was investigated at three time-points (24, 48, and 72 h). Results showed that knockdown of AGPAT1 and AGPAT6 in BuMECs decreased the cell population (P < 0.05) in the culture medium compared to control group ( Figure 4A and 4C). Moreover, the mRNA levels of two proliferation-related genes (TP53 and CyclinD1) were further determined after the AGPAT1 or AGPAT6 knockdown using qRT-PCR. The results showed that increases (P< 0.05) in the mRNA expression level of the TP53 gene were observed after the AGPAT1 ( Figure 4B) and AGPAT6 ( Figure 4D) silencing. For the CyclinD1, the signi cant level was observed after the AGPAT1 ( Figure 4B) and AGPAT6 ( Figure 4D) knockdown.

Knockdown of AGPAT1 and AGPAT6 promotes BuMECs apoptosis
Using the Annexin V-EGFP Kit, we observed the apoptosis rate of the BuMECs after AGPAT1 and/or AGPAT6 gene knockdown. The results showed that the apoptosis rate was decreased in BuMECs after AGPAT1 or/and AGPAT6 silencing at 48 h ( Figure 5A, 5B, and 5C). The mRNA expression analysis showed that both AGPAT1 and AGPAT6 had lower expression levels than that of the control group ( Figure  5D; P<0.05). Moreover, expression of four marker genes (BAX, FAS, BCL-2, and Caspasse6) related to apoptosis was determined using qRT-PCR in BuMECs following treatments with siAGAPT1 and siAGPAT6. The expression level of BAX, FAS, and BCL-2 was 0.765, 0.430, and 0.619 in the siAGPAT1 group, respectively, decreasing by 23.50%, 37.00%, and 28.10% in the control group; correspondingly, the expression level of Caspasse6 gene was 1.410, increasing by 41.0% in the control group ( Figure 5E). Similar results were observed in the siAGPAT6 group ( Figure 5F) and siAGPAT1/6 group ( Figure 5G).

Discussion
In the present study, a total of 13 AGPAT genes were identi ed in the river and swamp buffaloes. The identi ed AGPAT genes were unevenly distributed on the proximate or the distal ends of 10 chromosomes in buffalo. The members of the buffalo AGPAT family were divided into four clusters based on their phylogenetic relationships. Interestingly, phylogenetic classi cation of the buffalo AGPAT family was also supported by conserved motif and gene structure analyses. Four highly conserved motifs were observed in the AGPAT family. The members in cluster I had three acyltransferase motifs and one EF-hand_1 domain. Two acyltransferase motifs were observed in Cluster II, III, and IV. The conserved protein domain analysis also supported this point. In addition, most closely related members in the same cluster of the AGPAT family harbored similar introns-exon structures (Additional le 3: Table S3).
Expression patterns of orthologous genes are often conserved and closely related to their function [30][31][32]. In the current study, we observed a total of 12 AGPAT genes orthologous between buffalo and cattle, suggesting that they might have a similar function in both species. A previous study reported that the AGPAT6 gene was the most abundant isoform in mammary gland tissue in cattle, accounting for 60% of all AGPAT mRNA, followed by AGPAT1 (18%) and AGPAT3 (10%) [13]. Another study revealed that the AGPAT6 was also highly expressed in buffalo mammary gland tissues [33]. Interestingly, our data also showed that both AGPAT6 and AGPAT1 were highly expressed in milk samples during lactation and mammary gland tissue in buffalo. These results suggested that the two AGPAT members played a dominant role in the milk fat synthesis pathway.
To further con rm this hypothesis, we explored the potential function of AGPAT1 and AGPAT6 genes on milk fat synthesis in the BuMECs and BoMECs. Knockdown of AGPAT1 or AGPAT6 genes signi cantly decreased the TAG concentration in the BuMECs and BoMECs (P<0.05). We observed that the TAG content in BuMECs and BoMECs decreased by approximately 60.8% and 61.40% after AGPAT1 knockdown, respectively ( Figure 3B). While, the knockdown of AGPAT6 resulted in an approximately 20% and 17% decrease in TAG content of BuMECs and BoMECs, respectively ( Figure 3C). The reasons for the discrepancy may be caused by the fact the AGPAT1 and AGPAT6 genes have functional differences in regulating milk fat synthesis. Our data revealed that the knockdown of AGPAT1 and AGPAT6 in both BuMECs and BoMECs signi cantly decreased in the mRNA expression levels of the fatty acid synthesis and desaturation-related genes (SCD), fatty acid activation-related genes (ACSS1), TAG synthesis-related genes (GPAM, DGAT1, and LPIN1), and lipid droplet formation-related gene (PLIN2). These results indicated that both AGPAT1 and AGPAT6 genes could decrease the TAG content by regulating the lipogenic genes. Moreover, we observed that four genes (SCAP, FASN, FADS2, and ACSL1) had differential mRNA expression levels in the BuMECs after AGPAT1 or AGPAT6 knockdown. They had upward trend expression levels after AGPAT6 knockdown compared to control group, which was opposite except for FADS2 in the knockdown treatment of AGPAT1. Evidence showed that the SCAP is a regulator of the transcription activity of FASN gene involved in fatty acid synthesis [34,35], ACSL1 plays a role in activating fatty acids destined for the TAG synthesis [36], and FADS2 plays a role in the desaturation of fatty acid synthesis [37]. More importantly, similar results were also found in the BoMECs. These results indicated that the AGPAT1 and AGPAT6 genes had functional differences in regulating the fatty acids synthesis. In other words, the AGPAT1 knockdown inhibited fatty acids synthesis, in contrast to the knockdown of AGPAT6 which promoted fatty acids synthesis. As it is well known that triacylglycerol composes one glycerol and three fatty acids, the knockdown of AGPAT6 alleviated the decrease in the TAG concentration to some extent. This explains why the reduction rate of TAG content by the treatment of AGAPT6 knockdown is lower than that of AGPAT1. In short, these ndings suggested that the role of AGPAT1 and AGPAT6 genes in determining the milk fat synthesis is not only re ected in the lipogenic genes of de novo pathway but also the effect on fatty acid synthesis.
An increasing amount of experimental evidence demonstrated that the proliferation and apoptosis of mammary epithelial cells had in uences on the development of mammary gland, milk secretion, and lactation [38][39][40]. Using the siRNA strategy, we found that both AGPAT1 and AGPAT6 knockdown could suppress the proliferation of BuMECs. Evidence showed that the stronger expression of TP53 and CyclinD1 were associated with the inhibition [41] and promotion [42] of cell proliferation, respectively. The results suggested that both AGPAT1 and AGPAT6 inhibited the BuMEC proliferation by upregulating the mRNA expression of TP53 and downregulating the CyclinD1 expression. We also observed that both AGPAT1 and AGPAT6 reduced cell viability. It is well known that the four genes (BAX, FAS, BCL-2, and Caspasse6) were associated with cell apoptosis [43][44][45]. Interestingly, our data revealed that the mRNA expression of the selected marker genes of cell apoptosis also illustrated the same results. All these results provided evidence that the AGPAT1 or AGPAT6 gene could promote the cell apoptosis of BuMECs through downregulating the BAX, FAS, and BCL-2 genes and upregulating the Caspasse6 gene at the mRNA expression level.

Conclusion
In conclusion, we identi ed 13 AGPAT genes in buffalo, of which 12 orthologous gene pairs were observed between buffalo and cattle. Our work further found that both AGPAT1 and AGPAT6 were highly expressed in milk samples in buffalo and cattle during lactation. Using in vitro assay, we recognized that both AGPAT1 and AGPAT6 genes are not just regulating the TAG synthesis in mammary epithelial cells but also affecting their growth. These ndings provided new insight about the members of the AGPAT family on how they regulate the TAG synthesis and growth of mammary epithelial cells in buffalo.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.