Systematic Analysis of Cotton Non-specific Lipid Transfer Protein Family Revealed a Special Group That Is Involved in Fiber Elongation

Non-specific lipid transfer proteins (nsLTPs) had been previously isolated from cotton fiber but their functions were unclear so far. Bioinformatic analysis of the tetraploid cotton genome database identified 138 nsLTP genes, falling into the 11 groups as reported previously. Different from Arabidopsis, cacao, and other crops, cotton type XI genes were considerably expanded and diverged earlier on chromosome At11, Dt11, and Dt08. Corresponding to the type XI genes, the type XI proteins (GhLtpXIs) all contained an extra N-terminal cap resulting in larger molecular weight. The research revealed that the expression of type XI genes was dramatically increased in fibers of tetraploid cotton compared with the two diploid progenitors. High-level of GhLtpXIs expression was observed in long-fibered cotton cultivars during fiber elongation. Ectopic expression of GhLtpXIs in Arabidopsis significantly enhanced trichome length, suggesting that GhLtpXIs promoted fiber elongation. Overall, the findings of this research provide insights into phenotypic evolution of Gossypium species and regulatory mechanism of nsLTPs during fiber development. HIGHLIGHT A specific group, type XI nsLTPs, was identified with predominant expression in elongating fibers of Gossypium hirsutum based on evolutionary, transcriptional, and functional analyses.


INTRODUCTION
Cotton is one of the important commercial fiber crops worldwide. Cotton fiber initiates from the outer epidermis of the ovules, followed by extensive cell elongation (about 1,000∼3,000 times) and cell wall synthesis (Basra and Malik, 1984;Kim and Triplett, 2001). In G. hirsutum, lint fibers start to develop prior to or on the day of anthesis, and fuzz fibers develop a few days later (Joshi et al., 1967;Stewart, 1975). Fiber quantity and length, to a large extent, determine cotton yield and quality of the resulting spun thread. Long-chain fatty acid (LCFA) functions as a regulator of endogenous ethylene biosynthesis in cotton ovules, which can maximize the extensibility of cotton fibers (Qin et al., 2007). However, it is unclear how LCFA is transported and accumulated on the outer epidermis of ovules to regulate fiber development. The nsLTP is one of the well-known protein families with bound lipids (Boutrot et al., 2008) such as LCFA. Detailed studies of gene expression and function of nsLTPs are paramount for addressing the above question.
The nsLTPs are widely distributed in the plant kingdom and are capable to bind to acyl groups, various phospholipids and other fatty acid groups with a broad binding affinity (Ostergaard et al., 1993;Kader, 1996;Kragelund et al., 1997;Sodano et al., 1997;Charvolin et al., 1999;Han et al., 2001). It has been revealed that most nsLTPs have a small molecular weight of about 9 kDa in higher plants (Kader, 1996). All known plant nsLTPs contain an N-terminal signal peptide which is characterized by an eight cysteine motif (ECM) backbone forming four conserved disulfide bonds that stabilize a hydrophobic cavity to enclose the lipid to shield the hydrophobic portions of the lipid (Carvalho and Gomes, 2007;Boutrot et al., 2008;Wong et al., 2017). Plant nsLTPs are deemed to be responsible for the shuttling of lipids across cytoplasm and between membranes and regulate the beta-oxidation of fatty acids in glyoxysomes and the intracellular fatty acid pools (Kader, 1996;Cheng et al., 2004). Multiple physiological and biological functions of nsLTPs have been suggested, including membrane and liposome biogenesis (Pyee et al., 1994), somatic embryogenesis (Lee et al., 2009;Chae et al., 2010;Edstam and Edqvist, 2014), pollen development (Zhang et al., 2010;Chen et al., 2011), stress resistance , defense (Molina et al., 1993;Schweiger et al., 2013), and signal transduction (Sarowar et al., 2009).
Cotton nsLTP genes have been isolated from fibers decades ago (Ma et al., 1995). GH3, Ltp3, Ltp6, and GhLTPG1 are proved to be specifically expressed in fiber cells and Ltp3 expression reaches to a maximum at the late fiber elongation stage (Ma et al., 1995;Han et al., 2013;Deng et al., 2016). Expression analysis showed that GhLtp6, GhLtp7, GhLtp8, and GhLtp11 are highly expressed during fiber initiation (Han et al., 2013). Regarded as seed trichomes, fiber initiation, and elongation are regulated by MYB genes as leaf trichomes (Guan et al., 2008;Pu et al., 2008;Machado et al., 2009;Walford et al., 2011;Huang et al., 2013;Wang et al., 2013;Liu et al., 2015). And cotton LTP3 was regulated by a MYB protein (Hsu et al., 2005). Thus it is supported that nsLTPs are involved in fiber development. However, the gene family members and the biological function during fiber development are largely unknown of these a few cotton nsLTPs cloned previously. The availability of genome sequence of G. hirsutum and insights into trichome developmental mechanism enable us to understand these scientific issues.
In Arabidopsis and rice, nsLTPs can be classified into several types (I, II, III, IV, V, VI, VII, VIII, IX, XI, and nsLTPy) based on the sequence similarity (Boutrot et al., 2008). Currently, 51 nsLTP genes have been identified in Arabidopsis (Boutrot et al., 2008;, 63 in rape and maize Wei and Zhong, 2014), 58 in sorghum (Wei and Zhong, 2014), and 52 in rice (Boutrot et al., 2008). In the present study, 138 nsLTP genes were identified in G. hirsutum that considerably expanded in Gossypium species during fiber evolution. Both transcriptional and functional analyses suggested an important role of GhLtpXIs in promoting fiber elongation.
Multiple alignments of the nsLTP protein sequences of ECMs were performed using Clustal X version 2.0 program. Phylogenetic trees were constructed with the method of Maximum likelihood or Neighbor-Joining using MEGA 6.0 in pairwise complete deletion and Amino Acid P-distance model. For statistical reliability, bootstrap tests were carried out with 1,000 replicates.
The downloaded CDS and genomic sequences of GhLtps were used to construct gene structures on Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/). Mapinspect software was used to generate chromosome location image of each GhLtp according to their starting position on chromosomes. The segmental duplication was defined as the following: (1) the length of aligned sequences cover >80% of the longer gene; (2) the identity of the aligned regions >80%; (3) only one duplication event was counted for tightly linked genes (Kong et al., 2013;Wei et al., 2013;Liu et al., 2014). Tandem duplication was defined on the basis of the criteria that tandem duplicated genes are located within 15 predicted open reading frames or within 30 kb of each other (Shin and Bleecker, 2003;Wang et al., 2011). The non-synonymous to synonymous substitution ratio (Ka/Ks) for each duplicated GhLtp gene pairs was calculated using KaKs_Calculator 2.0 software.

