Introduction

Lipids play a vital role in plant growth and development. They can maintain cell function and mediate cell signaling associated with stress responses. Plant non-specific lipid transfer proteins (nsLTPs) (6.5–10.5 kDa in size) are able to transfer phospholipids and fatty acids between membranes in vitro1,2. They are abundantly present in various plants, representing up to 4% of the total soluble protein2. These peptides are structurally characterized by an eight cysteine motif (8 CM) backbone with the general form C-Xn-C-Xn-CC-CXC-Xn-C-Xn-C, and the cysteine residues are linked by four disulfide bonds to stabilize a tertiary structure of a hydrophobic cavity2,3. Moreover, almost all nsLTPs carry an N-terminal secretory signal peptide (21–27 amino acids in length) in their nascent polypeptides, indicating that they are secreted proteins2,4. In previous studies, plant nsLTPs have been identified to participate in a range of biological processes, including cuticular wax and cutin synthesis, biotic and abiotic stress resistance, plant signaling, seed maturation and sexual reproduction2,4. In some cases, nsLTPs are considered pathogenesis-related (PR) proteins and constitute the PR-14 family5. The defensive role of plant nsLTPs against fungal, bacterial and viral pathogens has been well verified by ample evidence2,6,7,8.

Plant nsLTPs have been well studied and identified in many plant species including both flowering and non-flowering plants2. They were initially classified into two types according to molecular mass: nsLTP1 (9 kDa) and nsLTP2 (7 kDa)1. The disulfide bond linkages of nsLTP1 at Cys1-Cys6 and Cys5-Cys8 differ from those of nsLTP2 at Cys1-Cys5 and Cys6-Cys8. Recently, a new classification according to sequence similarity and intervals of 8 CM was proposed by Boutrot et al.3. The system categorized 267 nsLTPs into nine types (Type I-IX) based on a genome-wide analysis of rice (Oryza sativa), wheat (Triticum aestivum) and Arabidopsis thaliana. Then, nsLTPs in other plant species were also grouped according to Boutrot’s method with slight modifications4,9,10,11,12, and novel types such as Type X9 and Type XI4 were identified.

Cotton provides the world’s most utilized natural fiber for the textile and garment industries. The genus Gossypium comprises approximately 45 diploid species and can be divided into eight monophyletic groups (each designated alphabetically as A through G, and K)13,14. The A- and D- genome diploids diverged from the same eudicot progenitor approximately 5–10 million years ago (MYA). Then, ancient hybridization between A and D diploids occurred, resulting in the generation of a clade of five allotetraploid species approximately 1–2 MYA15. Gossypium hirsutum is one of the descendant allotetraploid species and may be derived from a spinnable fiber capable A genome species (G. arboreum) and a non-spinnable fiber capable D genome species (G. raimondii)15,16. Cotton polyploidization has altered the function of the A and D ancestral diploid genomes and conferred emergent properties17. More than 90% of world cotton production is attributable to the cultivation of G. hirsutum15. However, a soil-borne vascular disease, Verticillium wilt, which is caused by Verticillium dahliae, is the most prevalent and lethal disease affecting cotton production. It is a critical plant quarantine disease and called the “cancer” of cotton crop18. Few G. hirsutum germplasms were found to be resistant to V. dahliae19. Comprehensive prevention and control measures are limited solely to improvements in genetic resistance19,20.

In previous studies, only a small portion of nsLTPs from cotton have been characterized2,21,22. With G. raimondii, G. arboreum and G. hirsutum genome sequencing completed14,15,23, an excellent opportunity is coming to initiate whole-genome annotation and to perform comparative genomic study in Gossypium. Until now, a genome-wide overview of the nsLTP family in Gossypium has yet to be reported. Thus, a systematic molecular evolution and expansion analysis of the nsLTPs in Gossypium is urgently required. In this study, putative nsLTPs were identified in G. raimondii, G. arboreum and G. hirsutum. We conducted a comprehensive study on the phylogenetics, genomic structure, chromosomal localization and expansion history to evaluate the molecular evolution, expansion history and expression profiles of the nsLTP family in Gossypium.

Results

Genomic identification of putative nsLTPs in G. raimondii, G. arboreum and G. hirsutum

