Genome-wide analysis of PRR gene family uncovers their roles in circadian rhythmic changes and response to drought stress in Gossypium hirsutum L.

Identification of PRR gene family in Gossypium spp.

The domain numbered PF00072 (Response receiver domain) and PF06203 (CCT domain) in the Pfam database are often found in plant light signal transduction factors (Sara et al., 2019). Firstly, genome sequence of G. hirsutum (NAU-NBI v1.1 assembly genome), G. arboretum (CRI-updated_v1 assembly genome), G. raimondii (JGI_v2_a2.1 assembly genome) and G. barbadense (ZJU_v1.1 assembly genome) were downloaded from the Cottongen database (http://www.cottongen.org), respectively. This study used the protein sequences of 5 Arabidopsis PRRs were as queries to search the four Gossypium spp. proteomes through the basic local alignment search tool (BLAST, v 2.10.0) with default parameters (E-value = 1 ×10³) for each identified gene (Altschul et al., 1990). PR (Response receiver domain) and CCT domains, the typical PRRs domains, were aligned and searched in HMMER 3.0 (https://www.ebi.ac.uk/Tools/hmmer/) (Potter et al., 2018). Next, sequences were searched and verified on the Conserved Domain Database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART (http://smart.embl-heidelberg.de/) (Letunic et al., 2002). Finally, the online site ExPASy Proteomics Server (http://www.expasy.org/) and Softberry (http://linux1.softberry.com/berry.phtml?topic=protcomppl&group=programs&subgroup=proloc) were used to analyze the physicochemical properties of the identified cotton PRR gene family, including amino acid number, nucleotide data, molecular weight, isoelectric point prediction and subcellular localization.

Chromosomal locations, duplications, and synteny analysis of PRR gene members

Chromosomal location information for PRR genes was obtained from general feature format (gff) files of each cotton genomic databases and genes were mapped on the chromosomes using TBtools (Chen et al., 2020). Then MCScanX (Wang et al., 2012a; Wang et al., 2012b) was used to determine and analyze cotton PRR duplication and collinearity, Circos (http://circos.ca/) software were used to conducted image showing gene location and gene homology relationship.

Phylogenetic analyses and gene structure organization of the PRR proteins in Gossypium spp.

To analyze evolutionary relationship, the PRR proteins sequence of various plant species including Arabidopsis thaliana (Initiative, 2000), Cocoa (Theobroma cacao) ((Argout et al., 2011)) and rice (Oryza sativa) (Yu et al., 2005) were downloaded from the Arabidopsis database TAIR10 (https://www.arabidopsis.org/), the plant genome database Phytozome 12 (http://phytozome.jgi.doe.gov/pz/portal.html) and EnsemblPlants (http://plants.ensembl.org/index.html), respectively. Multi-protein sequence alignment of the PRR proteins were aligned using MEGA7.0 (Sudhir, Glen & Koichiro, 2016), and constructed a phylogenetic tree using neighbor-joining (NJ) method with the bootstrap 1000. Finally, the evolutionary tree is visualized and beautified by the online software iTOL (https://itol.embl.de/) (Letunic & Bork, 2019). Location information of PRR members were obtained from gff files using SeqHunter1.0 (Ye, Wang & Dou, 2010) and the gene structures were displayed by the online software Gene Structure Display Server (GSDS 2.0) (http://gsds.gao-lab.org/) (Guo et al., 2007), and we performed motifs analysis on the online software MEME (http://meme-suite.org/) (Bailey et al., 2009) with following parameters: the maximum number discovered for the motif is 10, and the other parameters are default values. The graphic display is based on the Amazing optional gene viewer section in the software TBtools.

To identify the cis-elements in the promoter sequences of the 16 PRR family genes in G. hirsutum, the 2,000 bp of genomic sequences upstream of the start codon for each PRR gene were submitted to the online site PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), and the results are displayed by the Simple Bio Sequence Viewer in TBtools.

Plant materials and treatment

The upland cotton (G. hirsutum) accession (Lumianyan 19, LMY 19) (Li et al., 2004), an early maturing variety, selected in this study were kept in our laboratory, planted in growth chamber (day/night temperature cycle of 28 °C light/25 °C dark with a 12-photoperiod), and samples were picked every 3 h from leaf in three-true-leaves stage. Germinated TM-1 cotton seeds were planted in the same photoperiod and temperature environment as LMY19, and treated with 400 mM polyethylene glycol (PEG6000) at the three-leaf stage from 7 am. TM-1 seedlings were divided into four treatment groups, treated with PEG for 0, 1, 3 and 6 h, respectively, and non-PEG treated seedling (treated with sterilized water) as control check at the same time point. Then we collected leaf samples at 0 h, 1 h, 3 h and 6 h after PEG treatment and non-PEG treated samples at each of the time point. Three biological replicates for each sample, the leaves from three seedlings as a biological replicate, and all samples were freezed with liquid nitrogen immediately and stored at −80 °C for qRT-PCR.

RNA isolation and qRT-PCR analysis

The RNA was extracted from the samples using the Rapid Universal Plant RNA Extraction Kit (Huayueyang Biotechnology Co. Ltd.), and the Prime Scrip First Strand cDNA Synthesis Kit (Takara) used for reverse transcription, SYBR Premix Ex Taq II. (Takara) kit used for real-time PCR experiment, qRT-PCR analysis was carried out using SYBR Green on the Roche LightCycler^® 480 II. The primers of PRR gene family were designed using Primer Premier 5.0 software and listed in Table S1, and the actin gene (AF059484) was selected as the internal reference gene (Zhang et al., 2013; Zhang et al., 2013b). The volume of the qRT-PCR reaction was 20 µL, and the amplification procedure was as follows: pre-denaturation at 95 °C for 30 s; denaturation at 95 °C for 5 s, annealing at 60 °C for 30 s, 40 cycles. Three biological and technical replicates were performed for the qRT-PCR tests. The relative gene expression levels were quantified by the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Expression patterns and pathway enrichment analysis of PRR members

RNA-Seq data of G. hirsutum TM-1were obtained from the SRA database (PRJNA248163) (Zhang et al., 2015), and the FPKM (fragments per kilobase per million reads) values were calculated by RNA-seq data downloaded from the database of cottonFGD (Zhu et al., 2017).The gene expression pattern of PRR genes were displayed by R/pheatmap with the expression values normalized by log₂(FPKM+1). The expression profiles of all 16 GhPRR genes at different time of PEG treatment were further analyzed using R/Mfuzz. The differentially expressed genes (DEGs) were identified by DEseq2 (Anders & Huber, 2010). All detected genes in each sample were used to identify significantly DEGs (—log2 Foldchange—>1, P < 0.05) and KEGG analyses of DEGs were conducted in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database for enrichment (Kanehisa et al., 2014), KEGG enrichment of DEGs was evaluated with KOBAS2.0 software (Xie et al., 2011) and bubble graph was displayed by R/ggplot2.

Genome-wide identification of PRR family genes in Gossypium spp

Based on multiple sequence alignment analysis, complete PRR genes were identified in four Gossypium species, including 16 in G. hirsutum (AD₁), 9 in G. arboretum (A₂), 10 in G. raimondii (D₅), and 9 in G. barbadense (AD₂). Additionally, we proceeded with PRR genes retrieved from plant genome database, with 5 in Arabidopsis (dicots), 5 in rice (monocots), and 6 in cocoa (dicot). All of them were renamed based on the homologous genes in Arabidopsis (Table S2). The number of PRR gene family in G. hirsutum (AtDt) was about twice as that in G. arboreum (A group) or G. raimondii (D group), it is consistent with the former one being tetraploid and the latter two being diploid. The basic information of PRR genes including protein sequence length, isoelectric points, and molecular weight in cotton were listed in Table S3 . The predicted GhPRR proteins ranged from 552 (GhPRR1a and GhPRR1b) to 795 (GhPRR3a) amino acids, with isoelectric points changed from 4.97 (GhPRR9a) to 8.42 (GhPRR3d) and molecular weight from 61.87 kDa (GhPRR1a) to 85.93 kDa (GhPRR3a).

Chromosomal locations, duplications, and synteny analysis of PRR gene members

In order to display the chromosome distribution of PRR genes, mapping them on the corresponding chromosome. Eight of GhPRR genes were located on chromosomes of At sub-genome while five of GhPRR genes were on that of Dt sub-genome and three GhPRR genes were present in different scaffolds (Fig. S1). We further conducted whole genome collinearity analysis of 44 identified PRR genes in cotton, and explored the locus relationships between At and Dt sub-genomes as well as with A and D diploid cotton genomes (Fig. 1A, Table S4). There are 34 orthologous gene pairs were resulted from whole genome duplication or segmental duplication among Gossypium spp. Whole Genome duplication or segmental duplication was suggested to be the main causes of PRR gene family expansion in cotton (Table S5).

Phylogenetic analyses and gene structure organization of the PRR proteins in Gossypium spp.

To investigate the evolutionary relationship of GhPRR proteins among mentioned seven species, phylogenetic tree was constructed (Fig. 1B). The PRR family of Gossypium was divided into five subgroups (clade I–V). There were 13 PRRs in Clade III (three GaPRRs, GbPRRs and GrPRRs respectively, four GhPRRs) and 11 PRRs (one GrPRR, two GaPRRs, four GbPRRs and GhPRRs individually) in clade IV. Clade I consisted of 9 PRRs (one GaPRR, two GbPRRs and GhPRRs singly, four GrPRRs), Clade V contained 7 PRRs (one GrPRR, two GaPRRs and four GhPRRs) and Clade II had 4 PRRs (one GaPRR and GrPRR respectively, two GhPRRs). GhPRRs were distributed throughout five subgroups (clade I–V), clade-I, clade-II and clade-IV containing PRRs from monocots and dicots simultaneously, illustrating that evolution of GhPRR genes in three clades occurred before the separation of monocots and dicots.

PRRs protein in G. hirsutum was also divided into five subgroups (Fig. 2A), consistent with phylogenetic analyses. The motif distribution indicated that the order, size, and location of the motifs in the same subgroup were similar, but there were significant variety between different subgroups. Among them, 37.5% of the family members have the same sequence of motif structure: motif 4_9_3_1_7_5_6_10_8_2, while Clade-I contains the least number of motifs with only 5 motifs. All members of the PRR gene family contain motif1, motif2, motif3, motif4 and motif6, which are the conserved motifs of PRR family. In addition, the gene structure analysis exhibited that the distribution of introns and exons were similar among different subgroups, and the functional elements PR and CCT were distributed in both end side of each gene (Fig. 2B). All of member contained three PR structure elements, and most member contain two CCT domains, except that two members of the Clade-I subgroup contain one CCT domain.

To further analyze the transcriptional regulation and potential function of the PRR genes, the cis-elements in the promoter region were predicted (Fig. 2C). The results displayed that there are abundant regulatory elements existing in the promoter region, mainly focused on light response elements (G-Box, GT1-motif and TCT-motif, etc.), hormone responsive elements: abscisic acid response (ABRE), MeJA-response (CGTCA-motif and TGACG-motif), gibberellin-responsive element (TATC-box, P-box and GARE-motif), and stress responsive elements: drought-inducibility (MBS), low-temperature response (LTR), etc. There are 16, 14 and 6 PRR genes containing response elements to light, abscisic acid and drought stress, respectively. Motif sequences are often the binding sites of some sequence-specific proteins (such as transcription factors), have important biological significance for important biological processes, such as RNA initiation, RNA termination, RNA cleavage, etc.

The expression pattern of PRR members under diurnal changes

A feature shared by many clock gene transcripts is that their abundance is subject to diurnal oscillation. To analyze the peak transcripts of GhPRR s under diurnal cycle, the relative expression levels of GhPRRs together with its related genes (GhFT (FLOWERING LOCUS T), GhCO (CONSTANS LIKE -2), GhLHY (LATE ELONGATED HYPOCOTYL) and GhCCA1 (CIRCADIAN CLOCK-ASSOCIATED 1)) during 24 h was detected by qRT-PCR (Fig. 3 and Table S6). The results showed that GhLHY-mRNA began to accumulate after dawn, and then mRNA of GhPRR genes began to reach the peak sequentially within a 24-hour period with multiple members at each peak. GhFT, GhCO, and GhLHY had the peak expression at 3 h after light. Subsequently, members inclade-II (GhPRR9a and GhPRR9b) and clade-IV (GhPRR7a, GhPRR7b, GhPRR7c and GhPRR7d) reached the expression peak after 6 h of light condition, and then members in clade III (GhPRR5a, GhPRR5b, GhPRR5c and GhPRR5d) and clade-V (GhPRR3a and GhPRR3c) at 9 h, another two members of clade-V (GhPRR3b and GhPRR3d) at 12 h. Finally, members (GhPRR1a and GhPRR1b) in clade-I reached expression peak after 3 h of dark. Additionally, the expression of GhLHY and GhPRR1b always showed an opposite trend during 24 h, it can be speculated that a mutual inhibition maybe exist between the two genes. These results indicated that expression of GhPRR genes has obvious rhythmic expression trend waves during 24 h.

Identification of drought-stress related PRR genes in G. hirsutum

To investigate the roles for PRR genes in response to drought, we investigated the expression profile of GhPRRs under polyethylene glycol (PEG) treatment at 1, 3 and 6 h from the published transcriptome data sets. All detected genes in each sample were used to identify significantly DEGs (—log2 Foldchange—>1, P < 0.05) among PEG_1 h vs CK, PEG_3 h vs CK, PEG_6 h vs CK groups, and the PEG_6 h group contains the most number of DEGs (Table S7), so we selected the group data at 6 h treated with PEG for KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis (Fig. 4A).The results revealed that the DEGs are mainly involved in circadian rhythm, photosynthesis, starch and sucrose metabolism, etc. (Fig. 4B). Six of GhPRR genes including three members in clade III (GhPRR5a, b, d) and two in clade-V (GhPRR3a and GhPRR3c) were involved in circadian rhythm pathway. The expression patterns of these GhPRR genes have high expression level at 6 h with PEG treatment (Fig. 4C and Table S8), further analyzed and divided into 3 clusters, three members of clade III (GhPRR5a-d) and two of clade-V (GhPRR3a and GhPRR3c) in Cluster1 exhibited the same expression trend (Fig. S2 and Table S9), suggesting these PRR genes are significantly induced by PEG treatment. .

To further prove the expression changes of these genes at different time of PEG treatment (0 h, 1 h, 3 h and 6 h), the expression level of all member of PRR family were detected by qRT-PCR (Figs. 5A–5P and Table S10). The expression of genes (GhPRR3a, c and GhPRR5a, b, d) at the sixth hour after PEG6000 treatment was significantly higher than that of the blank control. All PRR genes displayed almost the similar expression changes compared with transcriptome data sets (CK, 1 h, 3 h, 6 h), and the correlation analysis between the transcriptome and qRT-PCR of GhPRR genes displayed by scatter plots, the result showed that the Pearson correlation coefficient log2 expression ratios calculated from qRT-PCR and RNA-seq of GhPRR genes was 0.78 (Fig. S3), suggesting the results are credible. It can be considered that genes mentioned above (GhPRR3a, c and GhPRR5a, b, d) maybe respond to drought stress.