Revisiting the Role of Transcription Factors in Coordinating the Defense Response Against Citrus Bark Cracking Viroid Infection in Commercial Hop (Humulus Lupulus L.)

Transcription factors (TFs) play a major role in controlling gene expression by intricately regulating diverse biological processes such as growth and development, the response to external stimuli and the activation of defense responses. The systematic identification and classification of TF genes are essential to gain insight into their evolutionary history, biological roles, and regulatory networks. In this study, we performed a global mining and characterization of hop TFs and their involvement in Citrus bark cracking viroid CBCVd infection by employing a digital gene expression analysis. Our systematic analysis resulted in the identification of a total of 3,818 putative hop TFs that were classified into 99 families based on their conserved domains. A phylogenetic analysis classified the hop TFs into several subgroups based on a phylogenetic comparison with reference TF proteins from Arabidopsis thaliana providing glimpses of their evolutionary history. Members of the same subfamily and subgroup shared conserved motif compositions. The putative functions of the CBCVd-responsive hop TFs were predicted using their orthologous counterparts in A. thaliana. The analysis of the expression profiling of the CBCVd-responsive hop TFs revealed a massive differential modulation, and the expression of the selected TFs was validated using qRT-PCR. Together, the comprehensive integrated analysis in this study provides better insights into the TF regulatory networks associated with CBCVd infections in the hop, and also offers candidate TF genes for improving the resistance in hop against viroids.


Introduction
Plants confront numerous biotic stresses in their natural habitat through a compendium of well-tuned surveillance systems including a sophisticated genetic defense mechanism. The successful defense response in plants relies on the timely and accurate detection of the invading pathogen and the subsequent induction of the components of the defense response pathways to ward off the pathogens. The transcriptional reprogramming is one of the key events in a plant's defense response and involves a myriad of complex molecular, biochemical, and physiological changes, occurring in a highly synchronized manner largely governed by transcription factors (TFs) [1]. The lineage-specific TFs constitute the repertoire of master regulators, which mediate the transcriptional regulations by binding a specific DNA sequence in response to an environmental stimulus and developmental cues in plants [1]. TFs are central to the modulation of gene expression and play distinctive roles in diverse biological processes such as developmental control, and the activation of signaling cascades in defense and stress responses, by intricately regulating the transcription of downstream targets with the interaction of co-regulators known as transcriptional regulators (TRs) [2]. Mounting studies suggested that several families of TFs play an important role in the regulation of the expression of the plant defense transcriptome [2,3]. For example, WRKY TFs, regarded as the core of the plant immune system, play an indispensable roles in pathogen defense and phytohormone signaling [4]. Members of MYB/bHLH TFs or their complexes regulate distinct cellular processes such as responses to biotic stress, hormone signaling and cell death [5,6]. Similarly, bZIP TFs play critical roles in numerous cellular functions, including the regulation of stress-responsive pathways and the pathogen defense [7]. A recent explosion in the field of genomics and transcriptomics has enabled the family-wide identification and characterization of a large repertoire of TFs involved in a diverse array of biological processes in plants, such as the MYB, ADP-ribosylation factor, WRKY [8][9][10]; NAC [11]; bHLH [12]; Dof [13] and bZIP [14]. These studies have opened new avenues for exploring and exploiting TFs, to improve disease resistance in plants against pathogens. Viroids are the smallest (246-401 nt) autonomous plant pathogens, with the lowest biological complexity currently known. Their genome is approximately tenfold smaller than viruses and consists of single-stranded, circular, highly structured RNA that lacks a protein-coding capacity [15]. They are classified into two taxonomic families based on their biological, structural, and biochemical properties; the Pospiviroidae and the Avsunviroidae, whose members replicate in the nucleus and chloroplast respectively [16]. Viroids are cosmopolitan in their distribution and are etiological agents of diseases affecting a broad range of plants, including herbaceous, agronomic, ornamental and tree crops. In hosts, they cause minor to severe symptoms leading to plant death, depending on the host genome complexity and the aggressiveness of the viroid species [15]. In addition to symptomatic infections, viroids also cause latent infections without any noticeable symptoms. The simplicity of the viroid genome and their high rates of per-base in vivo mutation compared to other nucleic acid-based pathogens makes them an intriguing model for studying the RNA structure-functional relationship [15].
Hop (Humulus lupulus), a member of the Cannabaceae family, is a dioecious, climbing, herbaceous, perennial plant that is native to Europe, Asia, and North America [17]. Hop plants are commercially grown for their cones (strobiles) that are rich sources of secondary metabolites such as essential oils, bitter acids, and prenylated flavonoids. These compounds serve as essential raw material for the beer production industry, providing distinctive bitterness, flavor and aroma. In addition, the hop plant possesses several medicinal and soporific effects and are used as edible delicacies. Diseases caused by viroids impose significant global challenges for the commercial production of hop [18]. This includes Hop stunt viroid (HSVd) Apple fruit crinkle viroid (AFCVd), Hop latent viroid (HLVd), and Citrus bark cracking viroid (CBCVd) [19,20], of which the CBCVd infection is the most aggressive. Though several non-hosts are identified for viroids, no durable natural resistant sources have been discovered in hop species. In hop, macroscopic symptoms of a viroid infection include leaf distortion and mottling, stunting, chlorotic or necrotic spots, epinasty, vein discoloration and clearing, malformation of flowers and cones, scaling and cracking of the bark, and death of the plant in the case of CBCVd infections. The current management practices for viroid disease include detection and eradication, as well as agronomic maneuvers.
Viroid pathogenesis is extremely complex and is influenced by both the host and viroid genomes, resulting in mild or asymptomatic to severe symptoms [17]. Our recent transcriptome analysis of the CBCVd-hop interaction has revealed a massive transcriptional reprogramming of the genes involved in the immune response, sugar metabolism, photosynthesis, and phytohormone signaling pathways [21]. Nevertheless, transcriptional reprogramming is governed by TFs, and in this context it is crucial to delineate the pivotal role played by TFs in finely coordinating diverse regulatory networks. The involvement of hop TFs in the viroid pathogenesis is still enigmatic. In this study, we systematically identified hop TFs, categorized them in various TF families according to their conserved domains and examined their phylogenetic relationships to those of the model plant A. thaliana. Furthermore, we analyzed the expression signatures of a subset of hop TF families during CBCVd infection using the digital gene expression and a qRT-PCR assay. To the best of our knowledge, this study represents pioneering attempt to comprehensively identify and classify CBCVd-responsive hop TFs and provides valuable information about their involvement in hop-viroid pathogenesis.

