Genome-wide identification of ABA receptor PYL family and expression analysis of PYLs in response to ABA and osmotic stress in Gossypium

Abscisic acid (ABA) receptor pyrabactin resistance1/PYR1-like/regulatory components of ABA receptor (PYR1/PYL/RCAR) (named PYLs for simplicity) are core regulators of ABA signaling, and have been well studied in Arabidopsis and rice. However, knowledge is limited about the PYL family regarding genome organization, gene structure, phylogenesis, gene expression and protein interaction with downstream targets in Gossypium. A comprehensive analysis of the Gossypium PYL family was carried out, and 21, 20, 40 and 39 PYL genes were identified in the genomes from the diploid progenitor G. arboretum, G. raimondii and the tetraploid G. hirsutum and G. barbadense, respectively. Characterization of the physical properties, chromosomal locations, structures and phylogeny of these family members revealed that Gossypium PYLs were quite conservative among the surveyed cotton species. Segmental duplication might be the main force promoting the expansion of PYLs, and the majority of the PYLs underwent evolution under purifying selection in Gossypium. Additionally, the expression profiles of GhPYL genes were specific in tissues. Transcriptions of many GhPYL genes were inhibited by ABA treatments and induced by osmotic stress. A number of GhPYLs can interact with GhABI1A or GhABID in the presence and/or absence of ABA by the yeast-two hybrid method in cotton.

In recent years, many PYL gene families have been characterized at genome-wide levels in rice, grape, soybean and other plants (Boneh et al., 2012;Kim et al., 2012;Bai et al., 2013;González-Guzmán et al., 2014;Tian et al., 2015;Fan et al., 2016;Gordon et al., 2016;Yu et al., 2016;Guo et al., 2017;Zhang et al., 2017). Furthermore, the structure properties of AtPYLs in Arabidopsis have been intensively investigated. It is found that an ABA-binding pocket, three α-helix, seven β sheets and four loops are conserved in AtPYLs Park et al., 2009;Santiago et al., 2012). Similar structures also exist in other plants. Furthermore, the expression patterns of PYL genes have been studied in tissues and in responding to exogenous ABA and diverse stresses in multiple plants (Saavedra et al., 2010;Boneh et al., 2012;Kim et al., 2012;Bai et al., 2013;González-Guzmán et al., 2014;Tian et al., 2015;Fan et al., 2016;Gordon et al., 2016;Yu et al., 2016;Chen et al., 2017;Zhang et al., 2017). These findings provide valuable information for researchers to further examine the functions of PYLs in ABA signaling in plants. However, knowledge about PYLs in Gossypium is scarce.
Cotton is the most important fiber crop and a key cash crop in the world. It provides the spinnable lint for the textile industry. It has been demonstrated that growth and development of cotton plants are seriously affected by a variety of environmental stresses such as drought, salinity and cold (Allen, 2010). These adverse stresses lead to significant decreases in the yield and quality of cotton fiber worldwide. Accordingly, it is of great importance to enhance the stress tolerance of cotton plants. One of the key strategies may be achieved via genetically engineering of PYL genes (Kim et al., 2014;Lee et al., 2015;Yu et al., 2016;Zhao et al., 2016;Chen et al., 2017;Liang et al., 2017). It is therefore essential for us to elucidate the functions and regulatory mechanisms of Gossypium PYLs. Here, we performed genome-wide and comprehensive analyses of the PYL family in G. arboretum (A 2 ) and G. raimondii (D 5 ), and their tetraploid species G. hirsutum (AD 1 ) and G. barbadense (AD 2 ). The expression patterns of PYL genes were studied in various tissues and in response to ABA and osmotic stress in G. hirsutum. Moreover, the interactions between individual GhPYLs and GhABI1A or GhABI1D were measured in G. hirsutum by the yeast-two hybrid method. These results may pave the way for further investigating the functions of cotton PYLs in the future.