The availability of G. raimondii, G. arboreum and G. hirsutum genome sequences makes it possible to identify all the nsLTPs in the three Gossypium species. The BLASTP program was utilized to search for candidate nsLTPs in cotton with the query sequences from Arabidopsis. Initially, 104, 104 and 182 protein sequences were identified in G. raimondii, G. arboreum and G. hirsutum, respectively. After screening and selection as described in the materials and methods, a total of 189 nsLTPs were confirmed and described (Table S2). Among them, G. raimondii and G. arboreum contain a similar number of nsLTPs (47 and 51, respectively), despite the fact that G. raimondii has a much smaller genome size (885 Mb/1 C) than G. arboreum (1,746 Mb/1 C). In G. hirsutum (2,173 Mb), 91 nsLTPs were identified, representing almost a two-fold increase over the number of nsLTPs in its diploid progenitors. We designated the genes identified in G. raimondii, G. arboreum and G. hirsutum as GrLTPs, GaLTPs and GhLTPs, respectively.

Phylogenetic analysis of the nsLTP family

To determine the evolutionary relationships of nsLTPs, phylogenetic analysis of the identified nsLTPs in Gossypium and Arabidopsis was completed with the MrBayes and PHYLIP tools (Fig. 1, Fig. S1). There were similar results with high support values from each method. According to Boutrot’s classification system, the nsLTP family in Gossypium was divided into 8 subfamilies (Type I, II, III, IV, V, VI, VIII and IX), and no Type VII nsLTPs were identified in cotton (Fig. 1, Fig. S1, Table S2). The member proportion was different in each subfamily (Fig. S2a). The Type I subfamily (33.33%) contained the most members, followed by Type II (23.28%), Type V (16.93%) and Type IV (11.64%). The least represented subfamily was Type IX (1.59%). A similar member distribution in each subfamily was found in each Gossypium species (Fig. S2b). Besides, the proportion of nsLTPs in Type I was 35.29% and 38.30% in G. arboreum and G. raimondii, respectively, whereas it decreased to 29.67% in G. hirsutum. The proportion of members in both Type II and Type VIII in G. hirsutum (25.27% and 7.69%, respectively) was higher than that in G. arboreum (23.53% and 3.92%, respectively) and G. raimondii (19.15% and 4.26%, respectively). Moreover, not all the subgroups were present in each Gossypium species, and no Type III and Type IX nsLTPs existed in G. arboreum.

Figure 1
figure 1

Phylogenetic relationships of the nsLTP family from G. raimondii, G. arboreum and G. hirsutum.

The phylogenetic tree was generated using the Bayesian method based on multiple alignments of the nsLTP sequences. The numbers in the clades are posterior probability values. The nsLTP subfamilies are indicated by different colors.

Full size image

Additionally, all the 189 nsLTP sequences were evaluated with OrthoMCL clustering. With the default stringency, 21 orthologous groups (OGs) in 8 subfamilies were shown (Table S3). Each subfamily contained one or more OGs, and each distinct OG was shared by only one specific subfamily. The OG distributions were similar to the phylogenetic classifications of the nsLTP family in the three Gossypium species. The protein structures were highly diverse in all the identified nsLTPs (Table S1). The amino acid lengths of the nsLTPs in the Type I, Type V and Type VIII subfamilies were relatively longer, while the proteins in Type II and Type III had relatively shorter amino acid lengths. A similar distribution in the molecular weight of the nsLTPs also existed.

Conserved protein motifs and exon/intron structure of nsLTPs

The main characteristic of plant nsLTPs is the presence of 8 CM in highly conserved positions. Using the WebLogo program, the sequence logos of the identified nsLTPs in Gossypium were generated to further confirm the conservation of amino acid residues (Fig. 2). Moreover, a variable number of inter-cysteine amino acid residues was displayed through multiple alignments, and 8 nsLTP types were therefore identified based on the typical spacing of the 8 CM (Table 1). Between the conserved Cys1 and Cys2 residues, Type II and VIII nsLTPs contained 6–8 residues, while Type V and IX contained 14 and 13 residues, respectively. Between the conserved Cys4 and Cys5 residues, Type I nsLTPs contained 19-20 residues, while the other types contained relatively fewer residues (8–12). Between the conserved Cys6 and Cys7 residues, Type III nsLTPs contained 12 residues, while the other types contained more residues (19–27). Between the conserved Cys7 and Cys8 residues, Type I nsLTPs contained more residues (13–15) than Type III and Type IX nsLTPs (6).

Table 1 Diversity of eight cysteine motifs in different types of nsLTPs identified in Gossypium.
Full size table
Figure 2
figure 2

Conserved domain analysis of the nsLTPs using the WebLogo program.

The height of the letter designating the amino acid residue at each position represents the degree of conservation. The numbers on the x-axis represent the sequence positions in the corresponding conserved domains. The y-axis represents the information content measured in bits.

Full size image

