Genome-wide investigation and expression pattern of PHR family genes in cotton under low phosphorus stress

Identification and analysis of cotton PHR gene family members

To obtain the members of PHR family genes, the amino acid sequences with conserved domain Myb DNA-binding (PF00249) and Myb_CC_LHEQLE (PF14379) of the PHR transcription factor family were searched from the four Gossypium genomes including Gossypium arboreum (CRI, version 1.0), G. raimondii (JGI, version 2.0), G. hirsutum(HEBAU, version 1.0) and G. barbadense (HEBAU, version 1.0) (Paterson et al., 2012; Wang et al., 2012; Li et al., 2014; Ma et al., 2021). The obtained sequences of PHR transcription factor were identified using SMART (http://smart.embl-heidelberg.de/) and CDD databases (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). The physical and chemical properties of the proteins were predicted using the online website ExPASY (https://web.expasy.org/protparam/).

Conserved domain and gene structure analysis of cotton PHR family members

The online website MEME (http://meme-suite.org/tools/meme) was used to predict the conserved motif of cotton PHR genes, and the number of motifs was set to 10, other parameters are the default settings. The genes structure of cotton PHR transcription factors were draw using online website GSDS (http://gsds.gao-lab.org/) based on cotton genome annotations.

Phylogenetic relationship of cotton PHR gene family members

The protein sequences of Arabidopsis thaliana PHR family gene were downloaded from Phytozome database (https://phytozome-next.jgi.doe.gov). The protein sequences of PHR family genes in cotton and Arabidopsis were compared by MEGA 11.0 software (Koichiro, Glen & Sudhir, 2021), and the phylogenetic tree was constructed by adjacency method (NJ). The bootstrap value was 1,000, the model is Poisson model. Finally, we showed the results using the online tool iTOL (Letunic & Bork, 2021).

Analysis of cis-acting elements of promoters of GhPHR gene family members

To understand the possible regulation and response mechanism of cotton PHR genes, the 1.5 kb upstream of genome sequences were obtained from each GhPHRs, and submitted to PlantCARE website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to predict cis-acting elements of promoters, and finally visualized by TBtools software (Chen et al., 2020).

Prediction of miRNA-target GhPHR genes

According to the principle of sequence complementarity, the regulatory miRNAs of GhPHR genes were predicted. The GhPHR target genes of miRNA were predicted using psRNATarget online software (https://www.zhaolab.org/psRNATarget/).

Differential expression analysis of GhPHR family members in roots under low phosphorus stress

Tissue expression specificity of GhPHR genes were analyzed by using RNA sequencing data of G. hirsutum ‘TM-1’ (Zhang et al., 2015) during different development stages and stress treatments downloaded from the NCBI Sequence Read Archive (PRJNA490626). To further screen GhPHR genes in response to low phosphorus stress, we generated and analyzed the genes expression in roots based on transcriptome data of a P-resistant accession under P deficient hydroponic conditions. The cotton seeds were grown in germination boxes containing quartz sand, and the seedlings were moved into half-strength Hoagland normal nutrient solution when the cotyledons had fully expanded. After 3 days, the nutrient solution was replaced with P-deficient and P-replete Hoagland nutrient solution, respectively. The roots of three cotton seedings were sampled at 0 day and 15 days for RNA-seq. The gene expression patten was drawn by Hemi 1.0 software (Deng et al., 2014) with log2 (1 + FPKM) values after averaging three replicates. The expression of six GhPHR genes were selected for further verification by qRT-PCR. All reactions were performed using the Roche LightCycler96 RealTime PCR System with three independent biological replicates, each with three technical replicates. The expression of GhUBQ14 was used as the internal control. Relative gene expression values were calculated using the 2^−ΔCT method (Schmittgen & Livak, 2008).

Identification of cotton PHR family gene members

To investigate the copy number variation of the PHR genes during cotton evolution, a comprehensive search was conducted for PHR genes across four cottons including G. arboreum, G. raimondii, G. hirsutum and G. barbadense. A total of 128 PHR gene sequences were detected in the four cotton species, 22, 23, 41, and 42 PHR genes were identified, respectively. The detailed information was listed in Tables S1 and S2. The results were further verified through the NCBI-CDD and SMART database. The results showed that the numbers of PHR genes in two diploid cotton species were almost similar as were those in two tetraploid cotton species. The number of PHR family genes in G. arboretum and G. raimondii were basically half in G. hirsutum and G. barbadense, indicating that the PHR family was relatively conserved, which conforms to the known evolutionary relationship in different cottons (Paterson et al., 2012).

The names of the PHR genes were determined according to the gene information of Arabidopsis and the locations on the chromosomes. The encoded protein of the PHR family genes in G. hirsutum contains 236~494 amino acid residues. The relative molecular mass is between 26.53 and 54.30 kDa, and the theoretical isoelectric point is between 5.50 and 9.82. Each of the family members contains Myb DNA-binding and Myb_CC_LHEQLE domain. Fourty-one PHR genes of G. hirsutum were distributed on 21 chromosomes except A05, A07, D01, D03 and D07 (Table 1). Subgenome A and subgenome D contained 22 and 19 sequences, respectively. The genes number in subgenome A was consistent with the genes number of GaPHR, and four GhPHR genes from subgenome D were missing compared to GrPHRs. This result indicated that subgenome D of G. hirsutum might have lost genes due to redundant gene functions during evolution.

Three GhPHR genes were observed on chromosomes A09 and A13, while D09 and D13 only contained two genes. The A01 and A03 chromosomes contained one gene, but no GhPHR gene sequence was contained in D01 and D03. These results showed that the GhPHR genes might have been lost and duplicated in the process of cotton evolution, indicating a strong correlation between subgenome A and subgenome D (Paterson et al., 2012; Wang et al., 2018).

Phylogenetic analysis of the PHR gene family in cotton

To explore the phylogenetic relationship of the cotton PHR genes, we constructed a phylogenetic tree with the neighbor-joining method using 41 G. hirsutum, 42 G. babardence, 22 G. arboreum, 23 G. raimondii and 13 Arabidopsis PHR amino acid sequences (Fig. 1). All the PHR genes can be divided into five subgroups. The PHR genes number of G. hirsutum and G. barbadense was basically twice of G. arboreum and G. raimondii in each subgroup. This result was consistent with the previous analysis and conformed to the conserved evolutionary relationship in cottons. Among these subgroups, the largest subgroup V consisted of 12 GhPHRs, 12 GbPHRs, seven GaPHRs, seven GrPHRs and six AtPHRs, showing that has expanded considerably in two allotetraploid cottons. In contrast, subgroup III only included four GhPHRs, four GbPHRs, two GaPHRs, two GrPHRs and one AtPHRs, indicating might be a highly-conserved clade.

Analyses of gene structures and protein motifs of PHR genes in G. hirsutum

Gene structure analysis of PHR gene family members showed that the gene structure of GhPHR members in the same subgroup was similar, and there was little difference in the number of exons of GhPHR members among different subgroups. Except for the GhPHR genes of class V which contained seven exons, most of the other GhPHR genes contained six exons (Fig. 2). Furthermore, these PHR protein sequences were submitted to MEME to discover conserved motifs. The adjacent clades carried similar motifs. Analysis of the conserved domains of GhPHR gene family members showed that Myb_DNA-binding and Myb_CC_LHEQLE were present in all GhPHR proteins and other motifs were functionally unknown motifs (Fig. 2). Among these, motif 4 occurred in class I, class II and class IV subgroups, motif 6 was unique to class II subgroup, motif 10 were unique to class IV subgroup, and motif 5 was unique to two of class I class II subgroup. It is worth noting that the GhPHR protein motif types of class II were different. Among them, there were eight motifs in GhPHR2, GhPHR15 and GhPHR35, only six motifs in GhPHR5 and GhPHR26, and seven motifs in the other five GhPHRs. These special conserved motifs may be the main factors enabling GhPHRs to participate in different biological functions.

Analysis of cis-acting elements in the promoter of GhPHR family genes

To further clarify the possible regulatory mechanism of GhPHR family genes under abiotic stress, the promoter sequences were analyzed by using PlantCARE database. The results showed 13 types of cis-acting elements, including light responsiveness, salicylic acid responsiveness, gibberellin-responsiveness, MeJA-responsiveness, anaerobic induction, auxin-responsiveness, abscisic acid responsiveness and so on (Fig. 3). In terms of composition and quantity, the GhPHR genes contained an average of 18 cis-acting elements, all of which had light responsiveness elements. Among them, GhPHR6 and GhPHR19 contained the most types of response elements (10 kinds), while GhPHR3, GhPHR5, GhPHR7 and GhPHR40 contained the least cis-acting elements (three kinds). In terms of element types, 17 genes contained gibberellin—responsiveness elements including GhPHR1, GhPHR12 and GhPHR34. Thirty-one genes contained anaerobic induction elements including GhPHR5, GhPHR19, GhPHR24, and 11 genes contained auxin-responsiveness elements including GhPHR16, GhPHR20, GhPHR36. The results showed that the GhPHR genes were not only regulated by light induction, but also played a role during drought, anaerobic and other stresses. Among them, the promoter regions of GhPHR3, GhPHR5, GhPHR7 and GhPHR40 contained less cis elements, which were only related to light, MeJA, abscisic acid and anaerobic induction.

Prediction of the regulatory miRNA of PHR gene family

The online software psRNATarget was used to predict and analyze the regulatory miRNAs of GhPHR genes. The regulatory combinations of 33 miRNAs and PHR genes were predicted (Table S3). It was found that 12 miRNAs could regulate 16 PHR genes. Unpaired energy (UPE) is the energy required to unlock the secondary structure of the target gene miRNA target site. A lower UPE value indicates that miRNA is more likely to bind or cleave the target gene. This study showed some of the predicted results with an expected value less than or equal to 5. GhPHR21 and GhPHR32 can be recognized by miR396 and miR7510a, miR2949 and miR7491 at the same time, respectively, which may be regulated by these two miRNAs. GhPHR19 may be regulated by the sequence cleavage of miR482, GhPHR4 and GhPHR24 may be regulated by the transcriptional inhibition of miR827, GhPHR17 and GhPHR37 may be regulated by the sequence cleavage of miR2948-5p, and GhPHR23, GhPHR25 and GhPHR40 may be regulated by the transcriptional inhibition of miR2948-5p. In this study, the regulation modes of interaction combinations were different. Approximately 2/3 of the regulation modes belonged to sequence cleavage and 1/3 belonged to transcriptional inhibition (Table 2).

Tissue expression analysis of GhPHR family genes

The expression of GhPHR family genes was investigated across different tissues and developmental stages of upland cotton. Most of these genes were expressed at varying levels across different tissues and developmental stages (Fig. 4). It was found that five GhPHR genes were the most commonly expressed in all tissues. A total of 13 GhPHRs were highly expressed in all tissues, indicating that these genes play an important role in all morphogenesis of cotton. Among them, GhPHR5 and GhPHR6 were highly expressed in leaves, while GhPHR13, GhPHR14, GhPHR26 and GhPHR27 were highly expressed in both roots and leaves. Combined with the cis-element structure of GhPHR promoters, it was speculated that they may play an important role in leaf photosynthesis. A total of eight GhPHRs were expressed across different tissues except fiber developmental stages, and the expression of other six GhPHR genes were low in different tissues. Altogether, the expression profiles of GhPHR genes showed that it plays a role in different tissues of cotton, among which GhPHR2, GhPHR5, GhPHR6 and GhPHR13 have obvious tissue expression specificity.

Expression pattern of GhPHR genes in roots under low phosphorus stress

The expression of 41 GhPHR genes in roots showed that most GhPHR genes were affected under low phosphorus treatment, except that GhPHR4, GhPHR7, GhPHR12, GhPHR19, GhPHR24, GhPHR33 and GhPHR39 were not detected (Fig. 5). The expression of GhPHR1 and GhPHR11 decreased under low phosphorus stress but that of GhPHR3, GhPHR6, GhPHR17, GhPHR18, GhPHR27, GhPHR30 and GhPHR38 increased. Among them, expression level of GhPHR17, GhPHR30 and GhPHR38 was significantly higher than that before stress treatment. In addition, GhPHR5, GhPHR13, GhPHR15 and GhPHR26 maintained a high expression level. It should be noted that the expression of GhPHR11 was significantly lower under low phosphorus stress than under normal phosphorus treatment, and that of GhPHR18 was significantly higher than that of normal phosphorus treatment (Fig. 5). Further, we selected six differentially expressed GhPHR genes for verification by qRT-PCR, showing a consistent trend (Fig. 6). This provided further evidence that the six putative genes were closely associated with low-phosphorus tolerance.