Gossypium hirsutum (Xuzhou 142), G. arboreum (Shixiya 1), and G. raimondii (CMD10) that were used in this study were grown in a climate-controlled greenhouse with a 16 h light and 8 h dark cycle at 30°C, as previously reported (Shi et al., 2006). For tissue-specific expression patterns, the roots, stems, leaves, and flowers that were collected from the same growth stage of G. hirsutum and G. arboreum were immediately frozen and stored in liquid nitrogen after harvest.

The cotton genome sequences were downloaded from the CottonGen website¹. The Arabidopsis genome sequence was acquired from TAIR 10². Other plant genome sequences that were used in this study were retrieved from Phytozome³. We used the HMMER software with default parameters and the conserved PIN transmembrane domain (PF03547.13) to search for corresponding protein sequences. The PIN sequences were further identified based on a homology search using the BLAST program. Multiple sequence alignments of all identified PINs in this study were performed using ClustalX (Thompson et al., 1997). A phylogenetic tree of deduced amino acid sequences was constructed using the neighbor-joining algorithm with default parameters and 1000 bootstrap replicates in MEGA 5.0⁴. The transmembrane topologies of the PINs were predicted using the TMHMM web server with corresponding parameters⁵.

For the exon–intron structural analysis of the PIN genes, the coding sequences were compared with their genomic DNA sequences and the structure diagrams were generated using the GSDS server⁶. The conserved domains were predicted using the Simple Modular Architecture Research Tool (SMART⁷) and Pfam⁸ (Hu et al., 2015). The chromosome positions of all PINs in G. hirsutum were retrieved from gene annotation files that were downloaded from the CottonGen website. The map of the PIN gene distribution along the chromosomes was drawn from the top to the bottom using the CIRCOS program (Zou et al., 2013). We used MCScanX software for the detection and evolutionary analysis of PIN gene synteny and collinearity (Wang Y. et al., 2012).

The plant materials that were harvested as described above were frozen in liquid nitrogen and then ground to a fine powder with a mortar and pestle using a modified method (Jin et al., 2013; Liu et al., 2015). The total RNA was extracted using the PureLink^TM RNA mini kit (Invitrogen, Lot no.1687455) according to the manufacturer’s instructions, and the cDNA was reverse-transcribed from 5 μg of the total RNA as previously reported (Ji et al., 2003). In the qRT-PCR experiments, each gene was run in three biological and three technical replicates with the following reaction parameters: 95°C for 10 min, followed by 40 cycles of 95°C for 10 s and 56°C for 30 s. A melting curve was generated from 65 to 95°C. SigmaStat software was used for one-way statistical variance analysis. The UBQ7 cotton gene was employed as the internal control for all related results.

The in vitro ovule culture was performed as previously reported (Xiao et al., 2016). N-1-naphthylphtha-lamic acid (NPA, Sigma) was added to the culture medium with different concentrations as in the previous experimental design. Fiber lengths from cultured ovules were manually measured under a bright field microscope.

The full length (from the start codon to the stop codon) of each of the GhPIN1a_Dt, GhPIN6_At, and GhPIN8_At genes was amplified from the cDNA using the following gene-specific primers, respectively: 5′ primer: ACGGGGGACTCTAGAGGAT CCATGATCACTTTAACAGATTTTTACC 3′ primer: TCATA TACCCAATAAAATGTAGTACAGCTCTGTCGACTGTGGTA CCTTA 5′ primer: ACGGGGGACTCTAGAGGATCCATGATA ACAGGGGGTGATTTTTACA 3′ primer: GTATTACATACTTT TAGGCATATGAAGCTCTGTCGACTGTGGTACCTTA 5′ primer: ACGGGGGACTCTAGAGGATCCATGATTTCCCTGG CAGATGTTTATC 3′ primer: TCACAATGCTAAAAGAAGATA GTAGAGCTCTGTCGACTGTGGTACCTTA.

The pQG110 vector was digested using BamHI and KpnI. The product of each gene was then cloned into the digested pQG110 vector driven by the constitutive cauliflower mosaic virus 35S promoter via recombinase vazyme according to the manual, and the construct was further transformed into wild-type Arabidopsis plants. Transgenic plants were generated via the floral dip method and were selected on solid half-strength MS medium plates containing 50 mg/mL of the appropriate antibiotics (Zhang Y.Z. et al., 2015). The trichomes in the transgenic plant leaves were examined under a dissecting microscope (Leica MZ APO).