The MEME motif search tool was employed to identify 20 distinct conserved motifs in the nsLTPs (Fig. 3). Based on the distribution of the predicted motifs, the 189 nsLTPs were categorized into 8 distinct subfamilies, which was consistent with the classification from the phylogenetic analysis. All the nsLTPs had either Motif 4 or 3, which represented Cys2. Motif 13, 11 and 14 corresponded to Cys1, and Motif 7, 18, 19, 6 and 15 mapped to Cys3Cys4. In addition, Motif 8, 9, 5, 1 and 12 corresponded to Cys7 and Cys8. Motif 1, 2, 6, 10 and 15 mapped to Cys5XCys6, and different residues were found in the central position of the Cys5XCys6 motif. Seven hydrophilic residues (Arg, Gly, Glu, Asp, Gln, Ser and Lys) and five hydrophobic residues (Tyr, Phe, Leu, Val and Met) existed at the X position of the Cys5XCys6 motif in the 189 nsLTPs. Moreover, some subfamily-specific motifs were identified. For instance, Motif 5, 7 and 10 were only present in Type I nsLTPs, and Motif 15 only existed in Type VI nsLTPs.

Figure 3
figure 3

Motif compositions (left panel) and gene structures (right panel) of some nsLTP members in Gossypium.

The conserved domains were identified using the MEME web server. Different motifs are represented by various colored boxes. The location of each motif can be estimated using the scale at the bottom. Gene structures of nsLTPs were predicted with the GSDS software. The exons are indicated by yellow boxes and the introns are indicated by black lines.

Full size image

The gene structure of the nsLTP family was also investigated (Fig. 3). Our results revealed low diversity in the distribution of intronic regions amid the exonic sequences. The intron patterns, formed by relative position and phase, were highly conserved within each phylogenetic group. The number of introns per gene varied from 0 to 2. No introns were identified in Type II and Type III nsLTPs, while some nsLTPs in other subfamilies were interrupted by 1–2 introns positioned 15 to 78 bp downstream of the codon, which almost encoded the Cys8 in the 8 CM.

Chromosomal localization and duplication of nsLTPs in cotton

Based on the genomic location of nsLTPs (Table S2), chromosomal distribution diagrams of the nsLTP family were generated for the three Gossypium species (Fig. 4). In G. arboreum, the 51 nsLTPs were distributed unevenly on 12 chromosomes, and the number of nsLTPs on each chromosome varied widely (Fig. 4a). Chromosome 6 had the largest number of nsLTPs with 11 members, followed by Chromosome 1 with 9 genes. In contrast, only 1 gene was located on Chromosome 5, and no nsLTPs were found in Chromosome 2. In G. raimondii, the 43 nsLTPs were located on 13 chromosomes, with Chromosome 8 containing the largest number of genes (9) (Fig. 4b). Two genes were present on Chromosome 1, 2, 3, 5 and 9, whereas only a single nsLTP was localized on Chromosome 4, 10 and 12 each. In addition, a total of 30 and 35 nsLTPs were mapped in the 9 At and 10 Dt subgenome chromosomes of G. hirsutum, respectively (Fig. 4c). Dt- Chromosome 6 contained the maximum number of nsLTPs with 9 genes, followed by At- Chromosome 2, Dt- Chromosome 1 and Dt- Chromosome 8 with 6 genes each, while only one gene was present on At- Chromosome 13, Dt- Chromosome 5, Dt- Chromosome 7 and Dt- Chromosome 10 each. Besides, no nsLTPs were located on At- Chromosome 1, 3, 10 and 12, and Dt- Chromosome 3, 4 and 12. Moreover, several nsLTP clusters were detected on chromosomes such as the top of Dt-Chromosome 6 and the bottom of At-Chromosome 2 in G. hirsutum.

Figure 4
figure 4

Chromosomal locations of nsLTPs from G. arboreum (a) G. raimondii (b) and G. hirsutum (c).

The scale represents megabases (Mb). Chromosome numbers are indicated above each vertical bar. The markers before the gene names indicate the nsLTP subfamily. The duplicated gene pairs are joined by black lines.

Full size image

Gene duplication events were investigated to elucidate the expansion pattern of the nsLTP family in Gossypium. In our study, 7, 5 and 6 duplicated gene pairs were identified in G. arboreum, G. raimondii and G. hirsutum, respectively (Table 2). In G. arboreum and G. raimondii, the duplication events were concentrated in similar subfamilies (including Type I, II and V), whereas in G. hirsutum the duplication events occurred in Type II, IV, V and VIII. This result revealed the strong expansion preference for some nsLTP subfamilies in Gossypium. Meanwhile, based on the sequence analysis and the chromosomal distribution, 2 and 2 paired genes in G. arboreum and G. hirsutum, respectively, were identified to be involved in segmental duplication events, while the other 14 pairs were linked to tandem duplication events.