Synteny and Ka/Ks analysis
The values of nucleotide substitution parameter Ka (non-synonymous) and Ks (synonymous) were counted based on the PAML program (Yang, 1997). The homologous genes were searched by the MCScanx software (Wang et al., 2012), and gene collinearity results were obtained by the CIRCOS program (Krzywinski et al., 2009). The syntenic maps of the PYL genes from G. arboretum, G. raimondii and G. hirsutum were constructed using the circos-0.69 ± 3 software package with default parameters (http://www.circos.ca).

Expression analysis of PYL genes in tissues and in response to ABA or osmotic stress
To monitor expression levels of PYLs in tissues, samples of roots, stems and leaves were obtained from TM-1 cotton plants grown in soil for 21 d. Flowers were picked in the morning at the first day of anthesis, and fibers at the secondary cell wall stage (about 23 d post anthesis) were sampled from the ovules. To determine PYL transcript abundances in responding to ABA or osmotic stress, cotton plants grew in liquid 1/2 MS medium (Murashige & Skoog, 1962) in a growth chamber (day/night temperature cycle of 28 • C/26 • C, 14 h light/10 h dark, and about 50% relative humidity) for 3 weeks. Then, the plant leaves were sprayed with 100 µM ABA or the roots were treated with 10% PEG6000 (PEG6000 was added in the liquid medium) for 3, 6, 12 and 24 h, respectively. The plants that were not sprayed with ABA or were not treated with PEG6000 were as controls. The roots for treatments and controls were collected, frozen in liquid nitrogen and stored at −70 • C. About 0.1 g of samples were ground to a fine powder with a mortar and pestle in liquid nitrogen. Total RNA was extracted from the powder using RNA Pure Plant Kit's protocol (TIANGEN Company). The cDNA was obtained by M-MLV reserve transcriptase synthesis system (Promega, Madison, WI, USA) following the instructions in the Promega kit (https://tools.thermofisher.com/content/sfs/manuals/superscriptIII_man.pdf).
Quantitative real-time RT-PCR (qRT-PCR) experiments were carried out applying the cDNA, SYBR Green Master mix, the specific primers of cotton PYL genes (Table S1, the primer sequences were submitted to NCBI), and an ABI 7500 real-time PCR system. Cotton UBQ7 was used as the internal control. Experiments were repeated at least three times.

Yeast two-hybrid assay
Full-length sequences of gene GhPYLs and GhABI1 were cloned into pGBKT7 and pGADT7 vectors, respectively, using primers listed in Table S2 (the primer sequences were submitted to NCBI). The resultant constructs were co-transformed into yeast strain AH109 (MATa,112,gal4 , gal80 , LYS2::GAL1UAS-GAL1TATAHIS3, GAL2UAS-GAL2TATA-ADE2, URA3::MEL1 UASMEL1TATA-lacZ, MEL1) according to the method described in page 18-21 in Yeast Protocols Handbook (Clontech Laboratories Inc., 2009). The cotransformants were plated on non-selective SD/-Leu/-Trp (synthetic dropout medium without Leu and Trp) agar plates and selective SD/-Leu/-Trp/-His/-Ade agar plates in the absence or presence of ABA. A cotransformed yeast spot was collected and diluted in sterile ddH 2 O. The optical density (OD) value for the first solution of yeast colony was set 1. Serial 1:10 dilutions were generated, and 2 µL of the dilution was dropped on an agar plate to obtain one spot. The interactions between GhPYLs and GhABI1 were observed after 4 d of incubation at 30 • C. The experiments were repeated at least three times.

Genome-wide analysis of PYLs in four cotton species
To investigate the PYL family in the cotton genomes, the 14 AtPYL gene coding sequences and amino acid sequences were applied as queries to search against the cotton genome databases. A total of 21, 20, 40 and 39 PYL genes were identified in the genomes of two progenitor diploid species G. arboretum and G. raimondii, and their derived tetraploid species G. hirsutum and G. barbadense, respectively. The PYLs in the two diploid species were named based on their orthologous similarity to the 14 AtPYLs according to the methods described by Mohanta et al. (2015) in other genes. Briefly, the first letter of the genus (upper case) and the first letter of the species (lower case) followed by PYL and an AtPYL number were used to name a Gossypium PYL. The number of a Gossypium PYL was the same to that of its orthologous AtPYL sharing the most similarity in protein sequences to the Gossypium PYL. When more than one Gossypium PYLs had the same ortholog in Arabidopsis, additional numbers followed by a hyphen were applied to distinguish among paralogs of the Gossypium PYLs. The small number after the hyphen means high similarity of a cotton PYL to its corresponding AtPYL. The PYLs of the tetraploid cotton plants were denominated based on their phylogenetic relationship with those in G. arboretum and G. raimondii; and the last letter A or D meaned that the PYL was derived from A or D genome (Table S3). Therefore, the Gossypium PYLs in the four species were named GaPYLs, GrPYLs, GhPYLs and GbPYLs, respectively. We noticed that the Gossypium PYLs named following AtPYL2, AtPYL4, AtPYL6, AtPYL9, AtPYL11, AtPYL12, and AtPYR1. Moreover, most of these AtPYLs possessed not only one ortholog in Gossypium. For example, 8, 7, 14 and 13 homologs of AtPYL9 were identified in genomes of G. arboretum, G. raimondii, G. hirsutum and G. barbadense, respectively (Table S3).
Analysis of the physical properties of the Gossypium PYL members revealed that these PYLs were highly conserved. Most PYLs had similar amino acid lengths, molecular weights (MWs), and theoretical isoelectric points (pI). Majority of PYLs in Gossypium possessed 177-222 amino acids. The MWs of the PYLs varied from 15.32 kDa to 32.22 kDa. The pIs of PYLs ranged from 4.73 to 9.51 with an average of 6.20. Most PYL proteins were predicted to locate in cytoplasm and/or nucleus (Table S3).

Phylogenetic and structural analysis of cotton PYLs
To explore the evolutionary relationship of the PYLs among G. arboreum, G. raimondii, G. hirsutum and G. barbadense, an unrooted phylogenetic tree for the 120 PYLs was generated ( Fig. 2A, Files S1-S2.). The PYLs can be divided into three subfamilies (I-III) based on the bootstrap values (>1,000) in the Neighbor-Joining (NJ) tree. Subfamily I had 41 members, which were orthologs of AtPYR1 (16) and AtPYL2 (25). Subfamily II consisted of 37 PYL genes. They were homologs of AtPYL4 (17), AtPYL6 (14), AtPYL11 (3) and AtPYL12 (3). Other 42 PYL members (homologues of AtPYL9) belonged to subfamily III ( Fig. 2A). We found that the Gossypium homologs of the same AtPYL frequently clustered closely, indicating their more closed relative relationship. Besides, several PYLs including GaPYL2-2, GrPYL2-2, GhPYL2-2A, GrPYL2-2D, GbPYL2-2A, GbPYL2-2D and GbPYL2-2D appeared to be distant clades from other PYLs in Gossypium ( Fig. 2A). This suggests that these PYL2-2 members had relatively distant phylogenetic relationship with other Gossypium PYLs. The exon/intron structures of the Gossypium PYL genes were studied. The results showed that the number of exons in the PYLs was 1-3. There was no intron in 66 PYL genes. Most of the remaining genes had two introns except that nine members had one intron (Fig. 2B). In addition, most Gossypium PYL genes clustered in the same subfamily had the similar number of exons and lengths of introns. For instance, vast majority of genes in subfamily I and II had only one exon, and 41 out of 42 genes in subfamily III had 3 exons and relatively long intron sequences (Fig. 2B). These results imply that the exon/intron organizations of cotton PYLs are closely related to the phylogenetic relationship of the genes.
To further clarify the diversity of motif compositions, the putative motifs in the 120 PYL proteins in subfamily I-III were analyzed by the MEME software. A total of 19 motifs designated as motif 1 to motif 19 were detected (Fig. 2C). The numbers of motifs in a PYL of subfamily I to III were 4-7, 4-7 and 3-6, respectively. Of the 19 motifs, motif 8 and 10 emerged in majority of PYLs in subfamily I, motif 3 and 5 were conservative among most PYLs in subfamily II; and motif 1, 2, 4 and 9 existed in overwhelming majority of subfamily III members. Most PYLs within the same subfamilies have very similar motif compositions and distributions, implying that PYLs within the same subfamilies probably share similar functions. Intriguingly, some motifs specifically existed in a particular subfamily or some PYLs. For example, motif 16 and 17 were only present in subfamily II, and motif 12 and 18 only belonged to members of subfamily III. Motif 11 was only seen in GbPYL2-3A and GbPYL9-4D. Motif 7 was only found in GhPYL2-4A and GbPYL6-2D ; and motif 15 was only observed in GbPYL6-2D . These findings suggest that these motifs might play specific roles in the PYLs of that subfamily or in those PYLs. The detailed mechanisms need to be experimentally examined in the future.

Conserved domains and amino acids of cotton PYLs
To better understand the structural similarity of Gossypium PYLs, the amino acid sequences of PYLs from G. arboreum, G. raimondii and G. hirsutum as well as those from the two diploid species and G. barbadense were aligned, respectively (Figs. S1 and S2). The results revealed that these PYLs had high sequence similarities. All the PYLs shared a similar helix-grip structure formed by seven β-sheets and three α-helices, and four identical conserved loops among the β-sheets and the α-helices (Figs. S1 and S2). Furthermore, many conserved amino acids (highlighted by red colour) were observed in the putative β-sheets, α-helices as well as loops (Figs. S1 and S2), consistent with the structure of AtPYLs in Arabidopsis Park et al., 2009). These conserved secondary structures and amino acids have been demonstrated to be essential for the functions of ABA receptors in Arabidopsis. For instance, the loops of CL1, CL2 and CL3 are essential for ABA bindings and the PYL-PP2C interactions (Zhang et al., 2015). Noteworthily, the conserved leucine (L) was replaced by methionine (M) in the CL2 loop, and the conserved arginine (R) was replaced by lysine (K) in the CL3 loop in GaPYL6-2, GbPYL6-2A and GbPYL6-2D (Figs. S1 and S2). Similarly, the conservative arginine (R) was replaced by methionine (M) in the CL3 loop in GhPYL9-5D (Fig. S1). These data hint that these PYLs may be significantly different in ABA bindings and interactions with PP2C from other Gossypium PYLs.

