A genome-wide analysis of the lysophosphatidate acyltransferase (LPAAT) gene family in cotton: organization, expression, sequence variation, and association with seed oil content and fiber quality

Background Lysophosphatidic acid acyltransferase (LPAAT) encoded by a multigene family is a rate-limiting enzyme in the Kennedy pathway in higher plants. Cotton is the most important natural fiber crop and one of the most important oilseed crops. However, little is known on genes coding for LPAATs involved in oil biosynthesis with regard to its genome organization, diversity, expression, natural genetic variation, and association with fiber development and oil content in cotton. Results In this study, a comprehensive genome-wide analysis in four Gossypium species with genome sequences, i.e., tetraploid G. hirsutum- AD1 and G. barbadense- AD2 and its possible ancestral diploids G. raimondii- D5 and G. arboreum- A2, identified 13, 10, 8, and 9 LPAAT genes, respectively, that were divided into four subfamilies. RNA-seq analyses of the LPAAT genes in the widely grown G. hirsutum suggest their differential expression at the transcriptional level in developing cottonseeds and fibers. Although 10 LPAAT genes were co-localised with quantitative trait loci (QTL) for cottonseed oil or protein content within a 25-cM region, only one single strand conformation polymorphic (SSCP) marker developed from a synonymous single nucleotide polymorphism (SNP) of the At-Gh13LPAAT5 gene was significantly correlated with cottonseed oil and protein contents in one of the three field tests. Moreover, transformed yeasts using the At-Gh13LPAAT5 gene with the two sequences for the SNP led to similar results, i.e., a 25–31% increase in palmitic acid and oleic acid, and a 16–29% increase in total triacylglycerol (TAG). Conclusions The results in this study demonstrated that the natural variation in the LPAAT genes to improving cottonseed oil content and fiber quality is limited; therefore, traditional cross breeding should not expect much progress in improving cottonseed oil content or fiber quality through a marker-assisted selection for the LPAAT genes. However, enhancing the expression of one of the LPAAT genes such as At-Gh13LPAAT5 can significantly increase the production of total TAG and other fatty acids, providing an incentive for further studies into the use of LPAAT genes to increase cottonseed oil content through biotechnology. Electronic supplementary material The online version of this article (doi:10.1186/s12864-017-3594-9) contains supplementary material, which is available to authorized users.

Phosphatidic acids (PAs) are a key intermediate in the biosynthesis of TAGs and the first acylation step catalyzed by GPAT can be partly bypassed by dihydroxyacetone phosphate acyltransferase in some tissues and organisms. Therefore, LPAATs are vital to PA biosynthesis in catalyzing the incorporation of acyl groups into the sn-2 position of the glycerol backbone [13]. In higher plants, two pathways of LPAAT catalysis control the metabolic flow of lysophosphatidate (LPA) into different PAs in diverse tissues. PAs are either dephosphorylated to form TAGs, or are used in the synthesis of phospholipids, which are important components of biological membranes. Therefore, LPAATs are crucial to the biosynthesis of membrane phospholipids and storage lipids in developing seeds. In fact, in higher plants, LPAATs play an essential role in improving the fatty acid component of seeds by reducing the proportion of saturated fatty acids (SFAs).
LPAATs are thought to be among the most stringent acyltransferases in terms of substrate specificity [14][15][16][17]. Studies on the acyl specificity of LPAATs in various plant species have shown that they have stronger affinity for 16-18 carbon monounsaturated fatty acids (MUFAs) than for SFAs [18]. There are at least two classes of LPAAT genes (classes A and B) in different plant species [14]. The characteristics of class A microsomal LPAATs, as defined by Frentzen [14], are ubiquitously expressed throughout the plant, and their enzymes show specificity for 18:1-CoA. These are typical features of enzymes that biosynthesize membrane glycerolipids. Class B LPAATs were first cloned and characterized from coconut (Cocos nucifer) [19]. In this subfamily, LPAAT genes are generally expressed in seeds, and they encode enzymes with a strong substrate specificity towards unusual acyls, such as erucic acid (22:1) and lauric acid (12:0), consistent with the seed microsomal activities and oil composition of the above species. However, recent studies have indicated that the characteristics of B-class LPAATs may differ in other species. For example, the B-class LPAATs in castor (Ricinus communi) showed specificity for 18:1-OH [20].
In a recent study, a rapeseed (Brassica napus) LPAAT was shown to be the rate-limiting enzyme catalyzing the conversion of lysophosphatidate (LPA) to PA in the TAG assembly process [21]. The expression of rapeseed LPAAT genes in Arabidopsis enhanced LPAAT enzymatic activities, which redirected more acyl-CoAs from PC and relieved feedback inhibition of TAG biosynthesis. Also, the overexpression of LPAATs increased seed oil content, thus reinforcing the utility of LPAAT genes as valuable biotechnological tools. Consequently, LPAAT homologs have been cloned from poached eggplant (Limnanthes douglasii) [22], A. thaliana [23], rapeseed [24], coconut [19], corn (Zea mays) [25], Indian cress (Tropaeolum majus) [26], Java olives (Sterculia foetida) [27], and castor [20]. The overexpression of a LPAAT gene from the SLCl-1 mutant yeast stain in A. thaliana and B. napus increased the incorporation of long-chain fatty acids into the sn-2 position of TAGs, leading to an increase in seed oil content by 8-48% [28]. However, thus far, no such a transgenic work has been performed using LPAAT genes from cotton.
To date, little is known about the natural variation of the LPAAT genes and their association with seed oil production in higher plants including cotton. Furthermore, none of the genes encoding for these enzymes has been characterized in cotton; therefore, their role in increasing cottonseed oil content is unknown. In addition, since fatty acids are essential for membrane biosynthesis, saturated very-long-chain fatty acids were found to promote fiber cell elongation in cotton [29]. However, the relationship between natural genetic variation in LPAAT genes and fiber development is also currently unknown. Near-isogenic lines (NILs) or backcross inbred lines (BILs) with the same genetic background but differing in oil content, fiber initiation or elongation can provide important genetic stocks to investigate LPAAT genes with relation to fiber development and seed oil content.
The genomic organization of the LPAAT family has been mostly studied in only Arabidopsis and castor, but not in other species including cotton. Upland cotton (G. hirsutum-AD 1 ) with higher yields and wider adaptations produces approximately 95% of the world cotton fibers and seed. Another cultivated tetraploid, i.e., G. barbadense-AD 2 , accounting for the remaining world cotton production, produces much longer, stronger and finer fibers with higher cottonseed oil content and lower protein content than Upland cotton. Therefore, BILs developed from an interspecific crossing between the two species provide important genetic stocks to address if LPAATs are associated with cottonseed oil content and fiber quality. The genomes of the two cultivated tetraploid cotton species were recently sequenced [30][31][32][33]. It is commonly accepted that, like other three wild tetraploid cotton species including G. tomentosum-AD 3 , G. mustelinum-AD 4 , and G. darwinii-AD 5 , both cultivated tetraploid species originated from a common ancestor of a natural hybrid between an extant diploid A-genome cotton (likely tree cotton G. arboreum-A 2 or levant G. herbaceum-A 1 ) and a wild D genome (likely G. raimondii-D 5 ) species [34]. In this study, we identified 40 LPAAT genes in the two allotetraploid species, i.e., G. hirsutum and G. barbadense and their likely ancestral diploid species G. arboreum and G. raimondii whose genomes were sequenced earlier [35][36][37]. After an insilico analysis of the phylogeny, genomic organization and gene structure for the LPAAT gene family, we then conducted a detailed analysis of transcription profiles of the LPAAT gene family by analyzing information in two RNA-seq datasets followed by a quantitative real-time RT-PCR (qRT-PCR) analysis. Sequence variations within each gene were further evaluated for development of polymorphic markers and association analysis with cottonseed oil content and fiber quality traits. One of the LPAAT genes with a sequence variation was used in a transgenic yeast study to confirm its effect on TAGs. The study represents one of the most comprehensive genomic approaches in plants that combine an in-silico bioinformatics analysis with transcriptomics at the RNA level, SNP identification and typing at the DNA level and their associations with seed oil and fiber quality.