Table 2 Ka/Ks analysis for the duplicated gene pairs.
Full size table

In addition, the Ka/Ks ratio of each duplicated gene pair was calculated to assess the molecular evolutionary rates (Table 2). In this study, 83.3% of the duplicated nsLTPs from cotton have mainly undergone purifying selection pressure after the duplication events (Table 2). As purifying selection apparently constrains the divergence of the duplicated genes, the functions of the duplicated nsLTPs might not diverge much during subsequent evolution. Moreover, the divergence times between the duplicated gene pairs were analyzed. In G. arboreum and G. raimondii, most of the Ks values were less than 0.17, and their corresponding duplication age might be less than 32.69 million years ago (MYA). Only one duplication event (GaLTP47/51) occurred approximately 194.74 MYA. In G. hirsutum, the Ks values were 0.03–0.13, and the duplication ages were estimated to be 5.84–24.50 MYA.

Orthologous gene analysis and synteny block detection

The evolution of nsLTPs was analyzed among the three Gossypium species. We identified 23 and 24 orthologous gene pairs within the At and Dt subgenomes, respectively, in G. hirsutum and their corresponding ancestral A and D diploid genomes (Fig. 5, Table S4). Of the orthologous gene pairs, most were distributed in Type I, II, IV and V (Fig. S3, Table S4). Besides, 24 GrLTPs have no orthologs in the Dt subgenome of G. hirsutum, and 31 GaLTPs have no orthologs in the At subgenome. Meanwhile, the Ks values of the orthologous gene sets (At vs A and Dt vs D) were also compared. The divergence times between the allotetraploid and its progenitor genomes were estimated to be 3.46–3.54 MYA (Ks peaks at 0.0180 and 0.0184, respectively) (Fig. S4). The similar divergence time suggested that the nsLTPs in At and Dt subgenomes might have evolved symmetrically.

Figure 5
figure 5

Distribution of nsLTP orthologous gene pairs and synteny block identification between the At subgenome in G. hirsutum and the diploid A genome in G. arboreum genome (a) and the Dt subgenome in G. hirsutum and the diploid D genome in G. raimondii (b).

The differently colored lines represent the subfamilies within the nsLTP family.

Full size image

The conservation and rearrangement of certain chromosome segments may play an important role in the adaptive evolution of the three Gossypium species. In our study, we discovered 3 conserved syntenic blocks between the At subgenome in G. hirsutum and the A diploid genome in G. arboreum, and 1 conserved syntenic block between the Dt subgenome in G. hirsutum and the D diploid genome in G. raimondii. This result indicated that the collinear relationships of the nsLTPs between the At subgenome and the G. arboreum genome were higher than those between the Dt subgenome and the G. raimondii genome. These conserved segments contain different numbers of nsLTPs, ranging from 2 to 3. According to the analysis, the genes situated on G. hirsutum At Chromosome 11, 4 and 5 were predicted to have corresponding orthologs on G. arboreum Chromosome 1, 10 and 7, respectively, whereas the genes on G. hirsutum Dt Chromosome 11 were found to originate from the conserved block on G. raimondii Chromosome 11 (Fig. 5). In particular, six paired nsLTPs (GhLTP64/GaLTP49, GhLTP65/GaLTP50, GhLTP38/GaLTP20, GhLTP39/GaLTP19, GhLTP82/GaLTP22 and GhLTP83/GaLTP22) were located in genomic regions with synteny between the G. hirsutum and G. arboreum genomes, while three paired nsLTPs (GhLTP32/GrLTP42, GhLTP33/GrLTP41 and GhLTP34/GrLTP40) were located in genomic regions with synteny between the G. hirsutum and G. raimondii genomes.

Analysis of orthologous gene expression patterns in cotton

In the current study, expression patterns of the orthologous genes were compared in G. raimondii, G. arboreum and G. hirsutum (Figs 6 and 7 and Fig. S5). Among the intraspecific orthologous genes, 5 paired genes (GaLTP47/51, GaLTP49/50, GhLTP10/86, GhLTP53/74 and GrLTP40/41) shared almost equivalent expression patterns in the tested tissues and leaves inoculated with V. dahliae (Fig. 6). However, this was not the case for GhLTP82/83, GrLTP13/14 and GrLTP21/22. The expression patterns of these duplicated gene pairs were strongly divergent. GhLTP83 and GrLTP14 showed predominant expression in stems compared with GhLTP82 and GrLTP13, respectively. GrLTP21/22 also showed reciprocal silencing in roots and stems. Apart from the tissue-specific expression biases, GhLTP82/83 and GrLTP13/14 showed similar expression profiles after V. dahliae infection, while GrLTP21/22 displayed reverse expression pattern in response to V. dahliae.