Expression Analysis of Cotton nsLTPs
The RPKM (Reads Per Kilobase of gene model per Million mapped reads) values of GhLtps were extracted from our RNA-seq data using ovules [0 and 5 days post anthesis (DPA)] and fibers (10-30 DPA) collected from G. hirsutum cultivars including HY405, CCRI8, and ND601 (Ma et al., 2018). Expression pattern analysis and differential expression analysis were subsequently applied with online tools Trend (http:// www.omicshare.com/tools/Home/Soft/trend) and Diffanalysis (http://www.omicshare.com/tools/Home/Soft/diffanalysis), respectively. The volcano plot was generated with Excel (Microsoft, United States) to display differentially expressed GhLtps. Transcriptome data of genome-wide gene expression in different organs and cotton species was obtained from cottonFGD (file name: fpkm.Ghir.NAU.txt.gz) and the FPKM (Fragments Per Kilobase of exon per Million fragments mapped) values of cotton nsLTPs were extracted for comparison between different tissues or organs and different Gossypium species. The heatmaps were generated with HemI version 1.0.

Plant Materials
The cotton cultivars (HY405 and ND601) were grown in the growing season in the field in Baoding, Hebei (Latitude 38 • 48 ′ N, 115 • 25 ′ E). Cotton ovules (0 DPA) and fibers (5, 10, 15, and 20 DPA) were collected at 14:00-16:00 and fast frozen in liquid nitrogen for qPCR. The cotton cultivar TM-1 were grown at 28 • C day/25 • C night with 10 h light/14 h dark cycles in greenhouse of Hebei Agricultural University and ovules and fibers were collected at 10 DPA followed by fast frozen in liquid nitrogen for gene amplification.
Arabidopsis including Col wild-type (WT) plants and transgenic plants were grown at 22 • C with 16 h light/8 h dark cycles in growth chamber at the North China Key Laboratory for Germplasm Resources of Education Ministry. Mature rosette leaves and cauline leaves were collected from flowering plants and frozen in liquid nitrogen before RNA extraction.

RNA Extraction and Quantitative Real-Time PCR Analysis
Frozen tissues were ground to a fine powder in liquid nitrogen and the total RNA was extracted with RN09-EASYspin Plant RNA purification kit (Aidlab). After quality and quantity analysis with gel electrophoresis and NanoDrop 2000, the first cDNA strand was synthesized using the PrimerScript TM RT Master Mix (TaKaRa). Then quantitative real-time PCR was performed with Fast Super EvaGreen qPCR Master Mix (US EVERBRIGHT R INC) on an ABI 7500 Real-Time PCR machine. Three biological replicates of each sample and three technical replicates of each biological replicate were used. The relative expression level of the target gene was calculated with the difference between the cycle threshold (Ct) of the target gene and reference gene ( Ct = Ct targetgene − Ct referencegene ) and corresponded to 2 − Ct . Primers used for qPCR were designed against gene specific sequences of CDS and listed in Supplemental Table S1.

Generation of Transgenic Plants Overexpressing GhLtpXIs in Arabidopsis
The full-length of GhLtpXIs coding sequences were PCR-amplified from cotton cDNA (RNA was extracted from ovules and fibers of TM-1) with primers listed in Supplemental Table S1, and the resulting products were cloned into pGreen 0229 harboring 35S promoter (Yu et al., 2004). After verification by nucleotide sequencing, the resulting vectors were transformed into Agrobacterium followed by transformation with floral dip into Arabidopsis WT (Col) plants. The transgenic plants were selected by BASTA and the first three mature rosette leaves were cut from the petioles after bolting for discoloration with gradient ethanol followed by observation of trichomes under microscope with 20× objective lens. The trichome length was measured using cellSens Standard (Olympus, Japan). Statistical analyses were performed using Excel (Microsoft, United States) with T-test (tail = 1, type = 3).

Identification and Classification of nsLTP Genes in G. hirsutum
To identify the entire collection of putative nsLTP genes in the G. hirsutum genome, local BlastP searches were conducted using 53 and 63 nsLTP protein sequences of Arabidopsis and B. rapa as queries on the whole genome of G. hirsutum. A total of 219 GhLtp genes were obtained. Each of the deduced protein sequences was analyzed with SingalP, after which 18 proteins lacking N-terminal signal sequences (NSS) were omitted. Since low molecular weight is a common character of nsLTP proteins, 46 mature proteins whose peptide length was more than 120 amino acids were not taken into consideration. Then 16 proteins lacking the Cys residues were excluded following manually scanning for the presence of the eight essential Cys residues and 1 protein without the M residue (starting residue) was removed. The remaining protein sequences were detected for the LTPAAI_LTSS structure using the Batch Web CD-Search tool. Subsequently, the GhLtp candidates were further verified according to the HMMER method. Finally, 138 genes of G. hirsutum coding the nsLTP were confirmed ( Table 1).
The sequence similarity method (Boutrot et al., 2008) was employed to classify the identified GhLtps and our results showed that the 138 GhLtps could be divided into 10 groups after multiple sequence alignments (Table 1 and Supplemental Table S2). The results showed that the ECMs were highly conserved in all of the 138 GhLtps, which could form four disulfide bonds to stabilize the tertiary structure of hydrophobic cavity (Figure 1). Different from rice, wheat, rape, and Arabidopsis, in which the majority of nsLTPs belongs to type I or II (Boutrot et al., 2008;, type XI is the largest group in G. hirsutum (Figure 1 and Table 1). The flanking amino acid residues of each CXC motif were conserved within the   members of the same group, except for type I. A more conserved consensus pentapeptide Thr/Ser-X1-X2-Asp-Arg/Lys and a more variable one Pro-Tyr-X-Ile-Ser were found in type I GhLtps as reported previously (Douliez et al., 2000;. As small molecules, nsLTPs lack supernumerary sequences at C/Nterminal or contain a short C-terminal tail. However, GhLtpXIs contained a relatively long N-terminal cap. The molecular weight (MW) and theoretical pI (isoelectric point) of each GhLtp were calculated and summarized in Table 1. Low MW is a common characteristic of plant nsLTPs and very few nsLTPs were found with a MW above 11 kDa in rice, rape and Arabidopsis. Interestingly, all GhLtpXIs displayed high MW (about 11-13 kDa). The high MWs of all the type XI members were due to the presence of supernumerary amino acid residues located at the N-terminal.

Phylogenetic and Sequence Analysis of nsLTPs
Type XI nsLPTs are rare in plants. None of them was identified in rice and wheat, and only two and three nsLTPs belong to this group in Arabidopsis and B. rapa, respectively (Boutrot et al., 2008;. On the contrary, type XI contained most members in G. hirsutum. To study the evolution of cotton nsLTP genes, this gene family was additionally identified and classified from G. arboreum, G. raimondii, Th. cacao and V. vinifera using the same method described above (Supplemental Table S3). An unrooted phylogenetic tree was subsequently built with nsLTP proteins from 8 species including G. hirsutum, G. arboreum, G. raimondii, A. thaliana, B. rapa, Th. cacao, O. sativa, and V. vinifera with neighbor-joining method (Figure 2). None of these proteins formed distinct monophyletic clusters. Type I was the largest group in plants and most sequences of this group formed a separated cluster in the tree. Cotton species have evolved a large number of type XI genes though this group is also identified in grape, cacao and rape. Two type I genes (Bra024938 and Tc11g016320) and five cotton type II genes (Ga3164g39880 and GhLtpII12-15) were close to type XI genes, indicating that cotton type XI genes maybe evolved from type I or II genes. Most importantly, only type XI sequences formed a specific cluster in the tree, indicating that these genes shared a common ancestor.
To further study the phylogenetic organization of the GhLtpXIs, a phylogenetic tree was constructed using Maximumlikehood inference from the alignment of respective 138 and 53 protein sequences of G. hirsutum and Arabidopsis. The results showed that all the nsLTPs can be divided into four clusters from A-D (Supplemental Figure S1). All the type I and III GhLtps belonged to cluster C. Type V, VIII, and IX GhLtps fell into cluster D, and cluster A contained all members of type VI and XI. Phylogenetic analysis had shown that the same type AtLtps constituted monophyletic groups except for the type II (Boutrot et al., 2008;. Consistently, GhLtpII12-15 were more distantly related to other type II GhLtps. GhLtpII13, GhLtpII14, and GhLtpII15 were more close to GhLtpXIs, while GhLtpII12 was more related to GhLtpVIs and GhLtpXI1/2/3. Gene structures of the GhLtps were obtained by comparing the predicted CDS with their corresponding genomic DNA sequences (Supplemental Figure S2). The results showed that only 16 GhLtps had introns. Among these GhLtps, six GhLtps were interrupted by two introns, and the others were interrupted by a single intron. According to the previous studies in some other plants, 35 out of 52 OsLtps, 25 out of 51 AtLtps, and 17 out of 63 BrnsLtps contain introns (Boutrot et al., 2008;. These results revealed that the percentage of GhLtps lacking introns was much higher among rice, Arabidopsis, B. rapa, and G. hirsutum, and the diverse distribution of intronic regions was quite low.

Chromosomal Localization and Gene Duplication of GhLtps
To determine the chromosomal distribution of GhLtps, the approximate position of each GhLtp was marked on the physical map of the 26 G. hirsutum chromosomes (Figure 3). It was found that the majority of the GhLtp genes located at the ends of the chromosomes. Chromosome Dt11 harbored the most GhLtps with 13 genes including 9 GhLtpXIs and none of the GhLtps located on chromosome At06 and Dt06. Since G. hirsutum is tetraploid with subgenome At and Dt, comparison was performed between each pair of homologous chromosomes. Ignoring three genes that are unable to locate, the GhLtps almost equally positioned on each genome and the number of genes on chromosomes At and its relative homologous chromosome Dt was uniform or close. Additionally, the distribution patterns of GhLtps were quite similar within each pair of homologous chromosomes, which would probably result from complex evolutional process including recombination, DNA exchanges and gene duplications after merge of genome A and D. These results showed that GhLtps distributed congruently within subgenome At and Dt.
Gene duplication events including segmental and tandem duplications were thought to be essential for the expansion of gene family in the genome (Maere et al., 2005) and thus gene duplication events of nsLTP family were investigated in G. hirsutum. A total of 32 duplication events were found, including 9 segmental duplication pairs and 23 tandem duplication pairs (Figure 3 and Supplemental Table S4). The duplication events were mainly concentrated in the same subfamily except for three pairs of tandem duplication (GhLtpII2/III3, GhLtpI28/III2, and GhLtpXI33/II12/VIII6). It was worthy to note that tandem duplication pairs distributed in the proximate locations of each paralogous blocks and the expansion of type XI group was largely contributed by tandem duplication. We subsequently calculated the nonsynonymous to synonymous substitution ratio (Ka/Ks) for each duplicated GhLtp gene pairs. Several duplication pairs of GhLtpXIs displayed larger Ks values and these genes constituted four gene clusters located in chromosomes 08 and 11 in subgenome At and Dt, implying an early divergence time of these GhLtpXIs. Additionally, most Ka/Ks ratios were below 1, except for five duplication pairs (Supplemental Table S4), suggesting that GhLtps had mainly experienced strong purifying selection pressure with limited functional divergence. These results revealed that GhLtpXIs   G. arboreum, G. raimondii, Arabidopsis, B. rapa, Th. cacao, O. sativa, and V. vinifera were used to construct the phylogenetic tree using a Neighbor-Joining method. Lines of different colors represent classification of nsLTPs. Different species were marked at the end of the lines. Greek numerals present the corresponding type of genes that cannot be displayed by lines. expanded early on chromosome At11 and Dt11 and the functions of the duplicated GhLtps did not diverge much during evolution.

Transcriptional Analysis of GhLtps
In this study, expression of GhLtps were analyzed in various organs/tissues (Supplemental Figure S3). It was worthy to point out that all the type nsLTPy genes showed specific expression in stamen, while other group members did not show any expression preference. Thirty four genes showed abundant transcription profile during fiber development and nearly 30% of them pertained to type XI.
Due to the irreplaceable economical value of cotton fibers, breeders have been making persistent efforts to develop various types of cotton (Gossypium spp.) with desirable fiber characteristics, which provides materials for investigation of the regulation of fiber development (Han et al., 2016). To further elucidate the function of GhLtps during fiber development, the expression of GhLtps in cultivars HY405, CCRI8, and ND601 was analyzed using our RNA-seq data. It was noted that 110 genes out of 138 GhLtps were expressed during fiber development (Supplemental Figure S4). The other 28 GhLtps with scarcely any transcripts covered each group. GhLtpI1/10/11 and GhLtpXI32 were highly transcribed during fiber development. And seven GhLtpXIs out of 18 GhLtps demonstrated high transcripts only in fiber initiation and early elongation stage. Further, the 138 GhLtps could be clustered into four groups based on the expression trend during fiber development (Supplemental Figure S5). The majority of GhLtps were expressed in the fiber initiation stage, and kept a relatively high expression level during the early stage of elongation. It was noticeable that the expression of nearly half of the GhLtps decreased gradually, suggesting important roles in cotton fiber initiation and elongation which are essential for lint fiber formation. Additionally, over 60% of GhLtpXIs displayed higher expression during fiber initiation and early elongation. It seemed that type XI GhLtps were important for fiber initiation and elongation. Further comparison of the transcription level was made between longer (HY405) and shorter (CCRI8 and ND601) fibers. The results between comparison of HY405 vs. CCRI8 and HY405 vs. ND601 were consistent ( Figure 4A). More GhLtps were significantly activated in cultivar HY405 with longer fibers, especially at 5, 15, and 25 DPA ( Figure 4A). Among these genes, GhLtpIs and GhLtpXIs occupied 54% proportion ( Figure 4B). The transcripts of GhLtpI26, V14, VI4, XI12, XI13, XI30, XI32, and XI33 dramatically accumulated in long-fibered cultivars when fiber started to elongate (Figures 4A,C), suggesting a possible contribution of GhLtps (especially GhLtpXIs) to fiber length improvement.
After its divergence from an ancestor shared with Th. cacao (Carvalho et al., 2011), the cotton lineage evolved into D-genome and A-genome. However, spinnable fiber only evolved in the A-genome and was further elongated after the merger of A and D genome (Paterson et al., 2012). In order to further clarify the function of GhLtpXIs in fiber development, the expression of cotton type XI genes was analyzed. Firstly, orthologs of GhLtpXIs were identified from D-genome and A-genome according to the sequence similarity, which was further confirmed by the syntenic relationships between diploid and tetraploid cotton species identified previously (Supplemental Figure S6) . Ten type XI genes in A and D genome lost during evolution and 4 pairs of non-reciprocal DNA exchanges were identified including Gr13g0054/GhAt13g0035, GhAt11g0237/GhDt11g0252, GhAt11g0234/GhDt12g1005, and GhDt12g1005/GhAt12g0916. Then the expression of orthologous gene pairs in ovules was compared at 10 and 20 DPA, respectively, using the downloaded transcriptome data. The results demonstrated that the transcriptional level of most GhLtpXIs, discarding 12 GhLtpXIs with scarce transcripts, significantly varied from their diploid progenitors (Figure 4D), among which the expression of 15 GhLtpXIs varied at 10 DPA and 12 varied at 20 DPA. It was notable that the expression of most GhLtpXIs was higher than their diploid progenitors except for GhLtpXI14-17. Especially, expression of GhLtpXI32 and GhLtpXI33 was extremely higher in tetraploid cotton than their orthologs in diploid cotton. Additionally, transcripts of GhLtpXI6 and GhLtpXI12 were dramatically abundant during fiber elongation in tetraploid cotton and GhLtpXI27, GhLtpXI28, and GhLtpXI30 increased sharply at 20 DPA.
The previous studies showed that polar lipids increased and reached the maximum level and were incorporated into the cell wall from 3 to 20 DPA (Wan et al., 2005;Edstam et al., 2013;Kumar et al., 2013). Our results suggested that GhLtps, especially GhLtpXIs play important roles during fiber elongation. Those differentially expressed GhLtps are likely to delivering lipids to the outer integument of cotton ovules (Edstam et al., 2013) and thus contribute to the variance of fiber length.

GhLtpXIs Are Involved in Fiber Development
Based on the analysis of the RNA-seq data, it is possible that the considerably expanded GhLtpXIs play important roles in different fiber developmental stages. On account of this, we selected nine genes to verify their expression by qPCR (Figure 5). The expression trend of the selected genes were consistent with RNA-seq data except for GhLtpXI12 and GhLtpXI33, which would be due to that samples at 5 DPA used for RNA-seq and qRT-PCR were different. All of these GhLtpXIs showed higher transcription levels in cultivar HY405 that produces longer fibers. The expression profile of each gene varied in different stages of fiber development, suggesting different roles of GhLtpXIs. GhLtpXI27, GhLtpXI28, and GhLtpXI30 exhibited highest transcripts at 0 DPA and decreased dramatically at 5 DPA, revealing a role in fiber initiation. The expression of GhLtpXI6/7 decreased gradually during fiber elongation. Noticeably, GhLtpXI6/7 maintained a high level in HY405 until 10 DPA when its transcription in ND601 has decreased to a relatively low level, and GhLtpXI14 and GhLtpXI32 demonstrated a constantly high level during fiber elongation in HY405, which would probably contribute to long fiber quality of the cultivar. GhLtpXI1, GhLtpXI12, and GhLtpXI33 showed a transcriptional peak during fiber elongation. The expression trend of GhLtpXI1 was consistent in both cultivars with a higher expression in HY405. Though the expression of GhLtpXI12 was comparable in two cultivars, its expression peaked earlier in HY405 at 10 DPA that is essential to fiber elongation. There is also an advanced expression peak of GhLtpXI33 in HY405 with a significantly abundant transcripts compared with that in ND601. It might be another important effecter on fiber length.
It is believed that the regulatory mechanism is shared by fiber and leaf trichome development (Kim and Triplett, 2001;Qin and Zhu, 2011;Lei et al., 2014), and Arabidopsis continues to serve as a useful experimental system for dissecting the mechanisms of cotton fiber development (Guan et al., 2014;Shangguan et al., 2016;Ma et al., 2018). Since GhLtpXIs considerably expanded in Gossypium species, several members with early divergent time were cloned and ectopically expressed in Arabidopsis to further verify the function of GhLtpXIs in fiber development. Expression analysis of T1 transgenic plants revealed that these cotton genes were successfully expressed in Arabidopsis ( Figure 6A). The trichomes of the mature rosette leaves were observed under microscope. The results showed that the trichome morphology was not affected by overexpressing of different GhLtpXIs, whereas the trichome length of all the transgenic plants was significantly longer than that of the WT plants (Figure 6B), and the improvement of trichome length was positively correlated to the expression of GhLtpXIs. These results suggested that GhLtpXIs function to promote trichome elongation. Cotton fiber which is derived from the epidermal cells is a single-celled trichome of seed and thus these GhLtpXIs might probably regulate fiber elongation.

Type XI nsLTPs Expanded Considerably in Gossypium Species
Plant nsLTPs are characterized by the ECM backbone forming a stabilized hydrophobic cavity and could be classified into 11 groups including I, II, III, IV, V, VI, VII, VIII, IX, XI, and nsLTPy according to the number of flanking amino acids within the conserved ECM domain (Kader, 1996;Carvalho and Gomes, 2007;Boutrot et al., 2008;. The nsLTP family has been identified in Arabidopsis, rice, rape, maize, and sorghum in previous studies (Boutrot et al., 2008;Wei and Zhong, 2014) and in G. hirsutum, G. arboreum, G. raimondii, Th. Cacao, and V. vinifera in this study (Supplemental Table S3). As an ancient species, V. vinifera genome contained 3 type XI genes, and this group is likely to disappear evolutionally because of the lack of type XI genes in rice and Arabidopsis. Although several type XI genes were identified in cacao, grape and rape, type I or II consistently possess most members in grape, cacao, rape, rice and Arabidopsis. It is interesting that cotton type XI nsLTPs harbor an extra N-terminal cap and are larger in molecular weight, which are different from other nsLTPs. These outcomes reveal that type XI nsLTPs specifically expanded in the Gossypium species (Figure 2).
Cotton type XI genes were close to type II genes in the phylogenetic tree and duplication pairs were found between GhLtpXIs and GhLtpIIs (Figures 2, 3 and Supplemental Table S4). The special protein structure of GhLtpXIs was only similar to GhLtpII13-15 (Figure 1). These findings imply that cotton type XI genes might diverge from type II genes. A large amount of tandem duplication events was found within GhLtpXIs (Figure 3), suggesting an essential contribution to the tremendous expansion of type XI genes.

GhLtpXIs Regulate Fiber Development
Considered as an ideal model for cell elongation, cotton fiber is the single-celled seed trichome developed from the outer integument of ovule and cotton fiber quality is quantitative traits controlled by multiple genes. Due to their economic importance and biological feature, cotton fiber development has been the subject of much scientific interest. Two decades ago, scientists have cloned nsLTPs from fibers and believed that they should play roles during fiber development (Ma et al., 1995(Ma et al., , 1997Liu et al., 2000;Orford and Timmis, 2000). Later, 12 nsLTPs were found highly expressed during fiber elongation in upland cotton fibers compared with fuzzless mutant ovules (Ji et al., 2003). However, little progress was achieved on how these small proteins affect fiber development. The complement of whole genome sequence of G. hirsutum (Li et al., 2015;Zhang et al., 2015) facilitates the study of gene families contributing to fiber development. In the present study, 138, 65, and 70 nsLTPs were strictly identified from G. hirsutum, G. arboretum, and G. raimondii, respectively. RNAseq data suggested that most GhLtps were transcribed during fiber development and some GhLtpXIs displayed much abundant FIGURE 5 | Temporal expression of nine GhLtpXIs in developing fibers. Ovules were collected on the day of anthesis and fibers were harvested on 5, 10, 15, and 20 DPA. Gene expression levels determined by qPCR were normalized to UBQ14 expression and shown as relative values to the maximal gene expression levels set at 100%. Error bars indicate SD of three biological replicates.
transcripts during fiber initiation and elongation (Figure 5 and Supplemental Figure S4). Significant transcriptional differences of GhLtps were found between cultivated upland cotton varieties with longer and shorter fibers. And many GhLtpXIs displayed preponderance of expression in longer fibers, which was likely to contribute to the variation of fiber quality (Figures 4, 5). The transcription of several GhLtpXIs was verified by qPCR with a higher level or advanced peak in long-fibered cultivar ( Figure 5). Furthermore, significant elongation of Arabidopsis leaf trichomes induced by overexpression of GhLtpXIs suggested functional roles of GhLtpXIs in promoting cell elongation (Figure 6).
LCFA regulates endogenous ethylene biosynthesis in cotton ovules to promote the extensibility of fibers (Qin et al., 2007). The basic function of nsLTPs is to transfer lipids. GhLtps were speculated to bind with various lipids to be responsible for the intracellular and intercellular movement of lipids including fatty acids and thus affect ethylene production and fiber elongation. A GPI-anchored lipid transport protein (GhLTPG1) has been identified to regulate cotton fiber elongation through mediating the transport of phosphatidylinositol monophosphates (Ostergaard et al., 1993;Deng et al., 2016), which supports the speculation. As a consequence of the transportation of lipids, the intercellular fatty acid pools will be regulated by GhLtps, which would cause a feedback regulation on reactions enrolling fatty acids, including LCFA biosynthesis and accumulation. In addition to transferring lipids, nsLTPs can also act as lipid sensors and lipid chaperones (Wong et al., 2017). Plant nsLTPs are supposed to bind a lipid molecular competing with elicitin and mediate signaling (Blein et al., 2002;Maldonado et al., 2002). When nsLTPs present part of a lipid (typically the hydrophilic head group) to another protein, they act as lipid chaperones to mediate signaling (Wang et al., 2005). Thus a hypothesis that GhLtps participant in the signaling of phytohormones involved in fiber development and other signaling pathways mediating fiber development could be proved when functional lipids were isolated as signaling molecules from fibers or the outer integument of ovules.

GhLtpXIs Might Comprise a Possible "Fiber Clade"
Shortly after the divergence from the same ancestor as Th. cacao at least 60 Myr ago (Carvalho et al., 2011), the genus Gossypium experienced an abrupt polyploidization with a maternal Agenome propagule resembling G. herbaceum and a pollen parent D-genome species resembling G. raimondii which diverged ∼5-10 Myr ago (Senchina et al., 2003). The nascent AtDt allopolyploid diverged into at least five species, of which two major cultivated species G. hirsutum and G. barbadense were domesticated independently to spawn textile industry and became a major oilseed. Although G. hirsutum has been well domesticated and bred for many cultivated varieties that vary in fiber quality, the evolution of fiber and regulatory mechanism of fiber development are unclarified. Our results suggested that the expanded type XI genes in G. hirsutum genome are involved in fiber development. Although GhLtpXIs were phylogenetically close to GhLtpIIs, they evolved to be different in protein structure and MWs. Members of this group contain an extra N-terminal cap (Figure 1). Protein sequences of GhLtpXIs were highly conserved within the ECMs and only varied at the N-terminal (Figure 1). Additionally, the divergence time of GhLtpXIs was earlier than other GhLtps. It is clear that GhLtpXIs form a unique group distinctive from other GhLtps.
The sequence of a G. hirsutum cultivar reveals many nonreciprocal DNA exchanges between subgenomes that may have contributed to phenotypic innovation that spinnable fiber evolved from the merger of A genome and D genome determining a fibered and a fibreless phenotype, respectively (Paterson et al., 2012). Noticeably, non-reciprocal DNA exchanges were found in cotton type XI genes and most transcribed GhLtpXIs in ovules displayed significant expression difference during fiber elongation compared with their orthologs in A or D genome, indicating a correlation between cotton type XI genes and fiber evolution. Therefore the cotton XI genes are speculated to comprise a possible "fiber clade" (Paterson et al., 2012) evolving from polyploidy, non-reciprocal DNA exchanges and duplication events, which would be supported by further sequence comparison of nsLTPs in different cotton species and their evolutional close species and molecular investigation on how these candidates regulate fiber development.

AUTHOR CONTRIBUTIONS
ZM, CM, and YY designed the experiments. CM and YY conducted the experiments and analyzed the data. ZL, LC, and YZ performed part of the experiments. XL helped in the bioinformatics analyses. LW and GZ helped in growing experimental materials. CM and YY wrote the manuscript. ZM and XW carefully edited the manuscript and analyzed the data. All authors revised the manuscript.

ACKNOWLEDGMENTS
We thank Xiyin Wang and Jinpeng Wang for their assistance in bioinformatics analysis. Our thanks also go to Hongping Wu for her help in data analysis. This work was supported by grants from the National Key Research and Development Program of China (2016YFD0101006), the Science and Technology Support Program of Hebei Province (16226307D), the Top Talent Fund of Hebei Province, the High School Science and Technology Research Program of Hebei Province (ZD2018057) and financial support for research from College of Agronomy of Hebei Agricultural University.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018. 01285/full#supplementary-material Supplemental Figure S1 | Phylogenetic tree of GhLtps. The amino acids of the mature proteins of cotton and Arabidopsis nsLtps were used to build the phylogenetic tree using Maximum likelihood method. The 10 types of nsLTPs and another five AtLtps are showed with different color background. Four clusters are marked against different colors.
Supplemental Figure S2 | Gene structure of the GhLtps. Exons are represented by green boxes and introns are showed with black lines.
Supplemental Figure S3 | Expression profiles of GhLtpXIs in different organs and tissues of G. hirsutum. Data shown were log 2 -transformed FPKM of each gene which was quantified using RNA-seq data downloaded from CottonFGD. Clustering of expression level is showed on the right. The color bar represents the relative expression level.
Supplemental Figure S4 | Expression of GhLtps during fiber development. Data shown were log 2 -transformed RPKM. Developmental stages of fiber are indicated in days post anthesis (DPA) above. The color bar represents the relative expression level.
Supplemental Figure S5 | Expression trend analysis of GhLtps. Number of genes is indicated below the relevant expression pattern.
Supplemental Figure S6 | Phylogenetic tree of nsLTPs in cotton species. The full lengths of mature protein sequences were used to construct the phylogenetic tree using a Neighbor-Joining method. The same color was used to present orthologous gene pairs and non-reciprocal DNA exchanges were labeled with different dots.
Supplemental Table S1 | List of primers used in this study.
Supplemental Table S2 | Diversity of the eight cysteine motif.
Supplemental Table S3 | Classification of nsLTP genes in cotton, Arabidopsis, rape, rice, cacao and grape.
Supplemental Table S4 | Ka/Ks analysis for the duplicated gene pairs.