Comprehensive Analyses of the PfGRF Transcription Factor Family and Its Response to Biotic and Abiotic Stresses

1 Introduction

Paulownia species, which are important fast-growing trees, are cultivated worldwide because of their high ecological, economic, and medicinal value [1,2,3]. Paulownia witches’ broom (PaWB) is an infectious disease caused by phytoplasma, with symptoms that include witches’ broom, stunting, short internodes, yellowing leaves, decreased leaf area, and death [3,4]. According to earlier reports, PaWB in China results in annual economic losses of billions of dollars [5,6].

The rapid development of high-throughput sequencing technology and the maturation of molecular biology research methods have enabled researchers to analyze uninfected Paulownia seedlings and phytoplasma-infected Paulownia seedlings in terms of mRNA expression, post-translational modifications, metabolomes, and epigenetic changes [7,8,9]. Many genes and proteins related to PaWB were identified in previous studies [4,7,8,9,10]. However, the molecular mechanisms underlying PaWB in adult Paulownia trees remain unknown. Thus, the molecular basis of PaWB will need to be more comprehensively characterized.

Transcription factors (TFs), which are encoded by the most important regulatory genes in plants, are involved in many biological processes influencing plant growth, development, metabolism, reproduction, and stress resistance [11,12,13,14]. Currently, more than 60 TF families have been identified in plants, of which the growth regulating factor (GRF) TF family is specific to plants, wherein it plays an important regulatory role affecting growth and development, flower organ development, and stress responses [15,16,17]. The first GRF TF identified in rice (OsGRF1) reportedly contributes to gibberellin mediated stem elongation [18]. Silencing OsGRF3, OsGRF4, and OsGRF5 expression via RNA interference retards growth, resulting in stunted rice plants [19]. The heterologous expression of maize ZmGRF in Arabidopsis thaliana (A. thaliana) leads to cell expansion and stem elongation, GA4 accumulation (3.7–5.7 fold), up-regulated expression of a GA-receptor gene (GIF), and down-regulated expression of a GA-insensitive growth suppressor DELLA protein-encoding gene. Thus, ZmGRF functions through the GA pathway [20]. Because an increasing number of genomes have been analyzed, GRF genes have been identified in several species, including A. thaliana [21], maize [22], and humans [23]. Therefore, the mechanisms mediating the effects of GRF on plant stress resistance should be clarified. Unfortunately, there are few reports describing research on GRF functions and their potential regulatory effects on Paulownia stress resistance.

In this study, on the basis of published Paulownia genome information, transcript and miRNA sequencing technology were employed to analyze gene expression changes in 10-year-old witches’ broom-resistant (WPF) and witches’ broom-infected (WPFI) P. fortunei (Seem.) Hemsl. plants growing under natural conditions. PaWB-responsive genes were identified in P. fortunei, and GRF family genes were further screened and analyzed. Notably, the GRF transcription factor family has not been systematically identified or functionally investigated in any species of Paulownia, representing a critical research gap. Hence, a bioinformatics analysis of the PfGRF gene family is important for further determining PfGRF functions. In-depth research on the response of PfGRF genes to PaWB will help clarify plant defense mechanisms and the GRFs related to disease resistance, thereby providing the basis for exploiting GRF TFs and the breeding of PaWB-resistant plants through genetic modifications.

2 Materials and Methods

2.1 Plant Materials

In this study, 10-year-old WPF and WPFI P. fortunei plants growing in the experimental field of Henan Agricultural University were used as experimental materials. Buds were collected from samples, with three biological replicates per sample. The collected buds were immediately frozen in liquid nitrogen and then stored at −80°C prior to subsequent analyses.

The phytoplasma in WPF and WPFI samples were determined following the method described by Yang et al. (2023) [7].

2.2 Transcriptome and miRNA Sequencing Analysis of P. fortunei Infected with Phytoplasma

2.2.1 Total RNA Extraction and Detection

Total RNA was extracted from all samples using an RNAprep Pure Plant Kit (TIANGEN, Beijing, China). The mass and concentration of the extracted RNA were determined using a 2100 Bioanalyzer (Agilent, CA, USA) and an RNA 6000 Nano Lab Chip Kit (Agilent).