Chromosomal distributions of cotton PYLs
The localizations of the Gossypium PYLs in chromosomes were determined. Generally, PYL genes were unevenly distributed on multiple chromosomes (Fig. 3). The 21 GaPYLs, 20 GrPYLs, 40 GhPYLs and 39 GbPYLs separately placed on 10, 11, 21 and 18 chromosomes, respectively. Four genes were located on each of the Achr10 and Atchr05 chromosomes. Many chromosomes, each possessed three genes. These chromosomes included Achr06 from G. arboretum, Dchr07 and Dchr08 from G. raimondii, Atchr11, Atchr12, Dtchr05, Dtchr11 and Dtchr12 from G. hirsutum; and At chr01, At chr10, Dt chr08, Dt chr10, Dt chr11 and Dt chr12 from G. barbadense (Fig. 3). A large number of chromosomes individually owned two genes. These chromosomes were Achr03, Achr07, Achr09 and Achr13 in G. arboretum, Dchr04, Dchr06 and Dchr11 in G. raimondii, Atchr08, Atchr09, Atchr10, Dtchr08, Dtchr09 and Dtchr10 in G. hirsutum; and At chr12, Dt chr05 and Dt chr09 in G. barbadense. Each of the rest chromosomes had one gene. The distributions of the PYL genes on individual chromosome were irregular. Some genes were located on the upper end of the chromosome arms, some placed on the lower end of the arms; whereas some lied in the region far from two end of the arms (Fig. 3). In addition, two genes in G. arboretum (GaPYL2-1, GaPYL2-3) and 6 genes in G. barbadense (GbPYL4-2A, GbPYL9-3A, GbPYL9-5A, GbPYL9-6D , GbPYL11A and GbPYL12D) were present in scaffolds.
We compared the positions of the orthologs among GaPYLs, GrPYLs and GhPYLs in chromosomes. Unexpectedly, only a few PYL homologs localized in their corresponding homoeologous chromosomes. Similar results also happened among GaPYLs, GrPYLs and GbPYLs (Fig. 3). These findings hint that many complex conversion events of PYLcontained homoeologous chromosomes or of PYLs occurred among different Gossypium species during evolution.