Southern blots were performed as described previously (Jin et al., 2013). 15 μg genomic DNA was digested thoroughly using BamHI or EcoRI. After purification, the digested DNA was loaded onto electrophoretic gels for 4 h, then blotted prior to hybridization and detected using the DIG High Prime DNA Labeling and Detection Starter Kit II (Roche) with GhPIN gene sequences as probe sequences.

Two diploid cottons, which were designated as AA genome and DD genome, diverged from a common ancestor approximately 2 to 13 million years ago (MYA) and developed genomes that differ approximately twofold in size. Approximately 1 to 2 MYA, these two genomes reunited in a common nucleus through allopolyploidization (Supplementary Figure S1), which led to the progenitor of the modern polyploidy cotton, G. hirsutum. This is the world’s primary commercial cotton species. Thus, allotetraploid cotton G. hirsutum contains two ancestral genomes, named the At and Dt subgenome. To identify all of the PIN proteins in G. hirsutum (AADD genome) and its two diploid ancestors, G. arboreum (AA genome) and G. raimondii (DD genome), we used the Arabidopsis PIN protein sequences to query the three reference genomes to screen out candidate PIN-like proteins in cotton. The collected PIN-like candidates were subjected to a further selection process based on their conserved domain using SMART. The PINs without the auxin efflux carrier domain were discarded. After this strict two-steps selection process, 17 deduced PINs were identified in G. hirsutum, along with 12 in G. arboreum and 10 in G. raimondii. Among the 17 PIN proteins identified in the G. hirsutum genome, seven originated from the At subgenome, nine from the Dt subgenome, and one was located on the scaffold (Supplementary Table S1). The lengths of the identified GhPIN proteins ranged from 127 (GhPIN8a_At) to 678 (GhPIN2_Dt) amino acids (aa) with an average length of 511 aa. The GaPIN proteins ranged from 124 (GaPIN3b) to 681 (GaPIN3) aa with an average length of 495 aa, and the GrPIN proteins ranged from 355 (GrPIN5) to 648 (GrPIN3) aa with an average length of 539 aa. We also analyzed orthologous PIN gene pairs in G. hirsutum and their corresponding diploid ancestors. Five genes from the A genome and two from the D genome were lost in G. hirsutum, whereas the orthologous genes of GhPIN8a_At and GhPIN8a_Dt in G. hirsutum were lost in its diploids, G. arboreum and G. raimondii, respectively.

To investigate the evolutionary relationships of PINs among different plant species, nine representative genome-sequenced dicotyledons that exhibited close evolutionary relationships to the three cotton species were used to construct a phylogenetic tree (Figure 1A). The evolutionary analysis showed that the PIN family genes in cotton are divided into long PINs (PIN1, PIN2, PIN3/4/7), short PINs (PIN5, PIN6, PIN8) and PIN9 (Figure 1B). According to the phylogenetic tree, the PIN family proteins can be divided into seven subgroups (Figure 1C). The PIN1 subgroup was most prominent in cotton. In contrast to most of the plant species analyzed, the PIN1 subgroup was extensively expanded in cotton, Populus trichocarpa and Glycine max. The PIN8 subfamily, contributing to auxin homeosatstic and facilitating the male gametophyte development in Arabidopsis (Ding et al., 2012), was extensively expanded in G. hirsutum, which may correspond to its better agronomic. PIN5 subgroup gene, which act antagonistically with PIN8, were missing in allotetraploid cotton, suggesting the loss of the gene occurred during the evolution of the G. hirsutum after polyploidization.

Compared with Arabidopsis, eight other plant species contained an additional PIN9 subgroup, which suggests a loss of this PIN in Arabidopsis during the evolutional course. Both PIN4 and PIN7 subgroups were only found in Arabidopsis, but none were found in the other species studied, indicating that these genes may have been acquired in Arabidopsis after it diverged from the common ancestor of the nine plant species used for analysis in the Figure 1A.