Figure 6
figure 6

Quantitative RT-PCR analysis for selected duplicated nsLTPs from G. arboreum, G. raimondii and G. hirsutum in different organs (a–h) and leaves after V. dahliae treatment (i–p).

Full size image
Figure 7
figure 7

Expression patterns of the orthologous gene pairs from G. arboreum, G. raimondii and G. hirsutum in different organs (a–h) and leaves after V. dahliae treatment (i–p).

Full size image

Among the interspecific orthologous genes, 10 pairs between G. hirsutum and G. arboreum and 6 pairs between G. hirsutum and G. raimondii were selected to evaluate their expression relationship (Fig. 7, Fig. S5). Among them, the 4 paired genes (GaLTP25/GhLTP69, GaLTP35/GhLTP42, GaLTP36/GhLTP86 and GaLTP36/GhLTP10) shared similar expression patterns, and both genes were up-regulated after V. dahliae infection. However, the other 12 pairs displayed strong expression bias towards one copy in different organs, and the responses to V. dahliae were drastically different (Fig. 7, Fig. S5).

Expression profile of nsLTP members in cotton

In the current study, a total of 78 genes (16 GaLTPs, 19 GrLTPs and 43 GhLTPs) from all the subfamilies were selected, and the expression patterns of these genes varied significantly in different tissues (Fig. 8). Generally, Type V nsLTPs were mainly expressed at high levels in roots and young stems, while nsLTPs in Type I, II, IV and IX were predominantly expressed in stems and leaves. Additionally, nsLTPs play an important role in the protective mechanisms against pathogens in plants. The nsLTPs from some plants have reported to show obvious antifungal activities against V. dahliae24,25,26,27, we thus investigated the expression patterns of the nsLTPs in response to V. dahliae to take a further step towards understanding the function of cotton nsLTPs against V. dahliae. As displayed in Fig. 8, most nsLTPs in Type I and V showed a marked increase in transcript levels after V. dahliae inoculation, and different genes responded to different stages of V. dahliae infection. Some Type II nsLTPs were also slightly up-regulated in response to V. dahliae stress, while some nsLTPs in Type III, IV, VI, VIII and XI showed a decline in expression levels after treatment with V. dahliae.

Figure 8
figure 8

Heat map representation and hierarchical clustering of GaLTPs, GrLTPs and GhLTPs across different tissues (a–c) and in leaves after V. dahliae treatment (d–f).

The color bar represents the relative signal intensity value.

Full size image

Discussion

The categorization of nsLTPs based on phylogenetic clustering provided comprehensive information about the gene family and facilitated further functional analysis. In the current study, the 189 nsLTPs identified in G. raimondii, G. arboreum and G. hirsutum were classified into 8 subfamilies (Type I, II, III, IV, V, VI, VIII and IX) through phylogenetic analysis, and no Type VII nsLTPs were found in Gossypium (Fig. 1, Fig. S1). It has been speculated that Type VII nsLTPs appeared specifically in monocots while Type IX nsLTPs were unique to dicots3,9, and our results in Gossypium further confirmed this viewpoint. In addition, G. arboreum lost all the Type III and Type IX nsLTPs (Fig. S2), suggesting that the evolution of plants not only involves gene retentions, but also is accompanied by gene losses and mutations4,28. Meanwhile, the proportion of nsLTPs in each subfamily indicated that Type I seemed to have contracted while Type VIII expanded in G. hirsutum compared with its diploid progenitors (Fig. S2b). The gene retentions and losses might be associated with the related functions during plant evolution29.

As for multigene families, the analysis of gene expression profiles often provides useful clues for functional assessment. The tissue-specific expression patterns of nsLTPs (Fig. 8) indicated their important roles in performing diverse developmental and physiological functions in cotton. Besides, Type I nsLTPs were reported to be involved in plant defense against phytopathogens4,10,30, and were classified as PR-14 family. Such a role is particularly consistent with our results (Fig. 8). Most Type I nsLTPs showed increased expression levels after V. dahliae infection. In addition, many Type V nsLTPs were also significantly induced after V. dahliae attack, suggesting that they may participate in the pathogen response. In a previous study, the expression of Type V nsLTPs was observed primarily in the vascular bundles, and they were deduced to be involved in defense signaling10. Additional work is needed to further characterize the role of Type V nsLTPs and their functional involvement in signal transduction. Moreover, the nsLTPs in other subgroups showed no expression changes or a decline of expression levels when treated with V. dahliae, suggesting that they may participate in other biological processes or respond to the attacks of other types of pathogens3,4,30,31,32.