Synteny analysis of PYL genes
It has addressed that gene duplication events including tandem and segmental duplications play key roles in expanding gene family during the evolutionary process (Cannon et al., 2004). To gain insight into the genetic origins and evolution of the Gossypium PYLs, we analyzed the homologous gene pairs of PYLs among G. arboreum, G. raimondii and G. hirsutum. A sum of 88 collinearity blocks were identified between G. Arboreum and G. raimondii, and each block had one gene pair. Ninety-one collinearity blocks with 93 homologous pairs were detected between the At-genome and Dt-genome of G. hirsutum. Moreover, two homologous gene pairs between chromosome At11 and Dt11 (GhPYL9-3A/GhPYL9-4D, GhPYR1-2A/GhPYR1-2D), and chromosome At12 and Dt12 (GhPYL9-5A/GhPYL9-5D, GhPYR1-3/GhPYR1-3D) were found in two individual collinearity blocks. Additionally, 395 homologous gene pairs distributed in 385 collinearity blocks among G. Arboreum, G. raimondii and G. hirsutum (Fig. 4, Table S4). Of these blocks, one harbored three homologous gene pairs (GrPYL9-1/GhPYL9-1D, GrPYL9-4/GhPYL9-4D, GrPYR1-2/GhPYR1-2D). The block was between chromosome D07 in G. raimondii and  (Table S4). No gene pair was implicated in tandem duplication. These results suggest that segmental duplications dominantly contribute to the generation of Gossypium PYLs during genetic evolution.