With the A-genome G. arboreum, D-genome G. raimondii, and the allotetraploid genome sequenced, we further explored the evolutionary history of the cotton PIN genes involved in gene duplication and loss. Our results showed that 12 GhPINs were evenly distributed on eight chromosomes, three from the At subgenome, and five from the Dt subgenome (Figure 2A). According to a whole genome analysis of gene duplications, we found that the whole genome duplication (WGD) and dispersed duplication, rather than the proximal duplication and tandem duplication, were the main driving forces of the GhPIN gene expansion (Supplementary Table S2). Specifically, nine GhPINs were produced via WGD, and three genes were produced via dispersed duplication among the analyzed genes.

Then, we analyzed the PIN gene conservation and loss based on the best match and the syntenic gene finder in MCScan. Three PINs were ultra-conserved homologous gene pairs in G. hirsutum. Of the orthologous gene pairs, six were lost from either the At or the Dt subgenome of G. hirsutum and two were lost from either G. arboreum or G. raimondii (Figure 2A). This indicated that the GhPIN genes experienced a higher frequency of genic sequence losses and only retained either the At or Dt homologous genes during the formation of G. hirsutum. Compared with one gene lost in the A genome and one lost in the D genome (Figure 2B), the rate of gene losses was higher in the allotetraploid cotton than in both diploid species. Additionally, one gene was only present in the Dt subgenome (Figure 2B).

To better understand the similarity and diversity of the gene structures and transmembrane topology of GhPINs, we constructed a phylogenetic tree with the deduced amino acids of GhPINs (Figure 3A) and investigated the exon/intron structures of GhPIN genes. In general, GhPINs possessed at least one intron and the structural patterns were highly conserved between long PINs and short PINs (Figure 3B). Although GhPIN6_At and GhPIN6_Dt genes belong to short PINs, both had five introns, including the two largest introns, in their open reading frame. The predicted transmembrane topology showed that GhPIN proteins, except for GhPIN1d_Dt, had the typical structure with two highly conserved hydrophobicity segments at the N- and C- terminus and a central hydrophilic loop within the two termini (Figure 3C). All GhPIN proteins possessed 8–10 transmembrane segments, except GhPIN8a_Dt, which presented six transmembrane helices. Notably, GhPIN8a_At possessed identical transmembrane segments as the other GhPINs, but these segments were evenly distributed across the entire protein, which resulted in the loss of the central hydrophilic loop.

Systematic analyses of gene expression patterns can help us understand the possible physiological functions of the cotton genes. First, we investigated the expression patterns of the GhPIN genes in different tissues, including the root, stem, leaf, flower, and fibers at various developmental stages. As shown in Figure 4, GhPIN transcripts were widely detected and exhibit a tissue-specific expression pattern in cotton. Genes from the PIN1, PIN2, PIN3, and PIN9 subgroups were predominantly expressed in vegetative tissues, but they were barely detected in the fibers, indicating that these genes may play specific roles in regulating the development of vegetative tissues. The transcripts from the PIN6 and PIN8 subgroup genes were abundant in diversified tissues. PIN1c_At from the PIN1 subgroup showed a root-specific expression, whereas the other four genes from the PIN1 subgroup were notably more highly expressed in the leaf (Figure 4A). Two PIN3 subgroup genes, PIN3_Dt and PIN3_At, were expressed in a leaf-preferential manner (Figure 4B). PIN6_Dt and PIN6_At, two genes from the PIN6 subgroup, were highly expressed in the 0 day post-anthesis (DPA) ovule (Figure 4C). Four PIN8 subgroup genes were significantly expressed in the stem (Figure 4D). The transcripts of PIN2_At genes were significantly higher in the root, whereas PIN2_Dt was specifically expressed in the flower (Figure 4E). PIN9_At was specifically expressed in the leaf.

