Identification 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 file 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 five representative mammals (Cattle, Goat, Sheep, Horse, and Human). Moreover, the motif analysis showed that a total of 10 conserved motifs were detected in the identified buffalo AGPAT genes (Additional file 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 file 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 file 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 first 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 file 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 file 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 influence 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 first 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 efficiency of AGPAT1 and AGPAT6 in the BuMECs was 95% and 76% (P <0.0001), while their interference efficiency 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 significantly 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 five 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 significant 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).