Analysis of Ka/Ks values of PYLs in G. arboreum, G. raimondii and G. hirsutum
To further investigate the divergence and selection in duplication of PYL genes, the nonsynonymous (K a), synonymous (K s) and K a/K s values were evaluated for the homologous  5). The average Ka/Ks value of all the 301 gene pairs was 0.17, less than 1. These data suggest that these genes were mainly under the purifying selection during evolution. In contrast, the ratio of Ka to Ks between GhPYL2-4A and GaPYL2-4 was 1.08, and that between GhPYL9-8A and GrPYL9-7 was 1.51, indicating that the two gene pairs may generate under positive selection. Besides, the Ka/Ks values of many gene pairs were 1 (Table S5), implying that these genes are under neutral evolution. , and a phylogenetic tree was constructed applying the neighbor-joining method and the MEGA 5.0 software. On the basis of the topologic structure, the PYL family members were classified into three subfamilies (I-III) (Fig. 6). Expectedly, vast majority of the Gossypium PYLs were closely clustered in a subfamily. Moreover, every subfamily contained PYLs from eudicots such as Gossypium, cocoa, castor and Arabidopsis, and from monocots like    B. distachyon and O. sativa, hinting that these PYLs generated before the divergence from eudicots and monocots. In addition, numerous PYLs from cotton species clustered more closely with those from cocoa than from other plants (Fig. 6), reflecting closer relationship of Gossypium PYLs with cacao PYLs relative to other plant PYLs.

Expression of GhPYL genes in tissues
To further understand the roles of cotton PYLs in diverse organs, the expression patterns of all the 40 GhPYL genes were monitored by qRT-PCR. Of these GhPYLs, 22 genes were preferentially expressed in the flower, 10 genes were dominantly expressed in the root. Moreover, three genes including GhPYR1-3D, GhPYL2-2A and GhPYL2-2D were highly expressed in the fiber. The transcripts of GhPYR1-3A, GhPYL4-3A, GhPYL4-3D, GhPYL9-1A and GhPYL12D were relatively abundant in the stem (Fig. 7). These results indicate that some ABA receptors may specially function in a unique tissue. Besides, more GhPYLs were predominantly expressed in the flower and the root, suggesting that ABA receptors may play important roles in fiber formation and in response to abiotic stresses in cotton.

Expression profiles of GhPYLs in responses to ABA and osmotic stress
The expression patterns of the GhPYL genes in roots were monitored after treatments with 100 µM ABA or 10% PEG6000 for different period of time. In general, the transcriptional levels of most GhPYLs were decreased in responding to ABA (Fig. 8), and moderately enhanced in response to osmotic stress (Fig. 9). After treatment with ABA for 3 h or 6 h, expression of many genes for example GhPYR1-1A, GhPYR1-3D, GhPYL9-5A significantly reduced, then increased slowly at 12 h or 24 h. Expression of some genes like GhPYL4-1D, GhPYL4-2A and GhPYL6-1D continually decreased in response to ABA. However, the transcriptional abundances of some genes increased upon ABA treatment. These genes included GhPYL2-2A, . By contrast, the expression of GhPYL9-2D was not influenced by ABA treatment (Fig. 8). These results indicate that GhPYLs play differential roles in perceiving ABA signals. The expression of most GhPYLs was upregulated by PEG treatments for 12 h and/or 24 h. By contrast, the transcriptional levels of GhPYR1-3, GhPYR1-3D, GhPYL2-1A and GhPYL2-4D were diminished whereas those of GhPYR1-2A and GhPYR1-2D were unchanged after exposure to high concentration of PEG6000 (Fig. 9). These results suggest that a great number of GhPYLs have different responses to ABA relative to osmotic stress in cotton.