Fiber length is one of the key agronomic traits that determine cotton yield and quality. First, we analyzed the fiber length in G. hirsutum and its diploid ancestors and found pronounced differences in G. hirsutum, G. arboreum, and G. raimondii (Figure 5A). The fiber length of G. hirsutum was more than 3 cm, whereas G. arboreum and G. raimondii showed a substantial reduction in fiber length (Figure 5B). To reveal the potential clues regarding whether the differential expression of PIN genes, which were generated during the origin and evolution of cotton, influenced the fiber development, we performed qRT-PCR to analyze the expression level of PIN genes during the different growth periods of the fibers from G. hirsutum and G. arboreum. Among the seven PIN genes analyzed, the PIN6 and PIN8 genes that originated from the At subgenome of G. hirsutum, rather than their homologs from the A genome of G. arboreum, were preferentially expressed in the fiber initiation and elongation stages (Figure 5C and Supplementary Figure S2). A comparative analysis of promoter sequences provides insights into the possible regulation mechanisms for the differential expression of PIN6 and PIN8 between G. hirsutum and G. arboreum. Compared with the GaPIN6 promoter, a large DNA fragment deletion was observed in the GhPIN6_At promoter. Additionally, a point mutation was identified in the GaPIN6 promoter, which resulted in the loss of an auxin response element (Figure 5D and Supplementary Figure S3). Although the GhPIN8_At promoter displayed a 98.1% similarity with the GaPIN8 promoter, we observed the loss of a pyrimidine box in the GaPIN8 promoter (Figure 5D and Supplementary Figure S4), which was reported to play an important role in regulating the expression of Ces3–4 gene from the At subgenome (Li F. et al., 2015).

We further compared the transcripts of the other 10 GhPIN genes with their corresponding orthologs in G. raimondii. In the qRT-PCR analysis, nine PIN genes showed no substantial differences in their transcriptional levels between the two cotton species, and only the expression level of PIN1a was approximately fivefold higher than its ortholog in G. raimondii (Figure 5E). Further promoter scrutiny revealed that several short sequences were missing in the GrPIN1a promoter (indels of 8–16 bp) compared with the GhPIN1a_Dt promoter, which resulted in the loss of a putative HD-ZIP binding site and a putative ZF-HD binding site (Figure 5F). Our results imply that the cis-element variations generated in these gene promoters during evolution may be coupled with the high-level expressions of the GhPIN6_At, GhPIN8_At, and GhPIN1a_Dt genes in allotetraploid cotton.

Increased expressions of PIN genes in longer fiber cells suggest that the auxin transport may promote fiber elongation. Therefore, we needed to confirm whether PINs play important roles in cotton fiber development. We utilized an auxin polar transport inhibitor to block the functions of the PINs and measured the fiber lengths with an ovule culture system in vitro. When NPA, an inhibitor of auxin polar transport, was added to the standard ovule culture media for 13 days, we observed a significant reduction in fiber elongation (Supplementary Figure S5A). Furthermore, the fiber lengths were reduced in a dose-dependent manner and the addition of 80 μM NPA to the culture media almost entirely blocked fiber growth (Supplementary Figure S5B). However, the application of different concentrations of NPA did not affect the final size of the cultured ovules (Supplementary Figure S5C). qRT-PCR analysis showed that the transcription levels of most PIN genes decreased after the application of NPA (Supplementary Figure S5D). These results suggest that the PIN proteins are essential for cotton fiber cell elongation.

To further investigate the functions of the GhPIN1a_Dt, GhPIN6_At, and GhPIN8_At genes identified through their differential expression patterns of the three cotton species, we generated three transgenic lines with GhPIN1a_Dt, GhPIN6_At or GhPIN8_At overexpressed in Arabidopsis that was driven by the constitutive 35S promoter (Figure 6A), and the expression level of GhPIN1a_Dt, GhPIN6_At and GhPIN8_At was 123, 48, and 87% of the Arabidopsis UBQ5 gene, respectively (Supplementary Figure S6). Meanwhile, the three transgenic lines used for phenotypic analysis were with single copy verified by southern blotting experiment (Figure 6B). In contrast to the transgenic lines carrying the empty vector, the three transgenic plants that carried the 35S::GhPIN1a_Dt, 35S::GhPIN6_At or 35S::GhPIN8_At constructs showed a considerable increase of number and length of trichomes, an organ similar to cotton fiber, in the leaves (Figures 6C,D). Furthermore, the calculated density of trichomes per leaf in the 35S::GhPIN1a_Dt, 35S::GhPIN6_At or 35S::GhPIN8_At transgenic lines increased to 250, 220, or 210% of the level in the leaf of transgenic line carrying the empty vector, respectively (Figure 6E). Compared with the transgenic line carrying the empty vector, the 35S::GhPIN1a_Dt, 35S::GhPIN6_At or 35S::GhPIN8_At displayed longer trichomes in the leaves (Figure 6F). These results suggest that the GhPIN1a_Dt, GhPIN6_At and GhPIN8_At genes possess biological activity for cell elongation in vivo, which strongly suggests that the high expression of the three PINs in G. hirsutum contributes to cotton fiber development.