Most nsLTPs within each subgroup shared common motif compositions, and the motif distributions among the different cotton species did not vary much (Fig. 3). Motif 5 and 7 corresponded to the highly conserved residues located in Pro-Tyr-X-Ile-Ser and Thr/Ser-X1-X2-Asp-Arg/Lys, respectively, and the two consensus pentapeptides were only found in Type I nsLTPs. Our result is consistent with previous studies4,33. Besides, the properties of the amino acid may determine the Cys pairing style, thus influencing the overall folding of the protein30,32,34. Generally, nsLTP1 has a hydrophilic amino acid in the CXC motif, whereas nsLTP2 contains a hydrophobic residue in the same position30,32. In our study, Motif 10 represented the pattern with a hydrophilic residue separating Cys5 and Cys6 in Type I, while Motif 1, 2, 6 and 15 corresponded to the Cys5XCys6 with a hydrophobic residue in the central position in the other subfamilies. Among the five hydrophobic residues in the CXC motif, Leu is the most frequent residue (75.40%), and this result is in accordance with previous studies3,4. Besides, Cys5-Val-Cys6, which corresponded to Motif 15, was only found in four Type VI nsLTPs.

The intron-exon pattern carries the imprint of the evolution of a gene family35,36. Previous studies have indicated the generality that some nsLTPs in Type I, III, IV, V and VI contained introns, while no introns were identified in Type II, VIII and IX nsLTPs3,4,32. However, our result showed some differences. In Gossypium, introns were found in Type VIII and IX, while no introns were identified in Type III (Fig. 3). The nsLTP family evolved in early diverged land plants, and new subfamilies evolved during land plant evolution2,32. Intron loss events were considered the main cause for the formation of novel nsLTP types, and also contributed to novel gene formation within the specific gene subfamilies32. Therefore, compared with the nsLTPs in other plants, Type IX and some Type VIII nsLTPs in Gossypium still possessed introns, while Type III genes have evolved with no introns harbored. The adoption of different nsLTP types could help evolved plants adapt to the stressful conditions on land gradually, so the nsLTP family was selected and expanded in land plants2,32.

The role of gene duplication in the genesis of evolutionary novelty and complexity has long been recognized37,38. The results in our study indicated that both tandem and segmental duplication events contributed to the expansion of the nsLTP family in G. arboreum and G. hirsutum, while tandem duplication contributed to the nsLTP expansion in G. raimondii (Fig. 4, Table 2). Previous studies have indicated that the ancient and recent whole-genome duplication (WGD) events occurred in Gossypium approximately 115–146 and 13–20 MYA, respectively, and the predicted hybridization and formation of allotetraploid cotton occurred approximately 1.5 MYA14,15,23. In the current study, most of the nsLTP duplication events in G. arboreum and G. raimondii might have occurred fewer than 32.69 MYA, and were considered the recent duplication events. The duplication event occurred 194.74 MYA in G. arboreum was therefore regarded as the ancient hexaploidization event. In G. hirsutum, the duplication events occurred approximately 5.84–24.50 MYA, suggesting that these duplication events might have occurred before the interspecific hybridization and subsequent polyploidization in G. hirsutum.

In addition, the orthologous genes in G. arboreum, G. raimondii and G. hirsutum were identified (Fig. 5, Table S4). The results revealed that 31 GaLTPs lost their orthologs in the At subgenome of G. hirsutum and 24 GrLTPs lost their orthologs in the Dt subgenome. These genes were mainly distributed in Type I, Type II, Type IV and Type V, and the proportions of these GrLTP and GaLTP genes distributed in each subfamily were complementary (Fig. S3). For example, in Type I, 66.67% of the GaLTPs had no orthologs in the At subgenome, whereas a relatively lower ratio (50%) of the GrLTPs lost their orthologs in the Dt subgenome. On the contrary, in Type V, only 33.33% of the GaLTPs lost their orthologs, while many more GrLTPs (75%) lost their orthologs. The ratios in other subfamilies such as Type IV, VI and VIII showed similar complementary phenomena, suggesting that the loss and retention of nsLTPs in G. hirsutum might not be random. During the process of polyploidization, chromosomal recombinations and rearrangements resulted in gene losses and gene retentions as duplicate orthologs39,40. The regular losses in conjunction with retentions might have avoided severe functional loss and redundancy of the genes in plants.