2.2.2 Database Construction and Sequencing

The Epicentre Ribo-Zero Gold Kit (Illumina, CA, USA) and TruSeq Small RNA Sample Prep Kits (Illumina) were used to construct strand-specific libraries (>200 nt) and small RNA libraries (<50 nt), respectively. An Illumina HiSeq 4000 system was used to detect mRNA via chain-specific library sequencing, whereas an Illumina HiSeq 2500 system was used to conduct a miRNA library sequencing analysis. Details regarding mRNA and miRNA were obtained from the data on the basis of biogenic analysis. Additionally, CPC, CNCI, and txCdsPredict software as well as the Pfam database were used to predict the coding potential of the identified new transcripts.

2.2.3 Transcriptome Sequencing and Analysis

The FPKM method was used to calculate gene expression. The Audic–Claverie algorithm was used to identify genes that were differentially expressed between WPF and WPFI. The criteria for identifying significant differentially expressed genes (DEGs) were as follows: false discovery rate < 0.05 and |log2(fold-change)| > 1. Finally, the DEGs and non-coding genes were functionally characterized via GO and KEGG pathway analyses.

2.2.4 Identification of miRNAs and Prediction of Their Target Genes

ACGT101-miR (Houston, TX, USA) was used to analyze miRNA data, whereas Target Finder was used to predict miRNA target genes. The functions of the predicted target genes were analyzed according to GO and KEGG analyses as previously described [21].

2.2.5 Correlation Analysis of miRNA and mRNA

A miRNA–mRNA association analysis was completed using Perl according to a published method [24].

2.3 PfGRF Transcription Factor Family in P. fortunei

2.3.1 Identification of GRF Genes in P. fortunei

