Genome-wide identification and expression reveal the involvement of the FCS-like zinc finger (FLZ) gene family in Gossypium hirsutum at low temperature

Plant material

We selected 86-1 (cold resistance) and Lumian 2 (cold susceptibility) plant cultivars as the materials for expression analysis under low-temperature stress (Ge et al., 2021). The stress temperature used was 12 °C. Seeds with uniform size were selected, and after disinfection with 75% alcohol, the seeds were soaked for 24 h until they turned white. Then, the seeds were sown into 10 cm × 10 cm pots filled with a mixture of vermiculite and sterilized farm soil (in a ratio of 1:2) at room temperature with a relative humidity of 70%. When the cotton plant was at the three-leaf stage, it was put into an artificial climate incubator (RGX-400P; Taisite, New York, NY, USA) at a temperature of 12 °C and a relative humidity of 70% with a 16 h light/8 h dark cycle, and the roots, stems and leaf tissues were collected at 0, 6, 12, 24 and 48 h under stress. Two materials were cryogenically treated with three biological replicates per time point. The samples were quickly frozen with liquid nitrogen and stored in a −80 °C freezer.

Identification and bioinformatics analysis of the FLZ gene family in cotton

The genomic and proteomic data of Gossypium arboreum, Gossypium raimondii, Gossypium hirsutum and Gossypium barbadense were obtained from the COTTONGEN (http://www.cottongen.org/) database. Protein sequences in cotton were identified by performing a hidden Markov model of the FLZ gene domain (PF04570) and a hidden Markov model (HMM) in Hmmsearch software (http://hmmer.org/) (Potter et al., 2018). Domain identification was performed on the GhFLZ protein sequence through the NCBI CDD database (Lu et al., 2020). Finally, ExPASy and PSORT were used to analyze protein sequence physical parameters and subcellular localization (Chang et al., 2013).

Phylogenetic and collinear analysis of the cotton FLZ gene family

Using the default settings of Clustal W in MEGA 8 software, the FLZ protein sequences of Arabidopsis thaliana, rice, Gossypium arboreum, Gossypium raimondii, Gossypium hirsutum and Gossypium barbadense were used for multiple sequence alignment. Based on the results of the sequence alignment, the neighbor-joining method was used to build a phylogenetic tree with the bootstrap value set to 1,000 (Kumar, Stecher & Tamura, 2016). The resulting phylogenetic tree was beautified with the online tool Evolview3 (https://evolgenius.info/).

Chromosomal location, gene structure and motif analysis of the FLZ gene family in Gossypium hirsutum

The chromosomal location information of FLZ gene family members was extracted from the Gossypium hirsutum genome annotation file, and the chromosomal location map of the FLZ genes was drawn by Mapchart software. The evolutionary tree was constructed for the FLZ gene family in Gossypium hirsutum by MEGA8 software, and the nwk file was obtained. Motif analysis was conducted by the MEME program (number of functional domains set to 10, minimum width set to 60, and maximum width set to 100) (Bailey et al., 2009). The XML file, the NWK file of the evolutionary tree and the GFF file of the gene structure were processed and visualized by TBtools software (Chen et al., 2020a).

Analysis of upstream cis-acting elements of the FLZ gene family members in Gossypium hirsutum

The 2,000 bp DNA sequence upstream of the FLZ genes in G. hirsutum was intercepted, and the possible cis-acting elements were predicted using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) and visualized using the R language ggplot package (Lescot et al., 2002).

RNA-seq analysis

Transcriptome data of organs (root, stem, leaf, pistil, stamen, calyx, petal and receptacle) in Gossypium hirsutum under low-temperature stress were downloaded from the NCBI SRA (Sequence Read Archive) database (Genome sequencing project accession: PRJNA248163) (Zhang et al., 2015). Data filtering (i.e., exclusion of reads of too low quality and reads that were too short, cutting of adapters, and trimming of polyX tails in 3′ ends to remove unwanted polyX tailing) and quality control were performed with fastp (https://github.com/OpenGene/fastp) software, and the resulting clean data were used for subsequent analysis (Chen et al., 2018b). The Gossypium hirsutum (TM-1) genome was used as a reference for read alignment (http://www.cottongen.org/), and String Tie was applied to quantify the aligned reads (Kovaka et al., 2019; Kim et al., 2019; Ramirez-Gonzalez et al., 2012; Liao, Smyth & Shi, 2014). FPKM (fragments per kilobase of exon per million fragments mapped) refers to the number of reads per kilobase mapped to an exon in every million reads on a map, and the FPKM method was used to assess gene expression. The expression heatmap was drawn with the R language pheatmap package (Kolde, 2019).

qRT–PCR analysis

According to the cDNA information of G. hirsutum, 5′ and 3′ primers were designed at the specific region of the gene sequence using Primer 5.0 software (Table S1). Root tissue cDNA was used as the template, and the expression of candidate genes was measured by qRT–PCR. Each sample was tested three times, and the internal reference gene was GhUBQ7. qRT–PCR was performed as previously reported (Zhao et al., 2021). A total RNA extraction kit (Tiangen, Sichuan, China) was applied. Reverse transcription was performed using an M-MLV RTase cDNA Synthesis Kit (TaKaRa, Kusatsu, Japan). qRT–PCR was performed using the Roche LightCycler® 480II System under the following conditions: 95 °C 15 s, followed by 40 cycles of 95 °C for 15 s, 55 °C for 15 s, and 72 °C for 15 s. The relative quantitative analysis was performed using the 2^–ΔΔCt method.

Identification of the cotton FLZ gene family

To systematically study the copy number changes of the FLZ gene family during cotton evolution, using the protein domain PF04570, a total of 168 proteins encoding FLZ were identified in four cotton genomes, and the search results were verified in the NCBI-CDD database (Fig. S1). G. arboreum, G. raimondii, G. hirsutum and G. barbadense were found to have 29, 28, 56 and 55 FLZ proteins, respectively. This result indicates that in the process of chromosome doubling and cotton evolution, the FLZ gene family does not exist as a result of gene loss or chromosomal rearrangement. We named the G. hirsutum FLZ proteins GhFLZ1–GhFLZ56 according to the chromosomal locations of the 56 proteins. The length of the open reading frame (ORF) of the G. hirsutum FLZ gene family members ranges from 270 to 1,215 bp, the encoded proteins contain 89–404 amino acid residues, and the differences between the different FLZ proteins are small. The relative molecular weights ranged from 10.18 to 44.64 kDa, and the theoretical isoelectric points ranged from 4.74 to 10.72, indicating that the physicochemical properties of the FLZ gene family members were not very different. The subcellular localization of the proteins showed that 46 were localized to the nucleus and 10 were localized to the chloroplast (Table 1).

To investigate the genomic distribution of G. hirsutum FLZ genes on chromosomes, we investigated the chromosomal location of GhFLZ. Fifty-five GhFLZ genes were distributed on 20 chromosomes of G. hirsutum, and one gene that could not be clearly mapped to any chromosome was named GhFLZ56 (Fig. 1). Subgroup A contained 27 FLZ genes, and subgroup D contained 27 FLZ genes. Previous studies suggested that G. arboreum and Redmond cotton were the donor species of the G. hirsutum A subgenome and D subgenome, respectively, and the number of FLZ genes in G. hirsutum subgroup A was two less than that in G. arboreum subgroup D. The number of FLZ genes in the subgroup was consistent with the number of FLZ genes in Redmond cotton. This indicates that the FLZ gene may have been conserved during the evolution of G. hirsutum, and tandem duplication and that segmental duplication events did not occur. In G. hirsutum, there are three sequences of the FLZ gene family members on chromosome A02 but only one sequence on chromosome D02. There are two sequences of the FLZ gene family members on A03 and A12, D03 and D12 each have three sequences, and the remaining FLZ gene family member sequences on the chromosomes are distributed in stacks. This indicates that the G. hirsutum FLZ gene family may have lost a certain FLZ gene in the process of evolution. However, on the whole, the A subgroup and the D subgroup still have a strong corresponding relationship, which is also in line with the evolutionary relationship of cotton.

Evolutionary analysis of cotton FLZ genes

To further understand the evolutionary relationship of G. hirsutum FLZ, we constructed a phylogenetic tree of the full-length sequences of 56 G. hirsutum FLZ proteins, 27 rice FLZ proteins and 16 Arabidopsis thaliana FLZ protein (Fig. 2A). According to the phylogenetic tree grouping results of Pseudomonas and rice, FLZ proteins were divided into five groups. Each group has at least two FLZ proteins from the monocotyledonous plants Arabidopsis thaliana and rice, indicating that the differentiation time of the FLZ gene family was earlier than that of monocotyledonous plants. Except for Group 3, the number of FLZ proteins of cotton in each subgroup is much greater than that of Arabidopsis and rice, which indicates that the FLZ gene family has undergone obvious tandem duplication in the process of evolution, resulting in tetraploid cotton containing more FLZ genes. Group 3 contains only four FLZ proteins of G. hirsutum but 7 FLZ proteins of rice, which indicates that different subgroups of FLZ may have different biological functions, and the number of different subgroups is very different. It is speculated that different subgroups are responsible for functions that may vary.

To further explore the phylogenetic relationship of the cotton FLZ gene family, an evolutionary tree was constructed using the protein sequences of four different cotton FLZ genes (Fig. 2B). This is consistent with the results of the previous analysis and in line with the evolutionary relationship of cotton. The results showed that the FLZ gene family members are relatively conserved in the evolution of cotton. Although Group 3 members are relatively few, they have always existed in the evolution of cotton, which indicates that they may play an important role in the development of cotton and the process of stress.

To further infer the evolutionary relationship of FLZ genes in G. hirsutum, intrachromosomal collinearity analysis of FLZ genes in G. hirsutum, rice and Arabidopsis was performed using Circos software (Fig. 3A). The results showed that there were a large number of collinear blocks among the three plants (Fig. 3A). The collinearity analysis identified six pairs of FLZ genes in G. hirsutum and Arabidopsis thaliana and 21 pairs of FLZ genes in G. hirsutum and rice. In addition, six G. hirsutum FLZ genes were homologous to both rice and Arabidopsis FLZ genes, indicating that plant FLZ genes may have evolved from a common ancestor of different plants.

According to gene number, chromosomal location and phylogenetic tree analysis, FLZ was found to be conserved in cotton evolution. To deeply study the evolutionary relationship of cotton FLZ, we chose G. hirsutum as the core species and constructed the collinear relationships between G. hirsutum and G. arboreum FLZ proteins and between G. raimondii and G. barbadense FLZ proteins (Fig. 3B). Among them, GhFLZ04 and GhFLZ20 did not have any collinear sequences with any G. barbadense FLZ proteins, GhFLZ19 and GhFLZ30 did not have any collinear sequences with any G. arboreum FLZ proteins, GhFLZ46 did not have any collinear sequences with any G. raimondii FLZ proteins, and GhFLZ43 did not have any collinear sequences with any G. raimondii or G. barbadense FLZ proteins. We found that, in addition to these six genes, the FLZ family gene sequences in the A subgenome of G. hirsutum were collinear with those of G. arboreum and G. barbadense, and the FLZ gene sequences in the D subgenome were colinear with those in G. raimondii and G. barbadense. These results suggest that FLZ was relatively conserved during the evolution of cotton.

Evolutionary tree, gene structure and motif analysis of the FLZ genes of G. hirsutum

The phylogenetic tree, gene structure and motif analyses were performed based on the full-length, CDS and protein sequences of the G. hirsutum FLZ genes (Fig. 4). G. hirsutum FLZ members were divided into five subgroups according to the phylogenetic tree results. All FLZ genes except GhFLZ32 contain two exons, and all contain an identical motif 1. The Group 1 subgroup contained motifs 3 and 9, the Group 5 subgroup contained motifs 3, 5 and 9, and the other subgroups did not contain any motifs. Although the GhFLZ32 gene contains 3 CDSs, it has no redundant motifs, which indicates that the structural differences between GhFLZ32 and genes in the same family may be caused by changes in gene function or errors in genome annotation. In the G. hirsutum FLZ gene family, most members of the same subgroup have similar motifs, lengths and structures, suggesting that they are functionally similar. The protein sequences of the same subgroup were highly conserved, but there were significant differences among the groups, especially for the sequences of Group 1 and Group 5.

Analysis of cis-acting elements of the G. hirsutum FLZ gene promoter

Transcription factors (TFs) can regulate plant growth and development and responses to stress, including responses to hormones and environmental factors, cell differentiation and organ development, by regulating gene expression. To further study the possible transcriptional regulation mechanism of the G. hirsutum FLZ gene family, the cis-acting elements in the upstream 2,000 bp promoter sequences of 56 GhFLZ genes were analyzed (Fig. 5). The results showed that these cis-elements were mainly divided into hormone and stress response elements. The cis-acting elements related to various hormones and biotic stress responses in the GhFLZ gene promoter were summarized, and the results showed that the GhFLZ gene contained at least eight related cis-elements in its promoter region, including five hormone elements and two stress-related elements. In contrast to the evolution results of GhFLZ, the distributions of these two types of elements in the closely related GhFLZ promoter sequences were not similar (Fig. S2). Each FLZ gene promoter contains different numbers and types of cis-acting elements, indicating that they may participate in different biotic and abiotic stress responses through different signaling pathways. Furthermore, of all cis-element types in the GhFLZ genes, SA-responsive elements were the least frequent, with most cis-acting elements being related to MeJA and ABA, indicating that the cotton FLZ gene family may exert its biological functions mainly through the MeJA and ABA pathways. We speculate that external environmental stress can induce the expression of the GhFLZ gene through its response to cis-regulatory elements and further improve the resistance of cotton to environmental stress.

G. hirsutum FLZ gene expression analysis

To elucidate the spatial expression pattern of the G. hirsutum FLZ genes, we analyzed the expression pattern of the G. hirsutum FLZ genes in eight tissues: cotton root, stem, leaf, receptacle, calyx, petal, stamen and pistil. G. hirsutum FLZ genes were all tissue-specifically expressed (Fig. 6A). Three GhFLZ proteins (GhFLZ19, GhFLZ46 and GhFLZ56) were expressed at low levels in all eight tissues, and the expression levels of these three proteins were not much different. Most of the genes were mainly expressed in the roots, stems, leaves, receptacles and calyx. This result indicated that except for GhFLZ19, GhFLZ46 and GhFLZ56, the G. hirsutum FLZ gene family members had strong tissue expression specificity and contained more complex functions. The expression pattern of each GhFLZ gene was shown to be tissue specific in this study, but this specificity may be related to more refined tissues than those examined here.

The expression analysis of FLZ family genes in G. hirsutum under low-temperature stress showed that the expression patterns of all G. hirsutum FLZ family genes can be divided into three categories (Fig. 6B). After low-temperature stress, only 12 genes (GhFLZ03, GhFLZ07, GhFLZ11, GhFLZ24, GhFLZ25, GhFLZ33, GhFLZ37, GhFLZ44, GhFLZ45, GhFLZ50, GhFLZ51 and GhFLZ55) showed evident changes in expression, while the remaining genes did not show visible changes. The expression of these 12 genes was induced by low-temperature stress; thus, these 12 genes may play corresponding roles in the response of G. hirsutum to low-temperature stress. Three genes with obvious expression level changes (GhFLZ50, GhFLZ51 and GhFLZ55) were significantly upregulated and reached maximum levels at 12 h. However, since RNA-seq was only performed on samples collected under low-temperature stress up to 12 h, we speculate that the expression of these genes may have continued to increase after 12 h of low-temperature stress. These results suggest that G. hirsutum FLZ family genes may play a role in the response of G. hirsutum to low-temperature stress.

qRT–PCR

The expression pattern of a gene is usually related to its function. Previous studies have shown that FLZ genes play an important role in plant responses to cold stress (Jamsheer et al., 2018a, 2018b; He et al., 2001). According to the transcriptome expression profile, we speculated that 12 genes may be involved in the stress response of G. hirsutum to low temperature. To understand the expression patterns of these 12 GhFLZ genes in G. hirsutum in response to low-temperature stress, we evaluated their expression patterns in low-temperature-tolerant and low-temperature-sensitive materials by qRT–PCR. Compared with before low-temperature treatment, the expression levels of all but 3 GhFLZ genes (GhFLZ37, GhFLZ50 and GhFLZ51) were significantly different at different time points after low-temperature stress (Fig. 7), indicating that these genes may be involved in cotton’s defense against low-temperature stress. Among them, four genes (GhFLZ11, GhFLZ25, GhFLZ44 and GhFLZ55) were also significantly different between the resistant and sensitive materials at the same time points. In summary, these four genes may all play a role in the response of G. hirsutum to low-temperature stress rather than having naturally high expression during the growth and development of cotton at this stage.

Due to gene expression levels being different in different tissues, the molecular biological functions of a gene product in different tissues are also different. To this end, we selected rhizomes and leaves under low-temperature stress for 24 h to further verify the tissue expression patterns of four genes (GhFLZ11, GhFLZ25, GhFLZ44 and GhFLZ55) in these two materials (Fig. 8). There was no significant difference in GhFLZ25 expression in any tissue of the two materials. The expression of GhFLZ11 in the leaves of the 86-1 cultivar was 1.5-fold greater than that in the Lumian 2 cultivar. The expression levels of GhFLZ44 and GhFLZ55 in the leaves of the 86-1 cultivar were three-fold greater than those in the leaves of the Lumian 2 cultivar, and there were no significant differences in the expression of these genes in other tissues. The expression levels of GhFLZ11 and GhFLZ44 in leaves were significantly higher than those in roots and stems, and the expression levels in stems were the lowest in both low-temperature-resistant materials and low-temperature-sensitive materials. The expression levels of GhFLZ25 in roots and leaves were higher than that in stems. The expression level of GhFLZ55 in leaves was the highest in low-temperature-tolerant materials, but its expression levels in roots, stems and leaves were almost the same in low-temperature-tolerant and low-temperature-sensitive materials, with no significant differences. These results indicate that G. hirsutum FLZ gene expression is likely upregulated in response to low-temperature stress and that the significant changes in G. hirsutum FLZ gene expression in leaves lead to the improvement of cotton cold tolerance. The results of this study show that the expression pattern of the G. hirsutum FLZ gene family is complex, and more in-depth research is needed in the future to analyze the important role of the G. hirsutum FLZ gene family in low-temperature stress. The results of this study also lay a foundation for us to further verify the function and molecular mechanism of GhFLZ proteins and analyze the molecular basis of their role in low-temperature tolerance.