Duplications of genomic content often occur during DNA replication and recombination through independent mechanisms38,41. Some duplicated genes might have been subjected to diversification of their functions due to the divergent changes in their protein sequences together with cis-regulatory elements37. The nsLTPs may have undergone independent duplication events followed by functional diversification in each species during speciation42,43. In rice and wheat, the nsLTPs may have suffered from a complex evolutionary selection mechanism including subfunctionalization43,44. The functional diversification of nsLTPs appears to be not random, but seems to have originated from the specific biological processes over the course of evolution. In the current study, five pairs of duplicated genes shared similar expression profiles (Fig. 6). These results indicated that these duplicated genes might retain some essential functions during subsequent evolution, and the similar expression patterns may be related to the highly similar protein architecture of these duplicated genes. However, three paired duplicated genes showed significant expression divergence (Fig. 6). The diversity might be caused by the significant variation in gene regulation after the duplication events45,46. The differential expression patterns of duplicated genes in Gossypium demonstrated that these genes might have experienced functionalization during the evolutionary process.

As a prominent mode of speciation in flowering plants, allopolyploidy can cause genomic changes, including chromosomal rearrangements and changes in gene expression45,46,47,48. Recent studies have indicated the unequal expression of orthologous genes in allotetraploids such as Gossypium13,49, Brassica47, Oryza minuta50 and the synthetic allopolyploid Arabidopsis51, and the gene expression levels in these allopolyploids deviated from the expression of their orthologous genes in the progenitors. In our study, the expression patterns of genes in G. hirsutum revealed considerable differences from the expression of their orthologous genes in its progenitors (Fig. 7 and Fig. S5), indicating that the functions of the orthologous genes were strongly differentiated over the course of evolution.

In summary, 189 nsLTPs were identified in G. raimondii, G. arboreum and G. hirsutum in this study. Comprehensive study of the nsLTPs in Gossypium provided some important features of the gene family such as gene structure, expansion and expression. The findings here provide researchers with novel findings about the molecular evolution and expansion history of the nsLTP family in Gossypium, and offer a good opportunity to further investigate the nsLTP family in plants, especially polyploid plants. Besides, the analysis of the cotton nsLTPs in response to V. dahliae might help guide future identification of disease-resistant genes and exploitation of wilt-resistant varieties.

Materials and Methods

Plant materials and V. dahliae spore treatment

The seeds of G. raimondii, G. arboreum (Shixiya1) and G. hirsutum (TM-1) were kindly supplied by National Medium-term Gene Bank of Cotton in China. Cotton seeds were surface-sterilized with 0.1% HgCl2 for 10–15 min and washed with sterile water for at least 5 times. Sterilized seeds were sown on the half-strength Murashige and Skoog (MS) solid medium. After 3-day germination, the seedlings were transferred to half-strength Hoagland nutrient solution at pH 6.0. Fresh nutrient solution was continuously aerated with pumps and air stones, and was changed every 3 days. Plants were cultivated in a growth chamber at 28 °C with a photoperiod of 16 h light and 8 h dark. For tissue-specific expression analysis, roots, stems and leaves were harvested from 3-week-old cotton seedlings, and quick-frozen in liquid nitrogen, then stored at −80 °C until RNA isolation. Three plants constituted an individual replicate and three biological replicates were sampled for RNA extraction.

V. dahliae isolate Anyang (ACCC no. 36207) is a popular V. dahliae strain isolated from Anyang, Henan Province, China. It was supplied by the Cotton Research Institute of the Chinese Academy of Agricultural Science and preserved at the Agricultural Culture Collection, Beijing, China. V. dahliae mycelia growing on potato dextrose agar (PDA) plates were transferred to Czapek’s medium (containing 2.0 g·l−1 NaNO3, 1.0 g·l−1 K2HPO4, 0.5 g·l−1 MgSO4 •7H2O, 1.0 g·l−1 KCl, 0.01 g·l−1 FeSO4•7H2O and 30.0 g·l−1 sucrose) on a shaker at 150 rpm at 25 °C for sub-culturing. The spore concentration was determined with a haemocytometer and adjusted to a concentration of 107 spores ml−1 with sterile distilled water prior to use. Three-week old plants were inoculated with V. dahliae spores using the root-dip method, and the leaves were collected 0 h, 6 h, 12 h and 24 h after V. dahliae inoculation to examine the expression patterns of the nsLTP family in cotton under V. dahliae treatment.

RNA isolation and quantitative real-time PCR (qRT-PCR)