Genome-wide identification of LPAAT genes and phylogenetic analysis
The whole genome sequence scaffolds from the ancestral diploids G. raimondii (D 5 ), G. arboreum (A 2 ) and their descendant tetraploids G. hirsutum (AD 1 ) and G. barbadense (AD 2 ) were used for a genome-wide search for LPAAT genes in Gossypium. As a result, we identified 40 LPAATs genes in the four genomes, including 8 Gr_LPAATs based on the sequence information of D 5 reported by Paterson et al. [35], 9 Ga_LPAATs in the draft A 2 genome reported by Li et al. [37], and 13 and 10 LPAATs on the A subgenome (At-Gh (Gb)_LPAATs) or the D subgenome (Dt-Gh (Gb)_LPAATs) in the draft AD1 genome [31] and AD2 genome [33], respectively. The further phylogenetic analysis of putative cotton LPAAT proteins included some previously characterized plant acyltransferases (Fig. 1). Based on their predicted protein sequences, the LPAAT genes were classified into four subfamilies; A-class LPAAT (LPAAT2/3) (11 genes), LPAAT4/5 (12 genes), B-class LPAAT (6 genes), and LPAAT1 (plastidial) (11 genes) (Table 1). Also, some close paralogous relatives of LPAATs were identified in the four Gossypium species.th=tlb= The LPAAT2/3 subfamily, also known as the A-class LPAAT family, encodes microsomal LPAATs that show a generalized expression pattern. This subfamily includes AtLPAAT2 from Arabidopsis [23]. Within the LPAAT2/3 subfamily, Gr6LPAAT2, Ga10LPAAT2 and two copies of GhLPAAT2 from tetraploid G. hirsutum (Dt-Gh9LPAAT2 and At-Gh9LPAAT2) encoded proteins with a high level of similarity to the cacao TcLPAAT2 protein, and were identified as putative orthologs from different genomes or subgenomes, i.e., homoelogs (Fig. 1). Interestingly, three putative LPAAT2/3 isoforms were identified in Gossypium. These three LPAAT isoforms appeared to have incomplete conserved domains (Additional file 1: Figure S2). Specifically, Ga4LPAAT3 and At-Gh1LPAAT3 had only one LPAAT-characteristic Box III (ΦFPEGTR-G, where Φ is a   [38], a motif that is well conserved among proteins with the ΦFVEGTR consensus sequence in the LPAAT2/3 subfamily. Gr4LPAAT3 had Boxes III and IV (ΦPΦΦPΦΦΦ), but an imprecise VLIPRTKG consensus sequence. The sequence around Box I (Φ-NHQS-ΦDΦΦ) was quite similar among the LPAAT2/3 proteins, conforming to the consensus sequence NHXSDIDWL, where X usually indicates an arginine (R) residue [38]. However, this consensus sequence was different (NHVSDSDTΦ) in Gr7LPAAT3. The same held true for sequences around Box II (G-ΦFIDR), where there was a wellconserved E(D)YLFLER motif in most members of the cotton LPAAT2/3 subfamily (Additional file 1: Figure S2), with an exception of Gr7LPAAT3 having a different sequence (ESIFLDR). The consensus sequences of Boxes III and IV were conserved among all of the LPAAT2/3 proteins. Twelve LPAATs were clustered into the LPAAT4/5 subfamily which was a closely related sister group with LPAAT2/3 ( Fig. 1), although there were some sequence differences around various LPAAT boxes. For example, there was an aspartic acid (D) residue near Box I in the LPAAT2/3 subfamily, but a glutamic acid (E) residue near Box I in the LPAAT4/5 subfamily. This variation may explain why these sequences were grouped into two subfamilies.
Six Gossypium gene sequences designated as a paralogous LPAAT group clustered together in the B-class LPAAT subfamily, which encodes LPAATs that are usually, but not always, seed-specific and shows substrate specificities for unusual acyl groups [14]. As the sister group to the B-class subfamily the LPAAT1 (plastidial) subfamily contained eleven LPAATs in two paralogous groups from two different branches. The LPAAT1 (plastidial) subfamily includes AtLPAAT1 in Arabidopsis [39] (Fig. 1). A close relationship between B-class LPAATs and plastidial LPAATs was revealed in the phylogenetic analysis, consistent with the results of a previous report [39]. The B-class LPAATs and the plastidial LPAATs shared the same conserved motif (FPEGT) around Box III, except for Ga4LPAATB which lacked Box III (Additional file 2: Figure S3).
In plants and other organisms, LPAAT activities have been detected in the endoplasmic reticulum [23], plasma membrane [39], and mitochondria [40]. Table 1 summarizes their predicted subcellular localizations of the putative LPAAT proteins in cotton, together with information on the predicted length, molecular weight (MW), and isoelectric point (pI). Each LPAAT subfamily was predicted to a specific subcellular location. For example, within the A-class LPAAT (LPAAT2/3) subfamily, 72.7% of the LPAAT proteins were predicted to localize in the endoplasmic reticulum; and LPAAT4/5, B-class LPAATs, and LPAAT1 (plastidial) were predicted to be in the plasma membrane, the mitochondrial inner membrane, and the plasma membrane, respectively (Table 1).
Structure and domain analysis of putative LPAAT genes in sequenced diploid and tetraploid Gossypium genomes We used the GFF (Generic Feature Format) files of the four Gossypium species to analyze the exon-intron structure of putative LPAAT genes. Figure 2 shows the exon-intron structure of each subfamily. The number and location of introns varied among subfamilies, but there were some common features. In the A-class subfamily, eight deduced genes included 11 exons, but the locations of introns differed. At-Gh1LPAAT3, Ga4L-PAAT3 and Gr4LPAAT3 were closely related to LPAAT2, but had only nine, six and five exons, respectively. Within this subfamily, homoelogous LPAATs had the same structure, such as Dt-Gh9LPAAT2, Gr6LPAAT2, At-Gh9LPAAT2 and Ga10LPAAT2. Members of the LPAAT4/5 subfamily had three exons, except for Ga5L-PAAT5 and Dt-Gb13LPAAT5 which had four and two exons, respectively. These similarities in gene structure were also observed among members of the B-class and LPAAT1 (plastidial) groups (Fig. 2). The results of the domain analysis indicated the presence of a highly conserved LPLAT_LCLAT1-like domain in the A-class and LPAAT4/5 LPAATs, and an LPLAT_AGPAT-like domain in the B-class and plastidial LPAATs (Additional file 3: Table S3). In conclusion, members belonging to the same subfamilies of the phylogenetic tree had a similar gene structure and conserved motif, consistent with the results of the phylogenetic analysis.

Chromosomal distribution of LPAAT genes on Gossypium genomes
After integrating chromosomes of the sequenced cotton genomes excluding G. barbadense as its LPAAT genes were not physically well mapped, we found that the LPAAT family members were distributed unevenly among chromosomes. Among the candidate LPAAT genes, eight GrLPAATs were located in six D 5 chromosomes (D 5 /c2, D 5 /c4, D 5 /c6, D 5 /c9, D 5 /c10, and D 5 /c13), and nine GaLPAATs to seven A 2 chromosomes (A 2 /c3, A 2 /c4, A 2 /c5, A 2 /c7, A 2 /c10, A 2 /c11, and A 2 /c12) (Additional file 4: Figure S4). Because there are many syntenic gene blocks between the chromosomes of G. raimondii and G. arboreum, the dotted lines in Additional file 4: Figure S4 link homoelogous LPAATs between A 2 and D 5 chromosomes, but the chromosomes linked may not be homoelogous between genes from the two subgenomes of the tetraploid species. For G. barbadense, three GbLPAATs matched to the At subgenome, and seven GbLPAATs matched to the Dt subgenome (Additional file 4: Figure S4). For G. hirsutum, five GhLPAATs were on the Dt subgenome, and eight GhLPAATs on the At subgenome (Fig. 3).

Analysis of transcript levels of GhLPAAT genes in two RNA-seq databases of tetraploid Upland cotton
To reveal a general pattern of gene expression for the LPAAT genes identified, we analyzed the transcript profiles of LPAAT genes in two RNA-seq datasets: one with transcriptomic information for two Upland BILs, i.e., NMGA-062 vs. NMGA-105 (with G. barbadense germplasm introgression), with differing fiber lengths but similar seed oil content, and the other containing transcriptomic information for Upland Xuzhou 142 vs. its fiberless and fuzzless mutant Xuzhou 142 fl (with a likely G. barbadense origin) [41] which differed in fiber initiation and cottonseed oil content [42]. The transcript abundances of all the 13 putative LPAAT genes were determined at the fiber initiation stage (0 DPA) and the elongation stage (3-10 DPA). Of the 13 genes in the sequenced G. hirsutum genome, two (At-Gh1LPAAT3 and Dt-Gh13LPAAT5) was not expressed. The other 11 genes were transcribed at both stages of fiber development and in all tissues represented in the two RNA-seq datasets. Based on the RNA-seq data of the two BILs, the transcript levels of the A-class LPAAT genes were higher than those of the genes in the other three subfamilies in both genotypes. Interestingly, the transcript levels of Dt-Gh9LPAAT2 and At-Gh9LPAAT2 were similar in maintaining a high expression level at three different stages, but the At-Gh1LPAAT2 had a low level of expression at 0 DPA fibers ovules, with the highest level observed in 10 DPA ovules in both genotypes (Fig. 4a). There was no significant difference in expression levels between the two BIL genotypes differing in fiber length but with similar seed oil content.
Similar patterns were observed in the other RNAseq dataset for Xuzhou 142 and its fiberless and fuzzless mutant Xuzhou 142 fl. The exception was that At-Gh13LPAAT5 showed higher transcript levels in Xuzhou 142 and Xuzhou 142 fl than in the other two BILs due to different genetic backgrounds (Fig. 4b).
The results indicate that the expression of LPAAT genes at the transcription level is not involved in the natural variation of fiber length between Upland and G. barbadense and fiber initiation as controlled by the genes in Xuzhou 142 fl.
Quantitative real-time RT-PCR analysis of transcript profiles of eight GhLPAAT genes Due to a high homology for genes within the same subfamily, we chose eight of the above eleven expressed LPAAT genes in Upland cotton that have homologs with its ancestral diploid species G. raimondii and G. arboreum to design gene-specific primers (Additional file 3: Table S1) for a further evaluation of the transcript levels. The eight genes included two A-class genes, two LPAAT4/5 genes, three plastidial LPAAT genes, and one B-class gene. Using qRT-PCR analyses, we further analyzed the gene transcript levels in roots, stems, leaves, petals, ovules, and fibers of Xuzhou 142 and the Xuzhou 142 fl. We analyzed ovules at different developmental stages, i.e., −3, −1, 0, 1, 3, and 5 DPA (fibers and ovules were not separated), and 10, 15, 20, and 25 DPA (fibers and ovules were separated) as cottonseed oil accumulates rapidly in ovules after 15-20 DPA. The results showed that the eight LPAAT genes were expressed in all the organs and tissues tested throughout the plant. However, seven LPAAT genes (except for Dt-Gh6LPAAT4)

Co-localisation of LPAAT genes with quantitative trait loci (QTL) for seed oil and protein contents
To analyse if any of the GhLPAATs is genetically associated with cottonseed oil content, a co-localisation analysis of GhLPAATs with quantitative trait loci (QTL) for seed oil and protein contents was performed. We first downloaded QTL for these two traits from two types of genetic populations, i.e., intraspecific G. hirsutum and interspecific G. hirsutum × G. barbadense populations [43]. As a result, 10 GhLPAATs on five pairs of homeologous chromosomes were located within a 25 cM region of cottonseed oil or protein QTL, including At-  Gh1LPAAT2 on chromosome A01/c1, At-Gh5LPAAT1 on A05/c5, At-Gh6LPAAT4 and At-Gh6LPAAT1 on A06/c6, At-Gh9LPAAT2 on A09/c9, Dt-Gh1LPAATB on D01/c15, Dt-Gh13LPAAT5 on D13/c18, Dt-Gh5LPAAT1 on D05/c19, Dt-Gh9LPAAT2 on D09/c23, and Dt-Gh6LPAAT4 on D06/c25 except for At-Gh13LPAAT5 on A13/c13 (Fig. 3). Interestingly, the respective cottonseed oil or protein QTL co-localised LAPPT genes on six chromosomes belong to three pairs of the homeologous chromosomes, i.e., LPAAT1 on A05 and D05, LPAAT2 on A09 and D09, and LPAAT4 on A06 and D06, while other two pairs of homeologous chromosomes (i.e., A01 vs. D01 and A13 vs. D13) carry cottonseed oil or protein QTL but with different LPAAT genes. These co-localised QTL include five QTL (on A09 and D09, A13 and D13, and D06) from the Upland population only and five (on A01 and D01, A05 and D05, and A06) from both Upland and interspecific populations (Additional file 3: Table  S4). However, since a 25-cM chromosomal region may contain several hundred genes [30,31], the colocalisation of a seed oil QTL with a LPAAT gene may not indicate a causative relationship between the natural variation of the LPAAT gene and seed oil content (see next section).
Sequence variation of LPAAT genes and its association with fiber quality and seed oil and protein contents Sequence variations in the predicted genes between the sequenced G. hirsutum (TM-1) and G. barbadense (3-79 and Xinhai 21) were further analyzed. The results showed that most of the sequence variations were detected between homeologous genes from the A-and D-subgenomes, while only a few single nucleotide polymorphisms (SNPs) were predicted in homologous genes between the two cultivated species or between 3-79 and Xinhai 21 within G. barbadense (Additional file 6: Figure S6). The results imply that the differences in seed oil and protein contents and fiber quality traits between the two species are unlikely related to the natural sequence variations in the LPAAT genes, which is verified by the following correlation analysis between LPAAT gene markers and the seed and fiber quality traits.
The above RNA-seq datasets were further used to identify SNPs in the LPAAT genes within Upland cotton. In comparison with the sequenced TM-1 genome, we identified 43 SNPs in the cDNA sequences of seven LPAAT genes in the NMGA-062 and NMGA-105 dataset and 63 SNPs in five LPAAT genes in the Xuzhou 142 and Xuzhou 142 fl dataset (Additional file 3: Tables S5  and S6). However, as expected, most of the sequence differences were between TM-1 and the pair of the two BILs in the SG 747 background or the NILs in the Xuzhou 142 background (i.e., BILs vs. TM-1 or Xuzhou 142 WT/fl vs. TM-1). Only eight SNPs between NMGA-062 and NMGA-105 BILs were detected in five LPAAT genes (Additional file 3: Table S5). Between Xuzhou 142 and Xuzhou 142 fl NILs, eight SNPs were identified in four LPAAT genes, confirming the non-near-isogenic status of the two genotypes used in this study, which was recently reported by Ma et al. [41].
For an association analysis between the existence of SNPs in the LPAAT genes (detected from a comparison between the two BILs, i.e., NMGA-062 and NMGA-105) and cottonseed oil and protein contents and fiber traits, primers were designed to amplify fragments containing these SNPs using single strand conformation polymorphisms (SSCPs), which were used to screen the BIL population of 146 lines derived from a backcross between Upland SG 747 and G. barbadense Giza 75. As a result, one polymorphic SSCP marker was developed from a pair of primers (Additional file 3: Table S2) designed for At-Gh13LPAAT5, named At-Gh13LPAAT5-342 (Additional file 7: Figure S7). A transition, A/G (or T/C), was located in the 342th nucleotide of At-Gh13LPAAT5 (Additional file 3: Table S6), but the mutation appears to be synonymous with no change in amino acid. To examine the relationship of the SSCP marker with seed oil and protein contents and fiber traits, a correlation analysis was performed in the BIL population of 146 lines tested in 2006, 2008 and 2009. The marker was found to be significantly associated with both seed oil and protein content in the BILs tested in 2008. The correlation between the presence of the SSCP marker and seed oil content was significantly negative (−0.281, P < 0.01) and significantly positive with seed protein content (0.245, P < 0.01) ( Table 2). The results indicated that the presence of the SNP allele reduced cottonseed oil and increased protein content. Since seed oil and protein contents are negatively correlated (−0.905, P < 0.01), which was also reported by others [44], the opposite effects of the marker on seed oil and protein contents are expected. Moreover, the SSCP marker had a negative correlation with only one fiber quality trait, i.e., fiber uniformity (−0.176, P < 0.05) ( Table 2). However, the above results obtained in 2008 were not confirmed by the results from the same BIL population tested in 2006 and 2009, as no significant correlations were detected between the SNP and seed oil and protein contents or fiber quality (Table 2). Therefore, the inconsistent association between the marker of this LPAAT gene and seed oil and protein contents or fiber quality indicates that the natural variation of this LPAAT gene may not affect seed oil and protein synthesis or fiber development, but more studies are needed. Overall, our results indicated that sequence variations in most if not all of the LPAAT genes are not associated with the natural variation (i.e., QTL) for seed oil content or fiber quality in cotton.

Content of TAG and fatty acids in At-Gh13LPAAT5 transgenic yeast strains
Because At-Gh13LPAAT5 is not located within the QTL region for the cottonseed oil and protein QTL on chromosome A13/c13 (Fig. 3) (Fig. 6a). The analysis of fatty acid composition revealed that the expression of At-Gh13LPAAT5 in the yeast led to an approximately 25-31% increase in palmitic acid (C16:0) and oleic acid (C18:1) when compared to the nontransgenic control strains (Fig. 6b), and no significant differences were detected between the transgenic yeasts with genes from the Upland and G. barbadense sources. As compared to the empty vector transformants, the total TAGs in the transgenic strains were also increased within the range of 16-29%, and transgenic yeasts with genes from both sources also showed similar results. The qRT-PCR analysis showed that the gene expression pattern was essentially similar to the concentration of TAGs among the three transgenic yeast strains tested (Fig. 6 c, d). These results indicate that expression of the foreign gene At-Gh13LPAAT5 in the yeast can enhance oil content and selectively incorporate fatty acids into TAGs, and the gene alleles from both the Upland and G. barbadense sources has a similar effort.

Discussion
The LPAAT gene family in Gossypium The results of this study revealed details of 40 LPAAT genes in the multigene family in four Gossypium species with sequenced genomes: diploids G. raimondii and G. arboreum, and their descendant tetraploids G. hirsutum and G. barbadense. The translated protein sequences were clustered into four groups (subfamilies), consistent with other previous reports [20]. Consistent with previous studies that whole genome duplication events occurred during the evolution of Gossypium [37], we observed that each GrLPAAT gene or GaLPAAT gene in each of the diploid species corresponded to two GhLPAAT genes in the tetraploid cotton belonging to one homeologous LPAAT group. The eleven putative plastidial LPAATs in cotton showed a high similarity to the previously characterized phospholipid/glycerol acyltransferase family proteins from coconut. Proteins in the A-class LPAAT (LPAAT2/3) group, the B-class LPAAT group, and the plastidial LPAAT group were shown to have LPAAT activity [20]. However, in Arabidopsis, members of the LPAAT4/5 subfamily did not show LPAAT activity in vitro [25]. In this study, At-Gh1LPAATB, Dt-Gh1LPAATB, At-Gb1LPAATB, Dt-Gb1LPAATB, Ga4LPAATB, and Gr2L PAATB, were predicted to be in the B-class LPAAT group in cotton. Previous reports suggested that the Bclass proteins are seed-specific isoenzymes, such as LdLPAATB and CnLPAATB isolated from poached eggplant and coconut, respectively [19,22]. However, recent studies have shown that B-class proteins are not restricted to plant seeds. RcLPAATB (formerly called RcLPATB) identified in castor, was shown to express at similar levels in other tissues/organs of the plant [20]. We obtained a similar result in that there were very low transcript levels of B-class LPAATs in ovules based on our RNA-seq analyses. This is possibly because the RNA-seq data were obtained at the fiber initiation stage, and not the seed oil accumulation stage. The qRT-PCR results showed that Dt-Gh1LPAATB did not show a typical B-class expression pattern, but was expressed throughout the plant.

Expression of LPAAT genes and their natural variations in relation to seed oil content and fiber development
In higher plants, LPAATs catalyze the incorporation of acyl groups into the sn-2 position of LPAs to yield PAs, the key intermediates in the synthesis of membrane phospholipids or TAG storage lipids [11]. Each cotton fiber is a single cell that can be used as a biological model to study cell expansion. During the rapid polarized expansion of the fiber cell, abundant membrane phospholipids are required to synthesize the cytomembrane so that the fiber cell can elongate from 10 to 15 mm at 5-10 DPA to 25-30 mm at 25 DPA. Cottonseed oil accumulates later in ovules after 15-20 DPA. Therefore, in developing cottonseeds, LPAATs maybe play a dual role in that they produce raw materials for the biosynthesis of the membrane as the cell elongates and also catalyze the conversion of LPAs to TAGs. To test this hypothesis, we analyzed LPAATs in Xuzhou 142 and Xuzhou 142 fl, which have a similar genetic background based on the origin of the mutant. Xuzhou 142 fl, the natural fuzzless-lintless mutant from the fibered wild type Xuzhou 142 with both lint and fuzz, not only Fig. 6 Content of TAG and fatty acids in At-Gh13LPAAT5 transgenic yeast strains. a Gas chromatography/mass spectrometry chromatogram of fatty acids from transgenic yeasts. b The content of fatty acids in transgenic yeasts. c The total TAG content in three transgenic yeast strains from two backgrounds G. hirsutum and G. barbadense, respectively. d Expression patterns of At-Gh13LPAAT5 in three transgenic yeast strains. ** and * indicate correlation at the 0.01 and 0.05 significant levels, respectively has few or no fibers, but also has higher seed oil content than that of Xuzhou 142 [42].
The qRT-PCR results revealed different expression patterns of LPAAT genes, even those in the same subfamily. These findings suggested that LPAAT genes have diverse developmental expressions. Based on their expression patterns, the LPAATs were grouped into two categories. The first category (At-Gh6LPAAT1, At-Gh1LPAAT2, At-Gh1LPAAT3, At-Gh13LPAAT5 and Dt-Gh1LPAATB) showed high expression levels at two significant phases, i.e., during the rapid phase of fiber elongation (10 DPA) in Xuzhou 142 and at the beginning of fatty acid accumulation in ovules of Xuzhou 142 fl (20 DPA). At 10 DPA, fiber cells require membrane phospholipids to form new membranes. It is possible that these genes were preferentially expressed during membrane biosynthesis in cotton. Interestingly, these genes showed a sharp increase in transcript levels in the 20 DPA Xuzhou 142 fl ovules, but slower increases in transcript levels in Xuzhou 142 at 20 DPA ovules and during other stages of Xuzhou 142 fl ovule development. This may be because the fuzzless/lintless mutant Xuzhou 142 fl seed has few or no fibers, and therefore, does not require large amounts of membrane phospholipids for new cytomembranes of fiber development. In that case, the LPAATs would synthesize TAGs for seed oil production, rather than membrane phospholipids for developing fibers. The second group of genes (At-Gh5LPAAT1, At-Gh6LPAAT1, Dt-Gh9LPAAT1 and Dt-Gh6LPAAT4) was expressed preferentially in fiber tissues, especially Dt-Gh6LPAAT4. Their expression levels were very low in ovules of both Xuzhou 142 and Xuzhou 142 fl. This result suggested that these three LPAATs in this group may be specific to membrane biosynthesis, but more studies are needed. However, it should also be pointed out that, the expression magnitudes at the RNA level might not necessarily reflect the expression of these genes at the enzymatic level. In this study, the expression patterns of the LPAAT genes between the two genotypes were overall similar, and no SNPs of LPAATs identified and developed from a comparison between the two lines were located within the two gene regions for the fiberless and fuzzless traits [41]. We therefore conclude that none of the LPAAT genes was directly involved with the fiberless and fuzzless mutation and its higher oil content in Xuzhou 142 fl. Therefore, a segregating population between the two genotypes is unnecessary to address the role of LPAATs in the molecular mechanism of the fiberless and fuzzless trait in this mutant.
To further test if LPAAT genes are related to cotton fiber elongation, a pair of BILs with different fiber length and similar seed oil content was used. Our results indicated that the expression levels of LPAAT genes were similar between the two lines, suggesting that the genetic differences in fiber elongation between the two NILs are not associated with the expression of LPAAT genes. Furthermore, a physical mapping of LPAAT genes with QTL for cottonseed oil content showed that 10 of the 13 LPAAT genes were located within QTL regions (i.e., within 25 cM) for seed oil or protein content. However, to further test if LPAAT genes are genetically related to cottonseed oil production and fiber quality traits, a backcross inbred line population (from which the two BILs were selected) of 146 lines differing in fiber quality traits and seed oil content and fiber quality were used. The results indicated that, of all the potential SNPs identified in the LPAAT genes, only one marker from one gene At-Gh13LPAAT5 showed a significant correlation with cottonseed oil content and fiber uniformity in only one of the three field tests, suggesting that sequence differences in most if not all of the LPAAT genes are not involved in the natural genetic variation for seed oil production and fiber quality traits in cotton.
Furthermore, we designed an experiment to see if the two LPAAT variants of the same gene conferred by a SNP (which differed in a synonymous base change) had a similar effect on oil content and composition in a yeast system. Both variants were supposed to produce the protein with the same sequence and therefore a similar biochemical activity. As expected, the transgenic yeast experiment showed that the total TAG and palmitic acid and oleic acid increased compared to the nontransgenic yeast control, similar to these obtained in other transgenic plant studies [28]. The results further showed that both sources of the At-Gh13LPAAT5 gene had a similar effect, proving the lack of association between the natural variation of the LPAAT gene and fatty acid content and composition.

Conclusions
A total of 40 LPAAT genes were identified and grouped into four distinct subfamilies in four Gossypium species, i.e., i.e., tetraploid G. hirsutum-AD 1 and G. barbadense-AD 2 and its ancestral diploids G. raimondii-D 5 and G. arboreum-A 2 . The detailed analysis of the sequence variation, QTL co-localisation and content of TAG in transgenic yeasts showed that natural sequence variations in the LPAAT genes are highly limited which are unlikely associated with the natural variations in seed oil and protein contents and fiber quality traits in cotton. However, the cotton LPAATs can increase oil composition and content as demonstrated in the transgenic yeast experiment. The results provide an important lead for further studies to elucidate the involvement of LPAAT genes in the natural variation of cottonseed oil content and a possible strategy in genetic engineering for the improvement of seed oil content and composition in cotton.

Materials
Xuzhou 142 (wild type), an Upland cotton cultivar and its fuzzless/lintless natural mutant (Xuzhou 142 fl) differing in fiber initiation, and two backcross inbred lines (NMGA-062 with fiber length of 32.58 mm and NMGA-105 with fiber length of 27.06 mm) differing in fiber length, were used for RNA-seq and tissue/organ quantitative real-time RT-PCR analysis. Xuzhou142 fl and its wild type also had significantly different oil contents and fatty acid compositions [42]. Xuzhou 142 and Xuzhou 142 fl were grown at the Institute of Cotton Research (ICR), Chinese Academy of Agricultural Sciences (CAAS), Anyang, Henan province, China. On selected days post anthesis (DPA), ovules were excised from developing bolls, and fibers were scraped from the epidermis of the ovules. Petals were sampled on the day of flowering, and roots, stems and leaves were also collected from seedlings at 2 weeks after germination. All samples were quick-frozen in liquid nitrogen and stored at −80°C until use. For each time-point or tissue, three biological replications were collected.
Since the other cultivated tetraploid cotton G. barbadense contains higher seed oil content and lower protein content than G. hirsutum, introgression lines between the two species were also used in this study. SG 747 (G. hirsutum) was first crossed with Giza 75 (G. barbadense), and the resulting F 1 was backcrossed with SG 747 twice followed by three generations of selfpollination to produce advanced backcross inbred lines (BILs). The two parents and 146 BILs including NMGA-062 and NMGA-105 were grown in two replications using a randomized complete block design in Anyang, Henan in 2006Henan in , 2008Henan in , and 2009. Crop management practices and boll sampling followed local recommendations [44]. Mature seed after ginning boll samples was determined for oil and protein contents using the near infrared reflectance spectroscopy [44,45], which were also used in other cotton studies [46,47]. Another nondestructive measuring techniques, i.e., nuclear magnetic resonance [48] was not used in the current study. Fiber quality traits were measured using High Volume Instrument (HVI) 1000 by ICR, CAAS, Anyang, Henan [49].
To study the gene function of the heterologous gene, At-Gh13LPAAT5 from both Upland cotton and G. barbadense sources, the methylotrophic yeast, Pichia pastoris was used. The P. pastoris strain GS115 and plasmid pPIC3.5 K were kept in our laboratory.

Prediction and cladistic analyses of LPAAT genes
The genome sequences of G. arboreum, G. raimondii, G. hirsutum and G. barbadense were downloaded from CottonGen database (https://www.cottongen.org/home). Putative LPAATs were identified in the PFAM protein family database using HMMER software version 3.0 [50], with the LPAAT domain (PF01553, corresponding to LPLAT_LCLAT1-like and LPLAT_AGPAT-like in NCBI CDD) as the search query [51], with an initial threshold value of E ≤ 10 −20 . For the cladistics analysis of LPAAT-like proteins, the representative LPAAT protein sequences from other species including the model plant Arabidopsis, and T. cacao, which shared an ancestor with cotton at least 60 million years ago [35], were downloaded from the NCBI website (http://www.ncbi.nlm.nih.gov/guide/). Amino acid sequences were aligned using Clustal X v. 2.0.11 (http:// www.clustal.org/) under the default settings, and were further refined by a visual inspection. The alignment outputs were used to construct cladograms using the Neighbor-Joining (NJ) method, as implemented in MEGA v. 5 (http://www.megasoftware.net/). The bootstrap consensus tree was inferred from 1000 replicates.

Analysis of LPAAT genes in two RNA-seq datasets
To study the expression of LPAAT genes in ovules and fibers during fiber initiation and elongation stages, the transcriptional profiles of LPAAT genes were analyzed using information from two transcriptome sequencing databases. One database contained transcriptomic information for 0, 3 and 10 DPA ovules for the following two backcross inbred lines (BILs) with significantly different fiber lengths but similar oil and protein contents in seed: NMGA-062 and NMGA-105. The other database contained transcriptomic information for −3 and 0 DPA ovules of the cotton cultivar Xuzhou 142 and its natural fiberless mutant Xuzhou 142 fl. These databases are generated in our own laboratory and can be accessed through NCBI under accession numbers SRP038911, SRP039385, and SRP056184. Transcriptome analyses were conducted with MeV software [55].

RNA isolation and real-time PCR analysis
Based on the LPAAT coding sequences, gene-specific primers for qRT-PCR were designed with Oligo 7 software (Additional file 3: Table S1). Because we attempted to study the differential expression of LPAATs in cotton fiber and seed oil, total RNA was isolated from ovules at the fiber initiation stage (−3, −1 and 0 DPA), ovules and fibers at the early fiber elongation stage (3-10DPA), and ovules and fibers at the seed oil accumulation stage (20-25DPA). Other cotton tissues including roots, stems, leaves, and petals were also sampled and analyzed for a comparison. Each sample possessed 3 biological replicates. The RNA Prep Pure Plant kit (Tiangen, Beijing, China) was used to extract RNA. Then, 0.5 μg purified total RNA was reverse-transcribed into cDNA using the Super Script First-Strand Synthesis System for RT-PCR (PrimeScript, Takara, Dalian, China), following the manufacturer's instructions. Because of high homologies among LPAAT gene sequences within the same subfamily, 8 typical LPAAT genes from four gene subfamilies were chosen. The transcript levels of the 8 LPAAT genes were normalized to the mean value of Histone3 (AF024716) used as an internal control. All qRT-PCR reactions were run on an iCycler iQ5 Fast Real-Time PCR System (Bio-Rad, Hercules, CA, USA) according to the manufacturer's instructions. Each 20 μL reaction mixture contained 10 μL SYBR Premix Ex Taq II (2×), 0.4 μL forward and reverse primers (10 μM), 2 μL diluted cDNA, and ddH 2 O up to 20 μL. The thermal profile used for all PCRs was as follows: 10 min at 95°C for DNA polymerase activation, 40 cycles of 15 s at 95°C, 30 s at 58°C, and 20 s at 72°C, and then a 5 min elongation step at 72°C. The default settings were selected for the melting curve analysis. The gene transcript levels were calculated using the 2 -△△CT method. Three biological replicates, each with three technical replicates, were evaluated.
Identification of single nucleotide polymorphisms (SNPs) for LPAAT genes and statistical analysis with cottonseed oil and protein content and fiber quality traits The two RNA-seq datasets were also used to identify SNPs for LPAAT genes. Assembled contigs (unigenes) were scanned for SNPs using SNP detection software SOAPsnp [56]. Putative SNPs were identified to design primers (Additional file 3: Table S2) using Oligo7 for developing single strand conformation polymorphic (SSCP) markers for the BIL population per the method of Lu et al. [57]. The SSCP markers were coded as "1" (i.e., G. hirsutum or G. barbadense allele) for presence and "0" for absence of a SSCP marker and used to perform a simple correlation analysis with the seed oil and protein contents in the BIL population using SPSS software (IBM, New York, USA).

Construction of the At-Gh13LPAAT5 yeast expression vector and yeast transformation
As a SSCP marker developed from the gene At-Gh13LPAAT5 was correlated to cottonseed oil content, the gene was used in this study. First-strand cDNA synthesis was carried out using ReverTra Ace qPCR RT Master Mix (Toyobo, Japan) from Upland TM-1, Xuzhou 142, Xuzhou 142 fl and G. barbadense 3-79 used as a template for PCR to amplify At-Gh13LPAAT5 using the primer pair At-Gh13LPAAT5F/At-Gh13L PAAT5R (5'-ATTATTCGAAGGATCCATGGAAGTTCC AAGTGCGAAA-3'/ 5'-CCGCCCTAGGGAATTCTTAA GCTCCCGACATGAACC-3'), and the PCR products were sequenced in Genewiz for validation of the sequences. Because the sequences from Upland Xuzhou 142 fl (which was likely derived from a natural hybrid between Upland Xuzhou 142 and an unknown G. barbadense based on Ma et al. 2016) [41] and G. barbadense 3-79 were identical which were different from the sequence from Upland TM-1 and Xuzhou 142 (Additional file 8: Figure S1), the cDNA sequences encoding the At-Gh13LPAAT5 from Xuzhou 142 fl (to represent the G. barbadense gene allele) and Xuzhou 142 were separately cloned into the EcoRI site of the expression vector pPIC3.5 K [58,59]. 50 μL competent cells of Pichia pastoris GS115 was transformed with 10 μg of SacI-linearized pPIC3.5 K that was recombined with At-Gh13LPAAT5 by LiCl according to the Invitrogen manual.
Gas chromatography/mass spectrometry (GC/MS) profiling and statistical analysis A GC/MS analysis was performed using a gas chromatograph (7890A, Agilent Technologies, USA) equipped with a flame ionization detector (FID) and an HP-FFAP capillary column (30 m × 250 μm × 0.25 μm) by OCRI-CAAS. High purity nitrogen was used as carrier gas. Inlet pressure is 25 psi,and sample volumes of 1 μL were injected with a split ratio of 30:1 using a hot-needle technique. The GC column temperature was programmed from 150 (initial equilibrium time, 1 min) to 230°C via a ramp of 5°C/ min and maintained at 230°C for 8 min. The injection temperature was 250°C, and the test temperature was set to 280°C [60]. Transgenic yeasts were induced for 72 h and internal standard gas chromatography was performed to determine the fatty acid components and the content of triacylglycerols (TAGs). The t-test was performed using the Microsoft excel.