Whole-genome data for P. fortunei were obtained from the NCBI database [2]. The Hidden Markov Model (HMM) file of WRC (PF08879) and QLQ (PF08880) were downloaded from the Pfam database (https://pfam.xfam.org/). The WRC and QLQ domain sequences were used as queries to screen the P. fortunei protein dataset for proteins containing these domains using the HMMER 3.0 program, with the threshold set at e < 1 × 10−5 [21]. Protein sequences encoded by A. thaliana GRF genes were downloaded from an A. thaliana database (https://www.arabidopsis.org) and then used as queries to search the P. fortunei protein dataset using BLASTP, with the thresholds set at e < 1 × 10−5 and 50% identity. The results obtained using these two methods were compared and analyzed to identify PfGRF family members, which were named according to genome information. The Conserved Domain Database (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and the Pfam database were used to confirm the identified PfGRFs, which were subsequently characterized in terms of their molecular weight (MW), isoelectric point (pI), and grand average of hydropathicity (GRAVY) value using the ExPASy server [23] (https://web.expasy.org/compute_pi/). Moreover, their subcellular localization was predicted using WoLF PSORT (http://www.genscript.com/psort/wolf_psort.html).

2.3.2 Phylogenetic Tree, Conserved Motif, and Structural Analyses of PfGRF Genes

The amino acid sequences of the confirmed PfGRFs were aligned using ClustalW. A neighbor-joining phylogenetic tree was constructed using MEGA 7.0, with 1000 bootstrap replicates [25]. Conserved motifs were identified using the online MEME software (http://memesuite.org/tools/meme) and the following parameters: arbitrary number of repeats, up to 10 motifs, and motif width of 6–200 [24]. A PfGRF gene structural analysis was performed using the P. fortunei genome database. The results of these analyses were visualized using TBtools.

2.3.3 Chromosomal Distribution and Collinearity Analysis of PfGRF Genes

The locations of PfGRF genes on chromosomes were determined using P. fortunei genome information. TBtools was used to find potential homologous gene pairs, identify syntenic chains, and determine the types of duplication mechanisms in the P. fortunei genome as well as for visualizing the results and calculating nonsynonymous (Ka) and synonymous (Ks) substitution rates among the PfGRF genes [24].

2.3.4 Phylogenetic Tree and Collinearity Analyses of GRF Genes in P. fortunei and Other Species

A. thaliana genome data were downloaded from the TAIR database (https://www.arabidopsis.org), whereas Oryza sativa genome data were downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/). An evolutionary tree was constructed using MEGA 7.0 software, while a collinearity analysis of GRF genes in P. fortunei, O. sativa, and A. thaliana were performed using TBtools [26].

2.3.5 Analysis of PfGRF Promoter Cis-Acting Elements

The genomic sequence 2 kb upstream of the start codon of each PfGRF gene was considered as the promoter region. The PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used to predict cis-acting elements in the promoter region [27].

2.4 PfGRF Transcription Factor Responses to Biotic and Abiotic Stresses

RNA sequencing (RNA-seq) data for the PfGRF genes in the WPF and WPFI plants were downloaded from the NCBI Sequence Read Archive (SRA) database (SRA accession numbers: SRR11787883, SRR11787894, SRR11787905, and SRR11787912–SRR1178792). A heatmap of PfGRF expression was generated using TBtools [26].

2.5 Y2H and BiFC Verified Protein Interaction

Methodological protocols for protein interaction validation via Y2H and BiFC, as described in Yang et al. (2023) [7]. All primers used are listed in Table S1.

3 Results

3.1 Transcriptome Sequencing Analysis of P. fortunei Infected with Phytoplasma

The results showed that phytoplasmas were detected in WPFI by quantitative real-time PCR, while no phytoplasmas were found in WPF (Fig. S1).

A total of 764,557,918 clean reads were obtained from the RNA-seq analysis of the following six cDNA libraries: 97.90% (WPF-1), 97.77% (WPF-2), 97.78% (WPF-3), 98.01% (WPFI-1), 98.08% (WPFI-2), and 98.20% (WPFI-3) (Table S2). The size distribution of mapped reads for the identified mRNA sequences is presented in Fig. S2. More than 70% of the mapped mRNAs were 0–1500 nt long (Fig. S2a), with 1 being the most common number of exons (Fig. S2).

To determine the transcriptional changes in WPF and WPFI, we identified DEGs through comparisons (7036 and 4897 genes with increased and decreased expression levels, respectively, in response to PaWB) (Fig. 1A). The main enriched KEGG pathways among these genes were biosynthesis of secondary metabolites, plant hormone signal transduction, plant–pathogen interaction, and starch and sucrose metabolism. Hence, the genes contributing to these metabolic processes are likely closely related to the occurrence of PaWB (Fig. 1B).

Figure 1: PaWB-related miRNA and mRNA regulatory networks in P. fortunei. (A) Scatter plot of differentially expressed genes. (B) Enriched KEGG pathways among differentially expressed genes. (C) Volcano plot of miRNAs in WPF/WPFI. (D) PaWB-related miRNA and mRNA regulatory networks.

3.2 Phytoplasma Infection-Responsive miRNAs in P. fortunei

A total of 161,080,220 clean reads were obtained for six miRNA libraries (Table S3). Clean and unique miRNA reads were mapped to the P. fortunei genome sequence. The number of miRNAs in the WPF and WPFI libraries was very similar, indicating that the phytoplasma infection had little effect on the classification of miRNAs in P. fortunei. The results of the correlation analysis of the six samples (Fig. S3a) reflected the high repeatability of the miRNA sequencing results, making them suitable for the downstream analysis.

According to the examination of the miRNAs in WPF and WPFI, 88 differentially expressed miRNAs were identified (38 up-regulated and 50 down-regulated) (Fig. 1C). On the basis of a KEGG analysis, the target genes were assigned to 34 pathways (Table S4), among which plant hormone signal transduction (ko04075) was the most common, followed by plant–pathogen interaction (ko04626), circadian rhythm–plant (ko04712), and ubiquitin mediated proteolysis (ko04120) (Fig. S3b).

3.3 Construction of PaWB-Related miRNA and mRNA Regulatory Networks

A miRNA–mRNA regulatory network consisting of 254 miRNAs and 1517 mRNAs was established. The regulatory relationships of four miRNAs and eight mRNAs related to the infection of P. fortunei by phytoplasma are presented in Fig. 1D. The results of an association analysis showed that PfmiR396 was involved in the P. fortunei response to phytoplasma and affected PfGRF expression. The PfGRF gene family, which consists of plant-specific TF genes, has been found in A. thaliana, rice, and soybean. Notably, it influences plant resistance as well as the regulation of plant morphology, growth, and development [21,22,24]. Therefore, the PfGRF gene family was selected for subsequent analysis, which aimed to provide fundamental insights for the genetic improvement of Paulownia species and the breeding of new varieties with enhanced stress resistance.

3.4 Identification of PfGRF Gene Family Members and Promoter Cis-Acting Elements

16 PfGRF family members were identified and designated as PfGRF1 to PfGRF16 based on their chromosomal locations in P. fortunei (Fig. 2A,B, Table 1). All of these genes contained a QLQ domain and a WRC domain (Fig. 2A). The number of amino acids encoded by these PfGRF genes ranged from 148 to 620, with molecular weights ranging from 16.9 to 67.2 kDa. Their pI values ranged from 6.11 to 10.14. Among the proteins encoded by these 16 PfGRF genes, 13 were basic proteins (pI > 7), whereas three were acidic proteins (pI < 7). The subcellular localization results indicated that the PfGRF proteins were localized mainly in chloroplasts (Table 1).

Chromosomal maps were constructed to visualize the locations of PfGRF genes on each chromosome. Sixteen PfGRF genes were mapped to 10 chromosomes (Fig. 2B), while the remaining three PfGRF genes were mapped to unassembled scaffolds. The PfGRF genes were unevenly distributed among the chromosomes, with eight on each of chromosomes 7 and 18, two on each of chromosomes 2, 9, and 16, and only one on each of chromosomes 4, 15, and 19. There were some apparent regional enrichment in the distribution of PfGRF genes on chromosomes. For example, eight PfGRF genes were located at the end of chromosome 7 and the beginning of chromosome 18.

To examine PfGRF transcriptional regulation, the cis-acting elements in PfGRF promoters were identified. The PfGRF promoter regions were revealed to contain light-responsive, phytohormone-responsive, and stress-responsive cis-acting elements as well as cis-acting elements involved in the regulation of plant growth processes (Fig. 2C). Specifically, a lot of the PfGRF promoters contained methyl jasmonate or salicylic acid-responsive cis-acting elements, which may be involved in biotic stress responses. These findings suggest that PfGRF genes may participate in the response to phytoplasma infection leading to PaWB.

Figure 2: Distribution of PfGRF genes. (A) Multiple sequence alignment of PfGRF genes. QLQ: QLQ domain binding sites. WRC: WRC domain binding sites. *: positions of complete conservation. (B) Chromosomal distribution of PfGRF genes. (C) Analysis of PfGRF promoter cis-acting elements.

3.5 Phylogenetic and Collinearity Analyses of PfGRF Genes

A phylogenetic tree was constructed for the GRF genes in P. fortunei, A. thaliana, and O. sativa. The genes were clustered into six groups (Fig. 3A). The duplicated gene pairs PfGRF5/PfGRF11 were clustered in group I and were evolutionarily closest to AtGRF1 and AtGRF2, implying that their functions may be similar to those of AtGRF1 and AtGRF2. In addition, PfGRF1 was clustered with AtGRF5, whereas PfGRF3 was clustered with AtGRF6, suggesting that they may have similar functions. The duplicated gene pair PfGRF1/PfGRF3 was clustered with AtGRF6 in group II, implying its functions may be similar to those of AtGRF6. Both PfGRF8 and PfGRF12 were clustered with OsGRF3, OsGRF4, and OsGRF5 in group III. Group IV contained PfGRF6, PfGRF9, and PfGRF13, which were distantly related to the GRF genes in O. sativa.

Collinearity analyses of P. fortunei, A. thaliana, and O. sativa were performed. Fifteen duplicated gene pairs involving 10 PfGRF genes and A. thaliana GRF genes were detected (Fig. 3C) as well as 9 duplicated gene pairs involving 7 PfGRF genes and O. sativa GRF genes (Fig. 3B). The results of collinearity analyses indicated that PfGRF genes were evolutionarily closer to AtGRF genes than to OsGRF genes.

Figure 3: Phylogenetic and synteny analyses of PfGRF genes. (A) Phylogenetic analysis of the PfGRF proteins in P. fortunei, A. thaliana, and O. sativa. A neighbor-joining tree was constructed using MEGA-X, with 1000 bootstrap replicates; P. fortunei, A. thaliana, and O. sativa are differentiated by colors. (B) Synteny analysis of P. fortunei. The red line indicates the collinear gene pair in the P. fortunei genome. (C) Collinearity analysis of GRF genes in P. fortunei, A. thaliana, and O. sativa. Collinear gene pairs are indicated by blue lines.

3.6 Conserved Motifs and Structure of PfGRF Genes

To further characterize the PfGRF gene family, a conserved motif analysis was performed, which identified 10 distinct motifs in the 16 PfGRF genes, among which motifs 1 and 2 were present in all family members (Fig. 4). Motifs 4–7 were specifically detected in PfGRF1/2/3/4/8/14/15/16. PfGRF genes in the same branch of the phylogenetic tree had a similar motif composition. In contrast, PfGRF genes in different branches differed regarding their motifs. This suggests that the differences in the conserved motifs may be a key factor associated with the functional diversity among PfGRFs.

The gene structure analysis indicated that 12 of the 16 PfGRF genes lacked untranslated regions (Fig. 4A). Introns can affect gene stability, with genes containing many introns forming variable spliceosomes during transcription. The number of introns in the PfGRF genes ranged from 0 to 5 (Fig. 4A). More specifically, PfGRF2, PfGRF3, and PfGRF16 had no introns, whereas PfGRF11 had five introns, suggesting it may be unstable during transcription. Both PfGRF2 and PfGRF16 lacked untranslated regions, but had a number of introns, indicating these genes may also be unstable during transcription.

Figure 4: Phylogenetic, conserved motif, and structural analyses of PfGRF genes. (A) Phylogenetic tree (left) and genetic structure (right). (B) Motif logos.

3.7 Effects of Salinity and Drought on PfGRF Expression

To further explore the roles of PfGRF genes in abiotic stress responses, expression levels of PfGRFs under salt and drought treatments were analyzed using transcriptome data (Fig. 5A,B). After the salt treatment, PfGRF1/6/7/8/9/10/13 expression levels were up-regulated (compared with the corresponding control level), suggesting that these genes are responsive to salinity stress.

Following the drought treatment, the PfGRF1/4/5/6/8/9/10 expression level were down-regulated, which were in contrast to the significantly up-regulated expression of PfGRF11/12/13/15, indicating that these genes respond differentially to drought. A transcriptome sequencing analysis of samples exposed to salt and drought conditions detected differences in PfGRF expression levels, suggesting that PfGRF TFs may have different functions during responses to drought and salinity.

Figure 5: Expression analysis of PfGRFs under biotic and abiotic stress. (A) Heatmap of PfGRFs genes expression in response to drought. (B) Heatmap of PfGRFs genes expression in response to salt. (C) Heatmap of PfGRFs genes expression in response to phytoplasma.

3.8 Analysis of PfGRF Expression in Response to Phytoplasma Infections

Considering the diversity and similarity of the PfGRF genes revealed by the analyses of gene structures and evolutionary relationships, the transcriptome data were used to examine PfGRF expression levels in P. fortunei infected with phytoplasma. Differentially expressed PfGRF genes in response to PaWB were identified by analyzing RNA-seq data from WPF and WPFI plants. The results showed that the expression of all 16 PfGRF genes was affected by the phytoplasma presence (Fig. 5C). PfGRF1, PfGRF3, PfGRF4, PfGRF5, PfGRF6, PfGRF11, and PfGRF16 expression levels were significantly up-regulated, which was in contrast to the significantly down-regulated expression of PfGRF2, PfGRF7, PfGRF8, PfGRF9, PfGRF10, PfGRF12, PfGRF13, PfGRF14, and PfGRF15. Notably, PfGRF14 expression levels were down-regulated by 4.93 times, respectively, in response to the phytoplasma infection.

According to the expression profile analysis, PfGRF1/3/4/5/6/11/16 act as positive regulators of PaWB, whereas PfGRF2/7/8/9/10/12/13/14/15 function as negative regulators.

3.9 Identification of PfGRF14-Interacting Proteins

Previous studies have shown that GRFs usually interact with GIFs to participate in regulating the size of the blades [28,29]. PfGIFa (Pfo04g014660) was confirmed to interact with PfGRF14 through both Yeast two hybrid (Y2H) and bimolecular fluorescence complementation (BiFC) assays (Fig. 6). It is speculated that PfGRF14 may contribute to the development of the small-leaf symptom during PaWB infection.

Figure 6: Y2H (A) and BiFC (B) respectively detect the interaction between PfGRF14 and PfGIFa. All panels share the same scale bar of 100 μm.

4 Discussion

The value of miRNA–mRNA regulatory networks is reflected by their utility for thoroughly analyzing RNA transcription data, extending the study of regulatory mechanisms to the network level binding model, accurately and comprehensively revealing the differential expression patterns of RNA related to the occurrence of PaWB, and elucidating PaWB development. In this study, gene expression profiles in P. fortunei were analyzed and a miRNA–mRNA regulatory network was constructed, with the generated data useful for research on the biological functions of RNA in P. fortunei, with implications for studies on the P. fortunei response to PaWB. According to earlier research, GRF TFs are associated with plant stress resistance [18,19,20]. Khatun et al. (2017) reported that tomato GRF TF family genes play a crucial role in plant responses to abiotic stress and hormones. Other researchers determined that the expression levels of some GRF genes are affected by cold stress, while also revealing the DELLA–GRF regulatory module [24]. Liu et al. (2009) found that transgenic plants overexpressing miRNA396 targeting GRF genes have a lower stomatal density and stronger tolerance to drought than wild-type control plants [28]. In our study, miR396 was upregulated in susceptible P. fortunei, and its target genes PfGRF14/15 were downregulated, suggesting an expression-level association with PaWB symptom formation; direct regulatory roles remain to be experimentally confirmed.

Because of a lack of mutants and methods for generating transgenic Paulownia plants, GRF functions and regulatory mechanisms are unclear. The study of GRF functions has relied on expression analyses and bioinformatics-based predictions of gene functions, mainly in woody plants. In A. thaliana, AtGRF1, AtGRF2, AtGRF3, and AtGRF5, which affect leaf size, can interact with the GIF-interacting factor that enhances the expression of functional [11]. Interestingly, PfGRF14/15 are closely related to AtGRF3/4/5, implying that they may contribute to leaf size regulation in P. fortunei. The reduced expression of these genes in phytoplasma-infected plants correlates with leaf atrophy symptoms but does not constitute proof of regulatory function.

Previous studies showed that GRF genes are involved in signal transduction pathways related to responses to abiotic stresses, including salt, drought, and exogenous reagents [18,19,20]. In Brassica rapa, BrGRF5 expression decreases in response to salinity [15]. In maize, ZmGRF4 and ZmGRF13 expression is significantly induced by saline and drought conditions, suggesting that ZmGRF genes may play a key role in maize responses to these abiotic stresses [15]. In cassava, MeGRF4 is responsive to low temperatures and salt stress [29]. These expression changes suggest potential involvement in stress adaptation, yet functional roles require further validation through genetic manipulation.

In conclusion, the findings of this study have enriched the available information regarding the GRF gene family and elucidated the response of GRF genes to abiotic and biotic stresses. Furthermore, the study data provide a foundation for future studies conducted to comprehensively clarify the functions of GRF genes in P.fortunei, identify disease resistance genes in Paulownia species, and breed new PaWB-resistant varieties.