Total RNA was extracted using an RNAprep Pure Plant kit (Tiandz, China), and first-strand cDNA was synthesized from DNase-treated RNA with a PrimerScript 1st Strand cDNA synthesis kit (TaKaRa). Gene-specific primers were designed based on their coding sequences (CDSs) and then synthesized commercially (Generay, Shanghai, China) (Table S5). qRT-PCR was performed with SYBR premix Ex Taq (TaKaRa) and the CFX96 Realtime System (Bio-Rad, France) by strictly following the manufacturer’s instructions. The qRT-PCR machine was programmed with 40 cycles and an annealing temperature of 60 °C. The cotton EF1α gene was used as an endogenous control for all the qRT-PCR analyses. The relative transcription levels were calculated using the 2−ΔΔCT method. Three technical replicates were performed for each sample. The expression profiles of the nsLTPs were clustered using Cluster 3.0 software52.

Sequence retrieval and structural analysis

All the nsLTP sequences of Arabidopsis (Table S1) were collected from TAIR (http://www.arabidopsis.org/index.jsp) and were used as queries by searching against the cotton genome database14,15,23 using the BLASTP program with default parameters (E_value ≤10) (http://cgp.genomics.org.cn/page/species/index.jsp). The conserved LTP domain of each putative LTP sequence in Gossypium was confirmed in the Conserved Domain Database in National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov). Then, the candidate nsLTP sequences were manually examined for the presence of the 8 CM, and proteins lacking the essential Cys residues were excluded. After that, the putative proline-rich proteins and the proteins without NSSs (checked with SignalP 4.0, http://www.cbs.dtu.dk/services/SignalP) were also removed. Additionally, the proteins with C-terminal GPI anchor signals (predicted with big-PI Plant Predictor, PSORT and PredGPI) were excluded53,54,55. The protein sequences of the 2S- albumins (At2S1 to At2S4)56 and alpha amylase inhibitor (RATI)57 were then BLAST-searched against the rest of the candidate nsLTPs to exclude possible inhibitors and cereal storage proteins. Moreover, the proteins with more than 120 amino acids at maturity were discarded. Through these searching and screening procedures, the rest of the nsLTPs were finally confirmed and used for the following analysis. Primary and secondary protein structures were predicted with ProtParam (http://web.expasy.org/protparam/) and SOPMA (http://npsa-pbil.ibcp.fr/cgi-bin/npsa_automat.pl?page=npsa_sopma.html). In addition, the nsLTPs were clustered using OrthoMCL to identify the orthologous groups (OGs)58. The inflation parameter was set to the default.

Phylogenetic analysis and sequence alignment

All the nsLTP sequences in this study were aligned using ClustalX version 2.159. The PHYLIP program (http://evolution.genetics.washington.edu/phylip.html) was used to estimate the maximum-likelihood phylogeny for all the nsLTP sequences in Arabidopsis and the three Gossypium species with the JTT model. In addition, MrBayes version 3.1.2 was used to conduct Bayesian analysis for the nsLTP sequences in the three Gossypium species60. Trees were visualized with Figtree version 1.4.0 (http://tree.bio.ed.ac.uk/software/figtree/).

Protein motif and gene structure analysis

The conserved motifs of the nsLTPs were investigated using the online MEME program61. The analysis was performed with a set of parameters as follows: the optimum motif width was set to ≥6 and ≤200; the maximum number of motifs was set to 2062. Sequence logos for the conserved domains were generated using the WebLogo tool (http://weblogo.berkeley.edu/). Genomic schematic diagrams of the nsLTPs were obtained by comparing the genomic sequences and their predicted CDSs using the GSDS tool (http://gsds.cbi.pku.edu.cn/).

Chromosomal mapping and gene duplications

The chromosome location information of the nsLTPs was searched in the cotton genome database. MapInspect software (http://www.plantbreeding.wur.nl/uk/software-mapinspect.html) was used to generate chromosomal distribution images for these nsLTPs in G. raimondii, G. arboreum and G. hirsutum. Gene duplication events were investigated as described in previous reports62. Orthologous groups of the nsLTP family among the three Gossypium species were generated using the OrthoMCL database. Conserved syntenic blocks were inferred by running the OrthoClusterDB tool available at GDR (http://genome.sfu.ca/cgi-bin/orthoclusterdb/runortho.cgi). The diagram was plotted using Circos63. Furthermore, the selection pressure for each orthologous gene pair was calculated by the Ka (non-synonymous substitution rate)/Ks (synonymous substitution rate) ratio64. The formula “t = Ks/2r” was used to estimate the divergence time, and a neutral substitution rate (r) of 2.6 × 10−9 was used in the current study.

Additional Information

How to cite this article: Li, F. et al. Genomic Identification and Comparative Expansion Analysis of the Non-Specific Lipid Transfer Protein Gene Family in Gossypium. Sci. Rep. 6, 38948; doi: 10.1038/srep38948 (2016).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.