Many GhPYLs interact with GhABI1A or GhABI1D
PYLs have been addressed to transduce ABA signals to downstream targets through selectively interplaying with clade A PP2C proteins such as ABI1 and ABI2 in plants (Cutler et al., 2010;Joshi-Saha, Valon & Leung, 2011). We therefore examined the interactions between GhPYLs and GhABI1A (Gh_A07G0123) or GhABI1D (Gh_D07G2383) in the absence or presence of ABA by yeast-two hybrid method. Twenty-five out of fourty GhPYLs were cloned because of high similarity of CDS sequences among different GhPYL members. In the absence of ABA, 9 GhPYLs individully interplayed with GhABI1A, and 8 GhPYLs respectively interacted with GhABI1D (Table 2). Among them, GhPYR1-1D, GhPYL6-2D, GhPYL9-1A and GhPYL9-5D displayed relative weak interactive signals with GhABI1D. GhPYL6-2A and GhPYL9-4D showed strong interaction signals with GhABI1D (Fig. 10). When supplied with 10 µM ABA in the medium, 17 GhYPLs could interact with GhABI1A, and 14 GhPYLs could interplay with GhABI1D. We observed that many interactions were ABA-dependent, and numerous interactions were ABA-independent (Fig. 10, Table 2). Moreover, the interactive intensities between GhPYL6-2A and GhABI1A as well as between GhPYL9-5D and GhABI1D were increased by ABA. Interestingly, the interaction signal between GhPYL4-2A and GhABI1A was slightly weakened by ABA. These results suggest that GhPYLs differentially or specifically bind to GhABI1 in response to stresses, reflecting the diverse interacting modes among GhPYLs, GhABI1 and ABA in cotton.  Fan et al., 2016;Yu et al., 2016;Zhang et al., 2017;Guo et al., 2017), indicating more complex ABA responses likely happen in Gossypium. Chen et al. (2017) identified 27 GhPYLs , significantly less than 40 members as described in this study. The main reason may be that the genome sequence database of G. hirsutum they chose to search homologs of AtPYLs was different from what we did. Chen et al. (2017) used the database from the website http://cgp.genomics.org.cn/page/species/index.jsp whereas we applied that from https://www.cottongen.org. Liang et al. (2017) reported that there exist 22, 22, 44 and 36 PYLs in G. arboretum, G. raimondii, G. hirsutum and G. barbadense, respectively, inconsistent with our results. The possible reason might be the standard we used to determine PYLs in Gossypium was slightly different from what they did. For instance, we also surveyed 44 putative GhPYLs. However, three of them (Gh_A10G0340, Gh_D10G0346 and Gh_D08G1226) lack the PYR_PYL_RCAR_like domain of PYLs, and Gh_Sca112916G01 encodes a too small protein. Accordingly, the four proteins were eliminated in this report.
Of the 40 GhPYLs we identified, 20 and 20 were from A and D subgenomes, respectively; and every member had its corresponding ortholog in a G. arboretum and a G. raimondii. Likewise, most GbPYLs individually possessed their corresponding orthologs in a G. arboretum and a G. raimondii. However, both GbPYL6-2A and GbPYL6-2D had high sequence similarity to GrPYL6-2. Similar cases occurred in GbPYL6-2A and GbPYL6-2D , GbPYL2-2D and GbPYL2-2D as well as in GbPYL9-6D and GbPYL9-6D . The three pairs of genes seemed to derived from GaPYL6-2, GrPYL2-2 and GrPYL9-6, respectively (Table S3). These findings suggest that gene replication events contribute to the generation of these GbPYL pairs in sea island cotton. Moreover, GaPYL9-4 had no corresponding homologous gene in GhPYLs, and some GaPYLs (GaPYR1-1, GaPYR1-2, GaPYL9-1,  and GrPYL4-2 did not share corresponding homologous genes in GbPYLs (Table S3), implying that gene loss events occurred during the evolution processes of Gossypium PYLs.
We observed that some AtPYLs (for example AtPYL9) had several homologues of GaPYLs, GrPYLs, GhPYLs or GbPYLs (Table S3), hinting that selective gene replications likely play essential roles in the generation of the Gossypium PYLs during evolution.
The PYLs of G. arboreum and G. raimondii exhibited similar physical properties to those of G. hirsutum and G. barbadense (Table S3). This suggests that the functions of the PYLs in the four cotton plants were conserved during the processes of genetic evolution. The predicted subcellular localizations of PYLs in Gossypium were cytoplasm and/or nucleus, in line with the locations of PYLs from Arabidopsis, rice, and soybean plants Bai et al., 2013;Tian et al., 2015). These data indicate that the functions of PYLs are conservative among different plant species.
We saw that each subfamily of Gossypium PYLs contained several members of G. arboreum, G. raimondii, G. hirsutum and G. barbadense. Moreover, the number of GhPYLs and GbPYLs was roughly sum of GaPYLs and GrPYLs. These findings imply that each diploid species already contained about 20 PYL members prior to formulating G. hirsutum and G. barbadense, and the GhPYLs and GbPYLs mainly came from the allopolyploidization (Wendel, Brubaker & Seelanan, 2010).
The organizations of intron/exon and the number of exons in the surveyed Gossypium PYLs were very similar to those in the orthologs in rice, maize, tomato, rubber tree and Brachypodium distachyon (González-Guzmán et al., 2014;He et al., 2014;Fan et al., 2016;Guo et al., 2017;Zhang et al., 2017). Moreover, all of PYLs in Gossypium had similar helixgrip structures consisting of seven β-sheets, three α-helices and four loops, consistent with PYLs in Arabidopsis, rice and other plant species González-Guzmán et al., 2014;He et al., 2014;Fan et al., 2016;Guo et al., 2017;Zhang et al., 2017). These results hint that these PYLs from different plant species have very similar functions and action modes.
Most Gossypium PYLs in subfamily I and II had no intron. However, some members in the two subfamilies such as GrPYR1-3, GhPYL2-4A, GrPYL4-2 and GbPYL4-1A contained one or two introns (Fig. 2), suggesting that these genes obtain new introns during evolution.
Colinearity results revealed that approximate 400 homologous gene pairs existed among PYLs from G. arboreum, G. raimondii and G. hirsutum (Fig. 4), and numerous PYL gene pairs like GbPYL2-2D and GbPYL2-2D had nearly identical amino acid sequences (Table S6). These results suggest that PYL family expands through segmental duplication events during evolution.
To uncover the homologous relationships of PYL gene family among different taxa, a phylogenetic tree was constructed based on the PYL protein sequences from four cotton species, T. cacao, R. communis, V. vinifera, B. distachyon, rice and Arabidopsis (Fig. 6). We observed Gossypium PYLs clustered more closely with cocoa PYLs than OsPYLs, in agreement with the evolutionary relationships among these plants.
Expression analysis results revealed that great majority of GhPYLs were expressed in various tissues like roots, stems, leaves, flowers and fibers. Notably, a large number of genes were highly expressed in the flower (Fig. 7). These results imply that ABA is implicated in modulation of reproductive development in cotton. Moreover, we noticed that some genes were preferentially expressed in roots, stems and fiber, pointing to the important roles of ABA in these tissues. PYLs have been documented to express in diverse tissues in various plants such as rice, maize and rubber tree (Tian et al., 2015;Fan et al., 2016;Guo et al., 2017). In rice, most OsPYLs were expressed in all tissues, OsPYL3 and OsPYL5 were predominantly expressed in leaves, and OsPYL1 in roots (Tian et al., 2015). In maize, PYL11 was upregulated in leaves, PYL6 and PYL10 in roots (Fan et al., 2016). In rubber tree, five genes were detected to express in all tissues tested, four genes were preferentially expressed in leaves, four in roots and one in flowers . These results were consistent with our findings in cotton, indicating the diverse biological functions of different PYLs in plants.
Changes in the expression levels of 40 GhPYLs were investigated in roots in responding to ABA or osmotic stress. The transcriptional abundances of many genes decreased after ABA treatments (Fig. 8), but increased post PEG treatments (Fig. 9), pointing to the different functional mechanisms of GhPYLs in response to ABA and osmotic stress. Tian et al. (2015) reported OsPYLs are differentially expressed after ABA treatment. Some AtPYLs have been addressed to be down-regulated by ABA treatment in Arabidopsis (Santiago et al., 2012). In B. distachyon, BdPYL11 is down-regulated under ABA and PEG6000 treatments while BdPP2CA4, BdPP2CA5, BdPP2CA6 and BdPP2CA8 are up-regulated (Zhang et al., 2017). ZmPYL4-11 was found to be down-regulated by ABA (Fan et al., 2016). Bai et al.