Identification of Hop Transcription Factors
The previously published hop genome and CBCVd-responsive transcriptome profile generated at 412 days post-inoculation (dpi) [21,22] were used for the global mining of hop TFs. The comprehensive hop TFs database was constructed based on a homology search of the available transcriptome datasets of hop against the plant TF database (https://plantgrn.noble.org/PlantTFcat/) [23], which resulted in the classification of hop transcripts into various TF families. To validate this classification, the Hidden Markov Model (HMM) profile for each of the identified TF families was downloaded from the Pfam protein family database (http://pfam.sanger.ac.uk/), and a HMMER search (http: //hmmer.janelia.org/) was performed to identify the presence of conserved family domains in identified TF sequences. To filter out false positives, we cross-checked all TF sequences using a batch CD-search in the NCBI conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) and SMART database (smart.embl-heidelberg.de). All non-redundant sequences encoding complete TF domains were considered to be putative TF genes. Subsequently, a nucleotide-protein (blastx) (https://blast.ncbi.nlm.nih.gov/Blast.cgi) search was performed using the nucleotide sequences of the identified hop TFs to examine their homology with related sequences in the NCBI protein database.

Protein Characterization, Conserved Motifs and Gene Ontology Identification
The ProtParam tool available on the Expert Protein Analysis System (ExPASy) proteomics server (http://web.expasy.org/protparam/) was used to compute the isoelectric points (pI) and molecular weight of the hop TFs. The publicly available web server CELLO v 2.5 (http://cello.life.nctu.edu.tw/) [24] was employed for the prediction of the localization of hop TFs.
The Multiple Expectation-maximization for Motif Elicitation (MEME) programme (http://memesuite.org/tools/meme) was used to statistically identify the additional conserved motifs of the hop TFs with the following parameters: maximum number of motifs: 10, and motif length set to 6-100 and motif sites to 2-120; The distribution of one single motif was "any number of repetitions", and the other search parameter was "search given strand only".
To gain insights into the molecular function, biological process and cellular component of the identified hop TFs, a Gene ontology (GO) annotation was performed using the Blast2GO command line tool v 4.1 [25]. The TF nucleotide sequences were BLASTX searched against the non-redundant protein database of plants in NCBI under the default parameters. The GO terms associated to hit sequences were obtained by mapping and InterPro annotation.

Phylogenetic Affinity of Hop TFs
To deduce the phylogenetic relationships of hop TFs, the amino acid sequences were aligned along with the reference TF protein sequences of A. thaliana downloaded from the plant transcription factor database (http://planttfdb.cbi.pku.edu.cn/) using the MUSCLE program [26] with default parameters. Based on the alignment, unrooted phylogenetic trees were constructed for each of the TF families using MEGA v 7.0 software [27]. The criterion used for the tree construction were the Neighbor-Joining method, the JTT model, and complete deletion. The confidence level of the monophyletic groups was assessed using a bootstrap analysis based on 1000 re-samplings. The tree was edited with the Figtree tool (http://tree.bio.ed.ac.uk/software/figtree/).

The CBCVd-Responsive Expression Datasets Availability and Analyses of Differentially Regulated TFs
To gain an insight into the digital expression profiles of hop TFs during CBCVd infection, we have used our recently published CBCVd-responsive transcriptome [21]. Briefly, we mapped the clean reads to the unigenes dataset and normalized to the number of fragments per kilobase of exon per million mapped fragments (FPKM) by expectation maximization (RSEM) protocol using the Trinity software package [28]. The obtained count value was exported to the DESeq2 R package [29] to determine differentially expressed gene transcripts (DEG) using the Benjamini and Hochberg approach [30] for controlling the false discovery rate (FDR). The expression of TF encoding unigenes with an FDR adjusted p-value ≤ 0.05 and at least a two-fold change (≥2 or ≤2) was considered to be significantly differentially expressed between the CBCVd-infected and control libraries. The TF gene expression was visualized using a heatmap with clustering drawn using the online Clustvis tool [31].

Identification of Orthologous Genes and Protein Interactions
The identification of orthologous genes in well-characterized model plants helps in predicting the putative functions of newly identified genes in non-model plants such as hop. Therefore, the A. thaliana genes that are closest (orthologous) to the hop TFs were determined using an eggnog-mapper [32] based on the eggNOG 4.5 orthology data [33].
TFs often interact among one another during transcriptional reprogramming. To dissect the plausible interactions of hop TF proteins, a protein-protein interaction (PPI) analysis was performed using the STRING v 11 (http://string-db.org) [34] database in the COG (Clusters of Orthologous Group) mode, using the best-assigned COGs of A. thaliana. The interactions with a confidence score of ≥ 0.7 and based on co-expression and experiment conditions were used to construct the PPI network. The top 10 interactions were identified using the cytoHubba package [35], and the interaction network was visualized using the Cytoscape v.3.7.0 [36].

RNA Extraction and qRT-PCR Validation of Selected Hop TFs
The total RNA from the hop CBCVd-infected and control leaf samples was isolated by Concert™ Plant RNA Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol. The RNA concentration and quality were assessed using a Nanodrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Waltham, MA, USA). The total RNA was treated with the TURBO DNA-free™ Kit (Invitrogen, Carlsbad, CA, USA) to remove genomic DNA traces, and first strand cDNA was synthesized using the Superscript III Kit (Invitrogen, Carlsbad, CA, USA) from 1 µg of total RNA, according to the manufacturer's protocol. The obtained cDNA was diluted 10-fold and used as the template for the qRT-PCR analysis.
Eleven pairs of specific primers were used to study the expression of the selected hop TFs in the CBCVd infected leaf tissues of the hop plants. These TFs were carefully selected based on their reported defense response roles in various stages of the plant-pathogen interaction (Table S1). TF gene-specific primers were designed using the PerlPrimer software [37], and the specificity of the primer pairs was assessed in silico by performing an NCBI primer blast [38]. The qRT-PCR was performed in a CFX Connect™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) using the TopBio SYBR master mix (TopBio, Prague, Czech Republic). The reactions were performed in 15 µL comprising of 7.5 µL of SYBR Green, 0.75 µL (10 µM) of forward and reverse primers, 3 µL of RNase free-water and 3 µL of cDNA. The PCR cycling regime was at 95 • C for 3 min, followed by 40 cycles of 95 • C for 15 s, 58 • C for 30 s and 72 • C for 30 s. The specificity of each primer pair was assessed using a melt curve at the end of the reaction. The threshold cycles (Ct) of each TF target were averaged for triplicate reactions, and the relative transcriptional changes in the gene expression levels (fold-change) were calculated through the comparative Ct (2 -∆∆Ct ) method [39] using DRH1 (DEAD-box ATPase-RNA-helicase) [40] as the reference gene. Three biological replicates were used for each sample.

Identification and Classification of Hop TFs
With the aim of defining the hop TF gene families, we performed a global examination of the hop transcriptome using the plant TF database (https://plantgrn.noble.org/PlantTFcat/), which identified 3,818 non-redundant TFs that were classified into 99 TF families based on their conserved domains ( Table 1). The blast analysis showed that most hop TFs shared the highest homology to Morus notabilis, with the similarity ranging from 69% to 99% (Table S2). A large proportion of the TFs identified in this study belonged to the well-characterized TF families (e.g., WRKY, AP2, bZIP, MYB, bHLH, and NAC). Since it is impossible to describe the exhaustive list of TFs and their plausible roles in a single communication, we have in this study adopted the selective method to characterize the role of the major TF families, viz. WRKY, AP2, bZIP, MYB, bHLH, and NAC, that play indispensable roles in the plant immune response [3,41,42]. Table 1. Details of the hop transcription factor (TFs) families identified in this study. The TFs were identified using a combination of the plant TF database (https://plantgrn.noble.org/PlantTFcat/) and HMMER searches using the assembled transcriptome sequences of hop plants infected with Citrus bark cracking viroid (CBCVd).

Protein Characterization, Conserved Motifs and Gene Ontology Identification
The ExPASy analysis of the CBCVd-responsive hop TFs showed a range of variation in their molecular weight (7.2-96.3 kDa) and isoelectric point (4.5-9.89) ( Table S2). The subcellular localization analysis of the hop TF proteins was predicted, and most of them were found to be localized in the nucleus (445, 95.91%), followed by the chloroplasts (11, 2.37%), cytoplasm (5, 1.08%) and mitochondria (3, 0.65%) ( Table S3).
The conserved and novel additional motifs in the hop TF protein sequence were identified using the online MEME analysis tool. The analysis revealed that all of the TFs possessed their characteristic conserved core family domains. Most TF members that were clustered in the same subfamilies displayed similar motif characteristic, suggesting that the protein architecture is conserved within a subfamily and might reflect functional similarities among the members ( Figure S1). However, there were also some differences observed in the members among the subfamilies. Few subgroups had common motifs for all of the members, whereas other subgroups (e.g., AP2, WRKY TFs) possessed special motifs which could be attributed to their functional diversity ( Figure S1).
The functional annotation of the 466 hop TFs revealed their association with a range of biological processes, molecular functions, and cellular component categories ( Figure 1). In the biological process category, a predominant portion of the TFs was found to be involved in the cellular process, metabolic process, and in biological regulation (386) (Figure 1). The most highly enriched molecular function categories were related to binding (423), transcriptional regulator activity (261), and catalytic activity (6) (Figure 1). The cellular component category included terms related to the cell (296), organelle (84) and membrane (7).

Phylogenetic Affinity of Hop TFs
To explore the phylogenetic relationships between the hop TFs and to clarify the evolutionary relations within the family, an unrooted phylogenetic tree was constructed, along with A. thaliana reference TF proteins. As is evident from Figure 2

Phylogenetic Affinity of Hop TFs
To explore the phylogenetic relationships between the hop TFs and to clarify the evolutionary relations within the family, an unrooted phylogenetic tree was constructed, along with A. thaliana reference TF proteins. As is evident from Figure 2, the phylogram distributed the hop TFs into different clades based on their evolutionary relationships. Sakuma et al. (2002) [43] have described a total of 12 groups for the AP2/ERF family in A. thaliana. In our study, the 140 hop AP2/ERFs were grouped into 11 groups (A1, A2, A4, A5, A6, B1, B2, B3, B4, B5, and B6) with Group B3 forming the major clade with 37 members (Figure 2a). In A. thaliana, WRKY TFs have been classified into 3 main groups (Group I-III) (Eulgem et al. 2005) [44]. Here, the hop WRKY TFs were distributed into three major groups (Group I, Group IIa, Group IIb, Group IIc, Group IId, Group IIe, and Group III; Figure 2f), where Group I formed the major clade with 18 members. In the case of the A. thaliana bZIP family, 11 groups have been previously reported (Jakoby et al. 2002) [7]. Our phylogenetic tree categorized the bZIP TF proteins into 10 major groups (Figure 2c

Differential Modulation of Hop TFs
To understand the intricate differential modulation of hop TFs during CBCVd infection, we analyzed the TF digital gene expression levels in CBCVd-infected and control samples using FPKM. The result revealed the massive differential modulation of hop TFs, with FPKM values ranging from 0.14 to 894.55 in CBCVd infected libraries (Table S4)

Differential Modulation of Hop TFs
To understand the intricate differential modulation of hop TFs during CBCVd infection, we analyzed the TF digital gene expression levels in CBCVd-infected and control samples using FPKM. The result revealed the massive differential modulation of hop TFs, with FPKM values ranging from 0.14 to 894.55 in CBCVd infected libraries (Table S4). We constructed an expression profile heat map based on the expression data in different libraries of hop TFs (Figure 3). Most of the TFs (440, or 95%) depicted a significant differential expression pattern (p ≤ 0.05, logFC ≥ 2 or ≤ −2). Of the total 466 TFs,

Identification of Orthologous Genes and Protein Interactions
Information regarding orthologous genes in well-characterized model species can help in predicting the putative functions of newly identified genes. Here, we surveyed orthologous genes from the well-characterized model plant A. thaliana, which helped us in predicting the putative functions of the hop TFs (Table S5).

qRT-PCR Validation of Selected Hop TFs
Based on the orthologous relationships and putative known functions, we selected 11 different CBCVd-responsive hop TFs for the qRT-PCR validation. The analysis revealed a high degree of correlation between the RNA-Seq digital gene expression and the qPCR expression statistics ( Figure  6), implying the high quality and reliability of our transcriptome datasets. In addition, the differential TF proteins are known to interact with each other and regulate the transcription of downstream targets. Therefore, to predict the PPI among the CBCVd-responsive hop TFs, a protein-protein interaction analysis was performed using the STRING database, using A. thaliana as a model system. The analysis revealed that most of the CBCVd-responsive hop TF proteins were involved in multiple interactions, with more than one CBCVd-responsive TF protein belonging to the same or different TF family ( Figure 5). For instance, NAC TFs were found to be actively interacting with MYB TFs; while AP2 TFs were largely involved in the interaction with WRKY and NAC TFs, suggesting a highly coordinated and programmed interaction of CBCVd-responsive hop TFs. Additionally, several TFs were predicted to be acting as hub genes (e.g., NAC83, MYB85, etc.) in the interaction network, reinforcing their crucial role in transcriptional reprogramming (Figure 5e,d) during CBCVd-infection in the hop. A detailed description of the functions of proteins represented in the interaction network is presented in Table S6.

qRT-PCR Validation of Selected Hop TFs
Based on the orthologous relationships and putative known functions, we selected 11 different CBCVd-responsive hop TFs for the qRT-PCR validation. The analysis revealed a high degree of correlation between the RNA-Seq digital gene expression and the qPCR expression statistics (Figure 6), implying the high quality and reliability of our transcriptome datasets. In addition, the differential expression profile of the selected hop TFs further strengthens their candidature in transcriptional reprogramming during CBCVd infection in the hop. For instance, the WRKY (WRKY68&40), bZIP (bZIP61) and NAC (NAC78) TFs depicted a significant up-regulation, which was consistent in both the RNA-Seq and qRT-PCR analysis. Similarly, the CUP-SHAPED COTYLEDON 3 (NAC), SHINE 2 (AP2) and WRKY7 TFs showed a similar trend of down-regulation in the transcriptome and qRT-PCR analysis.

Discussion
The rapid, large-scale and coordinated intracellular and tissue-wide reprogramming of gene expression at the transcriptional level constitutes a crucial step in mounting the immune response against invading pathogens in plants [46]. The fine-tuned, complex regulatory system consisting of TFs and co-regulatory proteins functions as important components in the tight regulation of the transcriptional reprogramming of genes associated with the plant defense response. The growing body of research over the past years has identified and characterized the intriguing roles of several TF gene families involved in the plant defense response. However, a larger proportion of these studies have been centered on model plant species such as A. thaliana [6,47], rice [48] and maize [49], while detailed systematic information on major TF families in a genomic context and their response

Discussion
The rapid, large-scale and coordinated intracellular and tissue-wide reprogramming of gene expression at the transcriptional level constitutes a crucial step in mounting the immune response against invading pathogens in plants [46]. The fine-tuned, complex regulatory system consisting of TFs and co-regulatory proteins functions as important components in the tight regulation of the transcriptional reprogramming of genes associated with the plant defense response. The growing body of research over the past years has identified and characterized the intriguing roles of several TF gene families involved in the plant defense response. However, a larger proportion of these studies have been centered on model plant species such as A. thaliana [6,47], rice [48] and maize [49], while detailed systematic information on major TF families in a genomic context and their response against viroid pathogenesis is still enigmatic in the hop. In this study, the CBCVd-responsive transcriptome data generated in our laboratory, combined with in silico approaches, have enabled us to conduct a systematic characterization of TFs and their regulatory roles in defense responses against CBCVd infection in the hop. The combined approach enabled the prediction of 3818 TFs, along with their sequence features and the putative phylogenies of the largest families in the hop, which is significantly higher than the number of TFs (ca. 2000) identified in A. thaliana [47], suggesting that the high rate and their propensity for parallel expansion patterns is a consequence of an adaptive response against selection pressure [50]. The comparison of the TF gene families of hop with other plant species has shown considerable variations in the count, which could be attributed to genome complexity, the succession of genomic rearrangements and the shaping of the gene regulatory network during the evolution of organism complexity [51].
Moreover, the phylogenetic comparative method analysis was employed to further interrogate the evolutionary correlations and classification of TF proteins, which resulted in the clustering of the identified TF family of the hop into several subgroups according to their phylogenetic affinity toward all well-classified TF subgroups in A. thaliana, which suggested that the in silico approaches of the mining of TFs from the hop transcriptome dataset is comprehensive and reliable. In general, gene duplication and differentiation play a major role in the origin of new genes and gene functions [52]. Importantly, the DNA binding domains of TFs are extremely crucial for their biological functions. Accordingly, the classification of the TF gene family based on the phylogenetic trees may reflect the role of domains in the gene. Though the roles of most hop TFs remain to be elucidated, it is likely that members of a given group/subgroup share common evolutionary origins and might have similar physiological functions. Intriguingly, we observed that few hop TFs (e.g., AP2 and bZIP family TFs) that belonged to the same subfamily were clustered in different clades in the phylogenetic tree (Figure 2a,b), which could be attributed to the occurrence of duplication and divergent events in hop TF genes. The previous observation of the clustering of some MYB TFs of Foxtail millet in different subclades [53] corroborated our analysis.
It is generally accepted that overrepresented amino acid motifs in TFs tend to represent functional regions that are evolutionarily conserved across or within specific lineages [54] with similar biological functions. In this context, we analyzed additional conserved motifs outside the family domain in the hop TFs using the online MEME program. The analysis revealed that the majority of the TFs were grouped into the same subfamilies and shared a similar motif composition ( Figure S1). The presence of highly conserved motifs among proteins of the same subfamily supports the notion that they are likely produced via gene expansion during plant evolution and that they are essential for the group-specific functions within subfamilies. However, a high divergence in the motif composition was also observed for TFs belonging to different groups, reflecting the complex protein architecture and diverse functions of hop TF proteins. Novel TFs are likely to emerge through the rearrangement into novel signature domains or through new combinations of existing domains [47]. Our results were in agreement with the previous reports of the identification of a considerable domain conservation within a TF subfamily but a divergence among the subfamilies [55][56][57]. In addition, the occurrence of some special motifs in a few subfamilies confirms the fact that domain-shuffling processes might have also have played a relatively significant role in the evolution of hop TFs. However, further studies are needed to confirm these observations.
Orthologues are defined as genes in different genomes that have been created by the splitting of taxonomic lineages and are likely to retain the same function [58,59]. It has been widely accepted that the functions of the newly identified genes could be putatively predicted according to their orthologous counterparts [60]. In this study, we predicted the functions of the identified hop TFs using the orthologous genes from the well-characterized model plant A. thaliana. A large proportion of the hop TFs was assigned to have putative functions in diverse biological processes such as hormonal signaling, developmental processes, and abiotic and biotic stress responses (Table S4).
In a given cell or organism, biological and physiological processes, including cell-to-cell interactions and metabolic and developmental control, are regulated by protein-protein interactions (PPI) between different TFs [61]. In particular, the expression of a diverse set of eukaryotic genes are regulated by relatively small numbers of TFs by their different combinations and under different stress conditions; the physical interactions among different combinations of TFs dictate and regulate the expression of specific genes [62]. Moreover, the analysis of interactions between transcription factors and other proteins also help in elucidating their role in different signaling cascades [63]. Several explicit evidences supported that the PPI network could be exploited for the manipulation of plants for disease resistance against pathogens [64,65]. Therefore, in this study, we have extensively investigated the putative protein-protein interactions (PPI) of CBCVd-responsive TFs to gain an understanding of the regulatory networks that modulate CBCVd-defense responses in the hop. In this study, the interaction network explicitly projected the complex interaction web of CBCVd-responsive hop TFs with diverse proteins involved in defense response, hormone signaling, and transcription. TFs act as sites of signal convergence in concert with other context-specific TFs and transcriptional co-regulators to determine the activation or repression of defense response pathways which ultimately establish sensory transcription regulatory networks that are required for plant immunity [3]. The PPI network explicitly revealed the interfamily interactions of hop CBCVd-responsive TFs (e.g., NAC TFs interacting with MYB and bHLH TFs ( Figure 5e); AP2 TFs interacting with WRKY TFs (Figure 5a)), reinforcing a complex, interconnected and finely-tuned transcriptional reprogramming of hop TFs in activating downstream defense pathways against CBCVd infection. A growing body of evidence suggests that WRKY TFs respond to biotic stresses through the auto-regulation and cross-regulation of other WRKY TFs, while interacting with a diverse array of defense responsive genes, forming an interaction hub [61,66]. The presence of a highly interconnected hub of several WRKY TFs, along with their dominant interactions with the mitogen-activated protein kinase 3 (MAPK3) (a component of the kinase module that plays a pivotal role in the transduction of diverse extracellular stimuli, including biotic stress, in plants) and the Non-Expresser of Pathogenesis Related Gene 1 (NPR1), a receptor for the plant defense hormone salicylic acid (SA) (Figure 5f), provide strong evidence for the auto/cross-regulation phenomenon of WRKY TFs in our PPI network. In addition, the explicit interaction of WRKY TFs with NPR1 suggests the crucial role of a SA mediated defense response in CBCVd-infected hop plants. The involvement of Brassinosteroid pathway (BR) genes in the viroid interaction has been previously reported in PSTVd-infected tomato [67] and potato [68]. In agreement with the above reports, we observed the involvement of BR pathway genes in CBCVd infection of the hop, as evidenced by the profound interaction hub of BIM1 & 2 with the hop bHLH TFs (Figure 5b). Together, the PPI interactions presented here are crucial to an understanding of the transcriptional regulatory networks that operate in triggering defense responses in hop against CBCVd infection, and they provide glimpses of putatively activated defense pathways and candidate proteins interacting with them.
The successful activation of the plant defense signaling involves the coordinated regulation of multiple transcriptional regulators, as well as DNA-binding transcription factors and their regulatory proteins, that rapidly reprogram the transcription in the plant cell in response to biotic stress [46]. To gain a better understanding of the intricate differential modulation of CBCVd-responsive hop TFs, we conducted a digital gene expression profiling using the FPKM derived from the hop-CBCVd transcriptome dataset. Our results illustrated the massive transcriptional reprogramming of hop TFs during CBCVd infection, with a higher percentage of TFs depicting a significant up-regulation ( Figure 4). It is worth noting that the AP2, bZIP, and WRKY families were found to be predominantly participating in the transcriptional reprogramming in response to CBCVd infection in the hop, as evidenced by their pronounced differential modulation (Table S3). The abundance of AP2, bZIP and WRKY TF families in the transcriptional reprogramming during viroid interaction has also been recently reported in tomato plants infected with Potato spindle tuber viroid (PSTVd) [69]. The AP2, bZIP, and WRKY TFs represent the largest repertoire of TF families, and considerable evidence suggested that they play pivotal roles in the regulatory networks controlling plant development, hormone signal transduction and disease resistance, in response to several plant pathogens [41,44,70].
In plants, the mitogen-activated protein kinase (MAPK) cascade plays a central role in the signal transduction of PAMP-triggered immunity (PTI) and Effector-Triggered Immunity (ETI) pathways, to activate defense-related genes and TFs: additionally, WRKY TFs play indispensable roles in this pathway. Consistent with the previous studies in A. thaliana and tomato plants infected with PSTVd [71,72], we observed the massive up-regulation of WRKY22 TFs, suggesting their involvement in the activation of the MAPK mediated defense signaling pathway in hop upon CBCVd infection. The intricate transcriptional regulation of stress-responsive genes is crucial to the plant's defense response, and WRKY TFs act as the master regulators for several defense responsive genes [4]. In this context, we observed an extensive up-regulation of several WRKY TFs (WRKY40, 68, 72, 75 and 31) in the CBCVd-responsive transcriptome. Previous studies have shown the crucial role of these WRKY TFs in the basal defense in tomatoes and A. thaliana [73,74] or in the resistance to Ralstonia solanacearum in tobacco [75]; Xanthomonas oryzae in rice [76]; and Phytophthora sojae in soybeans [77]. Therefore, we strongly believe that the suite of WRKY TFs plays an indispensable role in triggering defense signaling cascades in response to CBCVd infection in the hop. In addition to the WRKY TFs, a strong up-regulation was observed for bZIP and NAC TFs families. The TGA7 (bZIP TF) transcription factors, which are implicated as regulators of pathogenesis-related (PR) genes in A. thaliana [78] and NAC78 (having an important role in defending against Alternaria alternata infecting Callery pears [79]) were found to be profoundly induced in our study.
The intricate network of phytohormone signaling pathways enables plants to activate appropriate and effective defense responses against pathogens by maintaining an appropriate balance between defense and growth [80]. Among other classes of phytohormones, salicylic acid (SA), jasmonic acid (JA), and ethylene (ET) are regarded as core immune phytohormones. The cross talk between phytohormone defense-signaling pathways are common and can be either mutually antagonistic or synergistic, resulting in negative or positive functional outcomes and thus providing the plant with a powerful regulatory potential. In line with the previous observation of the high up-regulation of genes involved in the SA and JA hormone signaling pathways in PSTVd infected tomato plants [69], we observed the significant up-regulation of TFs involved in these phytohormone signaling pathways in the CBCVd-infected hop. For instance, WRKY22 and 27 showed a significantly high expression, which is known to have dominant roles in the SA-JA cross talk during pathogen defense. The CaWRKY27 TF in pepper (Capsicum annum) positively regulates the stress resistance response to R. solanacearum infection through the modulation of SA, JA and ET-mediated signaling pathways in tobacco (Nicotiana tabacum) [81], whereas AtWRKY22 mediates the cross talk between SA and JA-dependent defense pathways against bacterial and fungal pathogens in A. thaliana. [82,83]. Similarly, ERF 109, which is shown to mediate the cross-talk between jasmonic acid and auxin biosynthesis in A. thaliana [84], was found to be up-regulated in our study.
It is well known that a TF can act as an activator or repressor in the transcriptional regulation of gene expression [3]. Consistent with this observation, we found that 84 TF genes (18%) were down-regulated in our study, suggesting their role as a repressor during CBCVd infection in the hop. It is well established that when a plant perceives a pathogen, it switches to produce compounds necessary for a defense response by strategically compromising the mode of normal growth and development through the massive transcriptional regulation of gene expression via TFs [85,86], which was also evident from the results of the present study. Intriguingly, the SHINE and CUP-SHAPED COTYLEDON 3 TFs, which are involved in the regulation of flower organs [87] and leaf development [88] respectively in A. thaliana, were down-regulated in our study, suggesting their plausible roles in leaf deformation and cone size reduction upon CBCVd-infection in the hop.
In conclusion, our results provide a global overview of the role of TFs in mounting the defense response against CBCVd-infection in the hop. However, a functional analysis based on experimental approaches is essential for the screening of suitable TFs and their involvement in the transcriptional network in the defense response against CBCVd infection in the hop, which could be further utilized for genetic engineering and breeding programmes to improve crop resistance against viroid pathogens.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4915/11/5/419/s1, Figure S1 Phylogenetic relationships and architecture of conserved protein motifs in hop TF genes. The phylogenetic tree was constructed based on the full-length sequences of hop TF proteins using MEGA 7 software. The motifs numbers 1-10, are displayed in different colored boxes and sequence information of the motifs are presented in the adjacent box. The sequence information for each motif is provided in the boxes. Table S1: List of qRT-PCR primers used for the validation of expression profiles of selected hop TFs.