Identification and functional analysis of PIN family genes in Gossypium barbadense

Planting and sampling of cotton material

We used PimaS-7 (low fiber strength) and 5917 (high fiber strength) G. barbadense resource materials that had been preserved by the Key Laboratory of Crop Genetic Improvement and Germplasm Innovation, College of Agriculture, Xinjiang Agricultural University, Urumqi, China. The fiber characteristics of the samples are provided in Table 1.

The test materials were planted in an experimental field of the Xinjiang Academy of Agricultural Sciences, Urumqi, China. Cotton fiber samples were collected at 0, 5, 10, 15, 20, 25, 30, and 35 days post-anthesis (DPA), immersed in liquid nitrogen, and refrigerated at –80 °C until use.

Identification and physicochemical properties of PIN genes

We downloaded PIN genomic data for G. hirsutum TM-1 (ZJU-AD1_v2.1_a1.0) (Hu et al., 2019), G. barbadense Pima3–79 (HAU_v2.0) (Wang et al., 2019), Gossypium arboreum (CRI-A2_v1.0) (Du et al., 2018), and Gossypium raimondii (JGI_221_v2.1) (Paterson et al., 2012) from the Cottongen database (http://www.cottongen.org). PIN family accession numbers were cross-referenced with Pfam database (https://pfam.xfam.org/) and the associated Markov model files were downloaded. A local database was constructed using HMMER 3.3.2 software to retrieve the initial files, and the amino acid sequences were reconstructed and loaded into the Pfam database for searching (Ji et al., 2022). All amino acid sequences that did not contain the desired conserved structure were deleted from the original sequence, and a complete PIN family member result was obtained. The chromosomal locations and DNA lengths were verified using the cottonfgd database (https://cottonfgd.net/) and ExPASy program (https://web.expasy.org/protparam/) was used to predict the length, molecular weight (MW), isoelectric point (pI), instability coefficient, and hydropathic index of the corresponding proteins.

Phylogenetic analysis and gene chromosome mapping and replication

The amino acid sequences of PIN genes were loaded into the MEGA-X software, and the neighbor-joining method was used to construct an evolutionary tree, which was visualized using the iTOL online tool (https://itol.embl.de/). PIN gene conserved domain positions were determined using the National Centers for Biotechnology Information (NCBI) browser (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and the Clustal Omega tool (https://www.ebi.ac.uk/tools/msa/clustalo/); Jalview software was used for sequence alignment and to visualize the results.

The TBtools online toolkit was used to obtain chromosome density information from genome annotations, and position the PIN genes on chromosomes. Amino acid sequences of PIN genes in the four cultivated cotton species were aggregated and collinear scanning was performed using MCScanX (Wang et al., 2012).

Gene structure, motif, and promoter cis-element analysis

The amino acid sequences of PIN family members were loaded into the MEME tool (https://meme-suite.org/meme/tools/meme); the first 2,000 bp of the start codons of PIN genes were loaded into the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to predict the motif structure and promoter elements of the PIN genes.

Transcriptome analysis of PIN genes and determination of endogenous auxin content

The steps for the transcriptome analysis were as follows. Download transcriptome Hai7124 and TM-1 data at NCBI (https://www.ncbi.nlm.nih.gov/), then remove PimaS-7 and 5917 transcripts in the laboratory Group data (Hu et al., 2019). Using the G. hirsutum and G. barbadense genomes as reference genomes, conduct sequence alignment with HISAT2 (Kim, Langmead & Salzberg, 2015), then use feature Counts (Liao, Smyth & Shi, 2014) to obtain the raw count of each gene in each sample. Finally, import the data into R and use edgeR (Robinson, McCarthy & Smyth, 2010) for expression analysis. The method for determining endogenous auxin (IAA) in each period of fiber development adhered to the instruction manual of the plant auxin (IAA) enzyme-linked immunosorbent assay kit for fiber development, and this method was used to measure 0–35 DPA endogenous growth of PimaS-7 and 5917 (Fig. S1). SPSS was used to analyze the correlation between the endogenous auxin content and fragments per kilobase of exon model per million mapped fragments (FPKM) value of the PIN family gene at each stage of fiber development.

RNA extraction and quantitative reverse-transcription polymerase chain reaction (qRT-PCR) analysis

Total RNA was extracted from the collected samples using a polysaccharide polyphenol plant extraction kit provided by Beijing Tiagen Reagent (Beijing, China), according to the manufacturer’s instructions. The ABM reverse transcription kit was used to reverse convert the extracted RNA into first-strand cDNA, and the reaction was conducted according to the manufacturer’s instructions. qRT-PCR analysis was performed as follows. The primer GB-UBQ7 was designed using the cotton UBQ7 gene as the internal reference gene, and the relative expression levels of the target gene were calculated using the 2^–ΔΔCt method. The fluorescence quantitative reaction system consisted of 10 µL EVAGreen Express 2× qPCR, 2 µL cDNA template, 1 μL each of forward and reverse primers (10 μmol/L), and water to a final volume of 20 μL. The reaction conditions were 40 cycles of 94 °C for 30 s, 94 °C for 5 s, and 60 °C for 30 s.

Transient silencing analysis of the GbPIN7 gene

A cotton gene silencing system was established based on the tobacco rattle virus (TRV) vector (Liu & Page, 2008). The 320 bp specific nucleic acid fragment of the GbPIN7 gene was subcloned and loaded into the target site of a pTRV2 vector. Agrobacterium containing pTRV1, pTRV2-GbPIN7, and pTRV2-CLA was mass cultured, followed by 1:1 configuration of the transformed bacteria. Working solutions of the experimental group (pTRV1+pTRV2::GbPIN7), positive control (pTRV1+pTRV2::CLA), and negative control (pTRV1+pTRV2 empty) were each injected into PimaS-7 and 5917 plants that had not yet grown true leaves. After inoculation, the plants were grown at a constant temperature of 25 °C in an artificial culture room under a 16 h/8 h light/dark cycle. After 2 weeks, the albino phenotype was observed and photographed. Leaves were collected from the negative control and experimental groups for qRT-PCR analysis.

Evolutionary analysis of GbPIN family members with respect to physical/chemical properties and chromosomal localization

A total of 138 PIN family genes were identified among the four major cotton genera, of which 45 were identified in G. barbadense, 46 in G. hirsutum, 24 in G. raimondii, and 24 in G. arboreum. The resulting phylogenetic tree (Fig. 1) had three primary branches; the PIN and PILs genes were clearly distinguished. Analysis of the conserved domains showed that the 138 genes were divided into seven subgroups: PIN1, PIN2/PIN3/PIN4, PIN5/PIN6/PIN8, and PILS1/PILS3 (Fig. 2). The PILS5, PILS6, and PILS2 results further refined the subfamily positions of each gene; the phylogenetic branches and conservative structure analysis results matched completely, indicating the reliability of the classification results.

The PIN family genes of G. barbadense, G. hirsutum, G. arboretum, and G. raimondii were distributed on 20, 21, 10, and 11 different chromosomes, respectively (Fig. 3). The GbPILS9, GbPILS14, and GbPIN4 genes in G. barbadense were identified as the extrachromosomal Scaffold gene structure, similar to the G. arboreum GaPILS7 gene. These results may have been caused by our gene assembly and annotation methods.

A total of 23 genes were identified using descriptive statistics for 45 sea-island cotton genes; these 45 PIN family genes are widely distributed on 20 different chromosomes; they are absent from only six chromosomes (A04, A06, A08, D06, D08, and D13). The lengths of the DNA sequences varied from 767 to 13,644 bp, and those of the coding sequences varied from 576 to 1,938 bp.

The amino acid sequence lengths ranged from 191 to 645; the pI values ranged from 4.85 to 9.65; protein molecular weights were distributed from 20.97 to 70.48 kDa; the hydropathic index ranged from 0.21 to 0.76, indicating that all members were hydrophobins. Among the 45 PIN family genes, 16 were unstable proteins, and 29 were stable structures. All members of the three PIN subgroups were localized on the cell membrane, whereas all members of the four PILS subgroups were localized outside the cell (Table 2).

Gene structure and promoter analysis of PIN family genes in sea-island cotton

The gene structure analysis results showed that among the 45 PIN family genes, all members of the three PIN subfamilies except GbPIN13 and GbPIN14 contained six motif structures, three of which were at the beginning of the sequence (Fig. 4B). The GbPIN13 gene contained only motif2 and motif6, and GbPIN14 contained only motif3 and motif4. Fewer motif structures were found in each member sequence of the four PILS subgroups, and the most conserved motif1 structure was absent from 10 genes including GbPILS3 and GbPILS5. This indicates that sequence conservation of PIN subfamily genes was significantly higher than that of PILS subfamily genes. In terms of the sequence structure of gene family members, only GbPILS21, GbPILS22, and GbPILS23 had no introns, whereas most genes had large numbers of introns. PIN subfamily genes had significantly fewer introns. Among the four PILS subfamilies, nine genes, including GbPILS11, GbPILS8, and GbPILS14, did not have untranslated regions (UTRs). Two genes (GbPILS7 and GbPILS9) contained only 3′ UTRs, and the remaining 34 genes had two UTRs (Fig. 4C).

To analyze the transcriptional regulation and potential functions of the cotton PIN family genes further, we predicted transiently acting elements in the promoter region (Fig. 4D). The results showed abundant regulatory elements in the promoter region, mainly those for the phytohormones gibberellin, abscisic acid, and auxin; as well as responses to abiotic, anaerobic, low-temperature, drought, hypoxia, and light stresses; and biological growth processes including circadian rhythm, seeds, zein metabolism, meristem expression, endosperm expression, and flavonoid biosynthesis. These results indicate that members of the G. barbadense PIN gene family are involved in plant growth, differentiation, regulation, and responses to various biotic and abiotic stresses.

PIN family gene collinearity analysis

Comprehensive genome collinearity analysis results for the four cultivated cotton species (G. barbadense, G. hirsutum, G. raimondii, and G. arboreum) are shown in Fig. 5. PIN gene family members had 22 sets of collinearity in G. barbadense and G. hirsutum, three sets of collinearity in G. raimondii, and no collinearity in G. arboreum. These results indicate that each member of this gene family in G. arboreum remains highly independent in diploid cotton material; after combining to form tetraploid cotton, they corresponded to each other on chromosomes A and D, but then were lost, possibly due to duplication of function. There were 63, 42, and 22 sets of collinear relationships between G. barbadense and G. hirsutum, G. raimondii, and G. arboreum, respectively. These results may simply correspond to numbers of genes; however, the lack of collinearity between G. barbadense and G. arboreum may indicate that G. arboreum PIN family genes were largely lost through recombination into the G. barbadense gene.

Expression profiling and endogenous auxin correlation analysis of PIN family genes

We analyzed the FPKM values of the PIN family genes in transcriptome data from Hai7124, TM-1, PimaS-7, and 5917, and found that approximately 40% of PIN family genes were not detected in these four materials. Analysis of the Hai7124 and TM-1 transcriptomes revealed that 21 PIN genes had significant tissue-specific expression, and multiple genes were differentially expressed at different stages of fiber development. The expression levels of GBPIN1/GHPIN1, GBPIN2/GBPIN2, GBPILS1/GHPILS1 and 10 other genes were significantly higher in roots and stems than in leaves, and those of 11 genes including GBPIN11/GHPIN11, GBPIN12/GHPIN12, and GBPILS22/GHPILS22 were significantly higher in stems and leaves than in roots; these genes may be involved in regulating the local accumulation of auxin in plant organs. Six genes including GBPILS18/GHPILS18/GbPILS18 and GBPILS19/GHPILS19/GbPILS19 were expressed at 0–35 DPA, during fiber development, and the overall expression trend showed slightly higher expression at 0–35 DPA (Figs. 6A–6D).

Endogenous auxin content in cotton fibers was significantly lower in PimaS-7 than in 5917 at 10, 15, and 20 DPA, and that in 5917 was significantly lower at 15 and 20 DPA than at 10 and 30 DPA, respectively. In both materials, auxin content first increased, then decreased, and then increased again. In PimaS-7, auxin content peaked at 5 and 30 DPA, whereas that of 5917 peaked at 10 and 30 DPA; auxin content was significantly higher in 5917 than in PimaS-7 (Fig. 7A).

Correlation analysis showed that the expression levels of only five genes were significantly correlated with endogenous auxin content in PimaS-7 and 5917 at different stages of fiber development. GbPIN1 and GbPIN7 were negatively correlated with auxin content, whereas GbPIN16, GbPILS5, GbPILS14 were positively correlated (Fig. 7). However, as shown in Fig. 6, the expression levels of five genes (GbPIN1, GbPIN16, GbPILS5, and GbPILS14) were very low at various stages of fiber development; only GbPIN7 had a high expression level, which differed between the two materials and was directly related to endogenous auxin content.

qRT-PCR analysis of GbPIN family genes

Further analysis of the PIN I subgroup showed that 8 of 10 genes in this subgroup had detectable transcripts in the transcriptome data. Therefore, qRT-PCR analysis was performed on these eight genes; none of them showed consistent patterns in the fibrous tissues from 5 to 30 DPA, and some differences were observed between the two materials (Fig. 8). Three genes (GbPIN3, GbPIN4, and GbPIN8) showed a trend of increasing expression from 5 to 25 DPA, followed by low expression (Figs. 8A, 8D, 8H). The expression levels of four genes (GbPIN1, GbPIN2, GbPIN5, and GbPIN6) increased at 5–15 and 20 DPA, and decreased thereafter in 5917, whereas those of PimaS-7 gradually increased from 5 to 25 DPA, and then decreased (Figs. 8A, 8B, 8E, 8F). Thus, the expression of these four genes peaked later in the low-strength fiber material than in the high fiber-strength material, suggesting that these genes influence secondary wall thickening. The expression trend of the GbPIN7 gene was relatively consistent between the two materials, gradually increasing until 20 DPA, and then slowly decreasing (Fig. 8G).

Functional verification of the GbPin7 gene

Because GbPIN7 gene expression was observed in both the transcriptome and fluorescence quantitative analysis results and significant differences were detected between the two research materials, we conducted further transient silencing analysis of the GbPIN7 gene. The leaves of positive plants change from green to white after 15 days (Figs. 9A and 9B). Silencing efficiency was confirmed in the experimental and negative control groups using qPCR (Fig. 9C). Following silencing, the plants were cultured for 45 days under the same water and fertilizer conditions. Plant growth was significantly accelerated in both PimaS-7 and 5917 materials (Fig. 9D). Endogenous auxin content was measured in the upper leaves and stems of the same plants; the results showed significantly reduced auxin content.