Genome-Wide Transcriptomic Analysis Reveals Insights into the Response to Citrus bark cracking viroid (CBCVd) in Hop (Humulus lupulus L.)

Viroids are smallest known pathogen that consist of non-capsidated, single-stranded non-coding RNA replicons and they exploits host factors for their replication and propagation. The severe stunting disease caused by Citrus bark cracking viroid (CBCVd) is a serious threat, which spreads rapidly within hop gardens. In this study, we employed comprehensive transcriptome analyses to dissect host-viroid interactions and identify gene expression changes that are associated with disease development in hop. Our analysis revealed that CBCVd-infection resulted in the massive modulation of activity of over 2000 genes. Expression of genes associated with plant immune responses (protein kinase and mitogen-activated protein kinase), hypersensitive responses, phytohormone signaling pathways, photosynthesis, pigment metabolism, protein metabolism, sugar metabolism, and modification, and others were altered, which could be attributed to systemic symptom development upon CBCVd-infection in hop. In addition, genes encoding RNA-dependent RNA polymerase, pathogenesis-related protein, chitinase, as well as those related to basal defense responses were up-regulated. The expression levels of several genes identified from RNA sequencing analysis were confirmed by qRT-PCR. Our systematic comprehensive CBCVd-responsive transcriptome analysis provides a better understanding and insights into complex viroid-hop plant interaction. This information will assist further in the development of future measures for the prevention of CBCVd spread in hop fields.


Introduction
Viroids are small, single-stranded, circular, highly structured, non-protein-coding infectious RNAs with genomes ranging in size from 250 to 401 nt [1]. The viroid genome is approximately tenfold smaller than the smallest RNA virus and they solely depend on their own RNA, host factors, and host enzymatic machinery for autonomous replication and movement [2,3]. In addition, viroids have intrinsic potential for high rates of per-base in vivo mutation among all the nucleic modulation of miRNAs that are involved in several hormone pathways and transcriptional factors involved in the regulation of metabolism, growth and development [18]. This observation directed our research towards uncovering the transcriptional reprogramming associated with CBCVd pathogenesis in hop. In the present study, we employed a high throughput transcriptome sequencing approach to gain a comprehensive understanding of the global alteration in gene expression resulting from CBCVd-infection in hop, which would further facilitate the development of effective measurements against viroid diseases.

Preparation of Dimeric CBCVd Construct and Inoculation of the Hop Plants
The full-length monomeric cDNA of CBCVd was amplified from total RNAs isolated from CBCVd-infected leaves of the Slovenian hop cultivar "Celeia" by reverse transcription-polymerase chain reaction (RT-PCR) using specific primers (Table S1). The purified RT-PCR product was cloned into T-Vectors (Takara Bio Inc., Shiga, Japan) and it was confirmed by sequencing. The dimeric construct was generated by digesting SacI termini of monomeric positive cDNA strand of CBCVd and cloning into the SacI site of pBlueScript KS (+) vector (Stratagene, San Diego, CA, USA). For 35S::CBCVd construct, the dimeric cDNA of CBCVd was subcloned into the pLV07 binary vector via intermediate vector pLV68 harboring a 35S expression cassette as described previously [28]. The CBCVd cDNA was immobilized on microcarrier gold particles (1 µm) using a calcium-mediated precipitation protocol [29] and biolistically inoculated on three-month-old hop leaves (cv. Osvald's 72) grown in a container (10 cm height with at least three shoots-stage) under natural light condition. The individual hop plant was inoculated five times with 250 ng DNA per viroid species and was placed in darkness for 24 h after covering with plastic bags. Subsequently, the mock and CBCVd-inoculated plants were grown under natural condition and inspected visually for symptoms development.

RNA Isolation, Detection and Quantification of Genomic RNA of CBCVd
Systemically infected and mock-inoculated leaves were harvested at 120 and 412 days post inoculation (dpi) from the shoot apex after appearance and assessment of disease symptom in CBCVd-infected hop plants. The high-quality total RNAs were extracted using Concert™ Plant RNA Purification Reagent (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA), followed by RNA purification and DNA contamination removal using DNA-free TM DNA Removal kit (Ambion, Thermo Fisher Scientific, Waltham, MA, USA), according to the manufacturer's instruction. The concentration of total RNA samples was measured at 260 nm absorbance using NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), whereas the integrity of RNA samples used for cDNA libraries construction was confirmed by Agilent Technologies 2100 bioanalyzer (Agilent Technologies, Colorado Springs, CO, USA).
The hop plant samples were examined for the CBCVd-infection using RT-PCR and dot blot assay. RT-PCR amplification for CBCVd detection was performed using an One-Step RT-PCR Kit (Qiagen, Valencia, CA, USA). Three microliters of RNA sample (adjusted to a concentration of 0.1 µg µL −1 ) was added to master mixture buffer containing 4 µL 5× One-Step RT-PCR Buffer, 0.6 µL (20 µM) of each CBCVd-specific forward and reverse primer (Table S1), 0.8 µL dNTP mix (10 mM each dNTP), 0.8 µL of One-Step RT-PCR Enzyme Mix, and 10.2 µL RNase-free water. Reaction mixtures were incubated in a PE9700 thermal cycler (Bio-Rad, Richmond, CA, USA) using the following reaction condition: 50 • C for 30 min for the reverse transcription step, followed by a PCR amplification step with initial denaturation at 95 • C for 15 min and 40 cycles of cDNA melting at 94 • C for 15 s, primer annealing at 58 • C for 30 s, primer extension at 68 • C for 45 s, accompanied with final primer extension at 68 • C for 7 min. The amplification of PCR products was confirmed by agarose gel electrophoresis.
Real time quantitative PCR (qRT-PCR) was performed to determine (+) and (−) vd-sRNAs accumulation of CBCVd in hop plants by following two steps to overcome challenges that are associated with highly stable viroid RNA secondary structure [30]. In the first step, 1 µg of total RNA was mixed with 20 µM of either CVdCS_PS or CVdCS_MS primer (Table S1) to make volume 13.4 µL. The mixture was denatured at 95 • C for 3 min and immediately chilled in ice. Subsequently, components of One-Step RT-PCR Kit (Qiagen, USA) were added to make final volume of reaction mixture to 20 µL and reverse transcription was performed, as mentioned above. In the second step, the quantification was performed using 3 µL of 10× diluted RT samples, which was added to 20 µL reaction mixture containing 10 µL 2× SYBR Green Real-Time PCR Master Mix (Invitrogen, USA), 6 µL RNase-free water, 0.5 µL of each 10 µM forward (CVdQRT_F), and reverse (CVdQRT_R) primers. qRT-PCR cycles were performed on an IQ5 Multicolor Real-Time PCR Detection System (Bio-Rad, USA) with following temperature profile: 94 • C for 4 min, followed by 40 amplification cycles 94 • C for 20 s, 61 • C for 40 s, and 72 • C for 30 s. In addition, the melting curve program (heating rate of 0.1 • C per s and a continuous fluorescence measurement) was performed to examine the specificity of the amplified product. The 7SL RNA gene product was used as an internal reference gene for normalization of expression level [31]. All qRT-PCR experiments were performed in triplicate using the independent CBCVd-infected (CI) and mock-inoculated (MI) samples that are mentioned above. The quantification of vd-sRNAs of CBCVd and data analysis were performed in accordance with MIQE guidelines [32].
The dot blot analysis of CBCVd was performed with full-length CBCVd [α 32 P] UTP-labelled probes by following the previously described method [30].

Transcriptome Sequencing, Assembly and Differential Gene Expression Profiling
Total RNA from four CBCVd-infected (CI) and four mock-inoculated (MI) in equal quantities were used for RNA sequencing. Total RNA (5 µg) was pooled from each sample in same batch from three individual leaves and was used for cDNA synthesis using cDNA Synthesis System (Roche, Basel, Switzerland), as per the manufacturer's protocol. The cDNA samples were sheared via nebulization into small fragments and they were used for library construction using Illumina TruSeq Stranded Total RNA Sample Prep with Ribo-Zero (plant) kit (Illumina, San Diego, CA, USA). The quality of library was evaluated using the Agilent 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA). The resulting libraries were then paired-end sequenced (2 × 100 bp) on an Illumina Hiseq™ 2500 platform (Illumia) using myGenomics (Atlanta, GA, USA) sequencing services. The paired-end RNA-seq reads were subjected to the removal of adapters, filtering of empty and low-quality reads, trimming of reads with ambiguous nucleotides based on PHRED quality scores (Q-score) using default parameters of Trimmomatic v0.30 program [33]. The trimmed reads shorter than 40 bp were dropped to eliminate the sequencing artifacts and the quality of reads were evaluated using FastQC tool [34]. The high-quality reads were de novo assembled using Trinity software package version v2.4.0 [35] with default parameters settings (K mer = 25) and termed as unigenes. The evaluation of assembly was performed using Bowtie2 aligner by mapping the filtered reads against unigenes. All assembled unigene sequences were compared with the hop transcriptome database of HopBase genomic resources repository (http://hopbase.org/) using MEGABLAST at typical cut-off E-value of 1.0 × 10 −5 , with similarity level and alignment length more than 95% and 100 bp, respectively. In order to calculate the expression level of each unigenes, the clean reads for eight experimental samples were mapped 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 in-built scripts in the Trinity software package [36]. The obtained count value was exported to DESeq2 R package [37] for determining differentially expressed gene transcripts (DEG) using the Benjamini and Hochberg approach [38] for controlling the false discovery rate (FDR). The expression of unigenes with FDR adjusted p-value ≤0.05 and at least a two-fold change (≥2 or ≤−2) was considered significantly different between CI and MI-libraries. The gene expression values were imported from DESeq2-normalized FPKM data sets and matrix distance for expression heatmap was calculated using the Euclidean distance and complete-linkage methods. A heatmap was constructed using the R statistics package heatmap3 [39].

Functional Annotation and Gene Enrichment Analysis
The assembled unigenes were aligned against the NCBI non-redundant (nr) protein database of Viridiplantae using BLASTX with a significance cut-off E-value of 1.0 × 10 −5 . Blast homology searches and homology-based functional annotations were carried out using Blast2GO Command line tools (Version 1.4.1) (Biobam, Valencia, Spain) [40]. Gene Ontology (GO) terms comprising of three functional groups, such as biological processes, molecular functions, and cellular components were assigned to the unigenes by Blast2GO program. The bidirectional best hit (BBH) method was used for KEGG (The Kyoto Encyclopedia of Genes and Genomes) pathway assignment of the assembled sequences using the online KEGG Automatic Annotation Server (KAAS; http://www.genome.jp/kegg/kaas/) to gain an overview of the biological pathways. The DEGs associated GO terms were enriched with respect to the GO terms that are associated to non-differentially expressed genes. The Fisher statistical test was performed to find enrichment of functional categories with Bonferroni's correction (FDR ≤ 0.05) using the AgriGO toolkit [41]. The cut off p-value less than 0.05 was used as qualifying parameter for GO terms enrichment analyses, which were visualized using ReviGO [42]. Similarly, KEGG metabolic pathway annotation and enrichment of the DEGs were performed using hypergeometric test equivalent to Benjamini and Hochberg's correction method with 5% FDR. The pathway visualization of DEGs was performed using MapMan tool (Forschungszentrum Jülich, Germany) [43]. The Log 2 FC values of DEGs were assigned to functional categories (or bins) by Mercator (http://mapman.gabipd.org/web/guest/mercator). In the case of expression data for duplicated gene identifiers, the lower value of fold-change was used for the analysis to avoid an overestimation of the data. The logarithm values of gene expression values were used to construct the gene regulatory modules following the methods as described previously [44]. The ortholog group assignment between DEGs of hop and Arabidopsis thaliana were performed using OrthoMCL [45] and selected DEGs with A. thaliana ortholog were subsequently used for network construction using NetworkAnalyst [46].

Validation of Differentially Expressed Genes by qRT-PCR
To confirm the results of transcriptome data, 12 randomly selected candidate DEGs were subjected to quantitative real-time PCR (qRT-PCR) analysis using designed specific primers (Table S1). Aliquots of the total RNA (5 µg) extracted for sequencing from leaves of CI and MI hop plants was treated with DNase I (Ambion, USA) and reverse-transcribed to first strand cDNA using Superscript ® III First-strand cDNA Synthesis kit (Invitrogen, USA), according to the manufacturer's instructions. PCR amplification was performed in an IQ5 Real-Time PCR Detection System (Bio-Rad, CA, USA) using 20 µL of reaction mixture containing 10-fold diluted cDNA, 10 µL SYBR Green Real-Time PCR Master Mix (Invitrogen, USA), 10 µM of forward and reverse gene-specific primers (Table S1) under following conditions 95 • C for 5 min, followed by 40 cycles at 95 • C for 15 s, 58 • C for 30 s, and 72 • C for 30 s. The relative expression levels (fold-change) of the selected genes were calculated by the comparative C t (2 −∆∆Ct ) method [47] and hop glyceraldehyde 3-phosphatedehydrogenase (GAPDH) gene was used to normalize the amount of template cDNA in each reaction as the internal control [48]. The fold change of each gene was determined by three independent biological replicates and based on that standard deviation was calculated.

Biolistic Inoculation of Hop with cDNA of CBCVd and Infectivity Confirmation
The hop plants were biolistically inoculated with infectious dimeric construct of cDNA of CBCVd ( Figure 1) and examined for incidence of infection at pre-dormancy period (120 dpi) and after the appearance of typical symptom at post-dormancy period (412 dpi) following dot blot analysis, RT-PCR, and qRT-PCR ( Figure 1). The positive hybridization signal was detected with minus strand-specific probe in nine out of fifteen CBCVd-inoculated plants at the pre-dormancy period ( Figure 1B). The RT-PCR product of approximately 228 bp was observed in samples that were positive to CBCVd infection in dot blot analysis ( Figure 1C). The strand-specific RT-qPCR (ssRT-qPCR) confirmed the result of minus strand-specific dot-blot hybridization and suggested the higher accumulation of minus multimeric strand as compared with plus polarity ( Figure 1D). The nine CI hop plants were visually inspected for onset of typical symptoms such as stunted growth, leaf malformation, and bine cracking during the 14-months ( Figure S1). After the dormancy period (14-month), dot-blot hybridization and ssRT-qPCR displayed the similar trend of excess of minus strand over plus strand in all nine symptomatic plants, which was corroborated with our previous studies suggesting that CBCVd utilizes minus strand as replicative intermediates [30].

Illumina Sequencing, De Novo Assembly and Functional Annotation of Unigenes
To compare the transcript profiles of hop in response to CBCVd-infection, total RNA was isolated from leaves of the apex of CI and MI prior to dormancy and post-dormancy of each individual plants and were pooled together in an equimolar amount to yield a final sample for sequencing with four biological replicates. Illumina RNA-Seq deep-sequencing run generated over 33.75 and 40.02 million raw reads in CI and MI libraries, respectively ( Table 1). The removal of adapter sequences, shortest reads, ambiguous regions, and filtering out low-quality reads at high stringency using Trimmomatic software, resulted in retention of 24.39 and 36.68 million high-quality reads in CI and MI samples, respectively. The sequence abundance profiles of biological replicated samples showed high correlation (Table S2), indicating high quality reliable gene transcript data. De-novo assembly with "Reduce" option to reduce redundancy in assembled unigenes and mapping reads back to contigs produced 27,094 unigenes from MI and CI combined libraries with length ranging from 90 to 2590. The N50, N75 value of de-novo assembly and average length were computed as 469 bp, 324 bp and 433 bp, respectively. The average unigene size of hop was longer than the average length of unigenes those identified in previous studied in Camellia sinensis (355 bp) [49], Spartina alterniflora (386 bp) [50], and Eucalyptus grandis (247 bp) [51]. The average GC content of assembled unigenes was 41.50% with normal distribution and were comparable to White clover (38.90%) [52], Spinach (42.50%) [53], Sophora flavescens (39.30%) [54].

Identification and Functional Classification of Differentially Expressed Genes
To identify the differentially expressed genes (DEGs) that are associated with CBCVd-infection in hop, gene expression levels in CI and MI samples were calculated using FPKM (expected number of Fragments Per Kilobase of transcript sequence per Million base pairs sequenced) [55]. The result of mapping of all reads utilizing non-redundant set of hop transcripts illustrated that the number of reads corresponding to each transcript ranged from 0.57 to 4716.55 (FPKM) in CI libraries, suggesting the extensive change of expression levels of hop transcripts (Table S3). The comparative transcriptome abundance analysis revealed significant differential expression level of 2209 unigenes (DEGs, p ≤ 0.05, logFC ≥ 2 or ≤−2), of which 1358 unigenes were up-regulated and 851 unigenes were found to be down-regulated in CI plants, which suggested the significant impact of CBCVd infection on host gene expression in leaves (Table S3). Approximately, 96.10% of DEGs (1297 up-regulated and 826 down-regulated) were annotated using BLASTx procedure against the nr-protein database of NCBI (Table S3). The hierarchical cluster analysis of the 2209 DEGs based on their FPKM values using the cluster distance method associated with complete linkage illustrated the relationship and degree of responses of DEGs in leaves of hop against CBCVd-infection ( Figure 3). The candidate unigenes in clusters II, III, and IX were found to be up-regulated, whereas genes in clusters I and V were found to be significantly down-regulated in leaves of CI as compared to MI ( Figure 3). Conversely, DEGs that are associated with cluster VII and VIII showed an intermediate pattern of expression. The annotation based on GO terms categorized 2209 DEGs into three classes: biological processes (19 terms, 45.42%), cellular function (10 terms, 31.37%), and molecular function (7 terms, 23.20%), whereas, 570 DEGs were not classified (Figure 2). Many of the DEGs are involved in cellular process, and catalytic activities. Furthermore, the GO enrichment analysis in this study provided the overview of statistically significant and relevant GO terms, which were altered as a result of CBCVd-infection in hop. A total of 32 GO terms were screened at p-value of 0.05 and statistically significant enriched GO terms, namely biological process, cellular component and molecular function was plotted in two-dimensional scatterplot in format of semantic representations of similar terms in closest coordinates (Figure 4). Among the various biological process categories, the GO terms linked to "metabolism", "response to stimuli", "cellular homeostasis" "transport", etc. were significantly enriched (Figure 4). These pathways are associated with defense responses and are typically induced in plant-biotic interactions [56,57]. Similarly, in the molecular function category GO terms related to catalytic activity, such as "transferase activity, kinase activity" and binding activity, such as purine nucleoside binding, anion binding, nucleoside binding, adenyl nucleotide binding, and ATP binding were significantly enriched (Table 3), which plays a pivotal role in signal recognition and signal transduction during host-pathogen interaction. In addition, GO terms related to cellular components, such as thylakoid, plastid thylakoid membranes, and photosynthetic membranes were significantly enriched in the cellular component category.   The alignment of DEGs to COGs database for orthologous genes resulted in the classification of DEGs into 23 functional categories with 1987 functional annotations due to multiple COG functions of some DEGs ( Figure S3B).
The mapping of DEGs in KEGG pathway assigned them to 7 KEGG pathway ( Table 2). Several other disease-resistance pathways, including plant hormone signals, biosynthesis of secondary metabolites, amino acid, nitrogen, nucleotide, fatty acid, and sugar metabolism pathway were also found to be enriched (Table 2).
Furthermore, DEGs were subjected to MapMan analysis to gain an unbiased, systematic overview of important biological functions and their coordinated response to CBCVd-infection in hop ( Figure 5). Plants respond to pathogen attack by deploying an increased demand of energy and biosynthetic capacity, which is required for the priming of them for early defense [58]. As expected, genes that are involved in metabolism overview pathways, such as "cell wall", "lipids", and "secondary metabolic" pathways were consistently up-regulated. A cohort of genes associated with secondary metabolism pathway (phenlypropanoids, flavonoids lignin and glucosinolates) were highly expressed, which were the most significantly enriched pathways in the KEGG analysis. The down-regulation of genes that are involved in photosynthesis (e.g., PSI, PSII) and chlorophyll (tetrapyrrole) biosynthesis could be correlated to the appearance of disease symptoms, such as chlorosis on CBCVd-infected leaf tissues in hop. The expression of sucrose and starch biosynthesis genes was found to be diminished, which could be attributed to a reduction of photosynthesis activity. In addition, the hormone and signaling pathway related genes changed significantly. The gene regulatory network analysis revealed potential regulatory elements of CBCVd-disease induction in hop ( Figure S4). The identified hub genes might play crucial roles in regulating the innate immune response in hop against CBCVd-infection, but further experimental studies are required to validate the results.

Validation of Expression Pattern of Candidate DEGs by qRT-PCR
In order to validate the results that were obtained from Illumina sequencing, twelve CBCVd-responsive annotated unigenes were selected randomly (based on expression patterns) and were subjected to qRT-PCR analysis using specific designed primers (Table S1). The result showed the consistent pattern of fold-changes of transcripts in qRT-PCR and RNA-Seq experiments ( Figure 6). For example, both qRT-PCR and RNA-Seq analyses showed that genes encoding MLP-like protein 423, AP2-like ethylene-responsive transcription factor, and pathogenesis-related protein 1 were significantly higher and consistent in CI as compared to MI leaves of hop. The expression of genes encoding phosphomethylethanolamine N-methyltransferase and auxin efflux carrier component was found to be suppressed in qRT-PCR analysis, which was similar to the transcriptome data analysis. The expression pattern of two genes encoding S-formylglutathione hydrolase and Homeobox protein BEL1 were found not to be significantly up-regulated or down-regulated in CI leaves, which was in accordance with qRT-PCR expression patterns. Nonetheless, the expression patterns of genes by qRT-PCR were consistent with the DEGs analysis, indicating the reliable result.

Discussion
CBCVd has evolved a broad host range across host genotypes such as citrus, cucumber, tomato and eggplant etc. [59]. Recently, it has been discovered in hop, where it causes the most aggressive symptoms after first dormancy and complete dieback in 3-5 years [17]. Soon after the discovery of new host of CBCVd, quite intriguing issues that remain to be clarified about the transcriptional response of hop during the symptom development process. In this context, a previous study was mainly focused on the precise measurement of target transcript of CBCVd-derived vd-sRNAs by small/micro RNA target prediction computational tools [6]. To disentangle the CBCVd-hop interaction, in the current study, we investigated the systematic and comprehensive changes in the transcriptome of hop leaves following infection by the biolistically inoculated virulent variant of CBCVd.
Plants have evolved series of defense mechanism via signal transduction pathway to defend against pathogens. The plant innate immune system such as pattern-triggered immunity (PTI) and effector-triggered immunity (ETI) is conserved and constitute the first line of defense [60]. The PTI or ETI is triggered by the activation of various membrane-associated receptor-like kinases (immune receptor) upon the perception of non-self components of pathogenic origins via MAPK (mitogen-activated protein kinase) and CDPKs (calcium-dependent protein kinases), signalling cascades that act downstream of receptor complexes. The transducing signalling events eventually lead to transcriptional reprogramming, induction of the PR proteins, and hypersensitive response by the generation of reactive oxygen species (ROS) [61]. Our results revealed the extensive activation of signaling pathways genes (e.g., MAPK3 and CDPKs) as well as prominent marker genes that are associated with innate immunity (e.g., PR proteins, 1,3-beta-glucanase, and ROS biogenesis genes), which is in agreement with the results of a previous study [27]. Since viroids (CBCVd) lack the coding capacity, thus we cannot assume the effectors-based triggering of immune response in the host plant. Nevertheless, our DEGs data illustrated that the hop possesses the effective perception system, which can easily discriminate and recognise CBCVd infectious RNA. In the animal system, the double-stranded RNA (dsRNA)-activated protein kinase (PKR) binds to viral dsRNA and invoke innate antiviral responses [62]. The previous study focusing on the RNA-binding activities of viral and host proteins in plants suggested that the PKR homolog in plant P58 IPK plays an important role in viral pathogenesis. The silencing of P58 IPK caused massive cell death in Nicotiana benthamiana and A. thaliana plants upon infection with the tobacco mosaic virus (TMV) [63]. Intriguingly, the P58 IPK gene was not found to be differentially regulated in CBCVd-infected samples. Therefore, the role of P58 IPK in triggering innate immunity against CBCVd-infection in hop is uncertain. Nevertheless, among several elevated protein kinases in our data set (Table S3) yet-to-be-discovered candidate player might play the dominant role in binding and perception of CBCVd and the further activation of defense response.
In plant cell, chloroplast is a dynamic organelle, which is not only involved in photosynthesis, but also actively participates in the synthesis of numerous compounds, interorganelle signaling, and regulation of cell-death programs as a part of plant immune system [64]. Several chloroplasts associated proteins are transcribed by nuclear genes and viroid replication in the nucleus cause triggering of chloroplast-based signaling cascades resulted in the impairment of the structural (thylakoid membrane abnormalities and paucity of grana) and functional changes in the chloroplast that is associated with the process of photosynthesis [64,65]. Several studies using microarray and RNA-seq approaches revealed that viroid infection could modulate the efficiency of photosynthetic rate by down-regulating the chloroplast-and photosynthesis-related genes [20,23,27]. In this study, we observed that many genes that are associated with chlorophyll metabolism and photosynthesis pathway such as photosystem I subunit l, ferredoxin-oxidoreductase, 2-oxoglutarate, and Fe(II)-dependent oxygenase superfamily, glyceraldehyde-3-phosphate dehydrogenase, ferredoxin, etc. were down-regulated, which might be correlated to symptom development (chlorosis, stunting or mosaic) in CBCVd-infected hop plants.
Numerous elegant studies have demonstrated that plant hormones play prominent roles in the plant life cycle and act as an important regulator of metabolic pathways that are related to plant growth, development, and abiotic/biotic stress responses [66,67]. Viroids or viruses infections are frequently associated with phytohormone alterations, which causes the disruption of host cellular physiology resulting in the developmental disorders [23,68]. In our study, the expression level of several genes that are associated with hormone signal transduction pathways showed the altered expression. In general genes involved in indole-3-acetic acid (IAA), ethylene (ET), gibberellins (GA), jasmonate (JA), brassisteroid (BA), and abscisic acid (ABA) were up-regulated. Similarly, genes involved in IAA, ET, and BA-biosynthesis and responses were induced in PSTVd-infected tomato plants [69]. Salicylic acid (SA) has recently received attention in rendering basal resistance of tomato plants against citrus exocortis viroid (CEVd) [70]. However, changes of expression of genes that are involved in SA metabolism or response were not observed. A similar result was observed in PSTVd-infected tomato plants, which showed the altered transcript levels for several genes involved in GA biosynthesis and BR signaling, but not other genes that are associated with SA or JA dependent pathways [23].
Translational reprogramming is a critical component of the cellular response, which is required for rapid cellular adaptation under different types of stress conditions, including pathogen infection [27]. The pathogen utilizes host translation apparatus to their own replication and intercellular movement, and often equated with pathogen strategy to circumvent the mRNA competitiveness [27]. Intriguingly, the marked up-regulation of genes encoding ribosome biogenesis proteins were observed, suggesting that CBCVd enhanced the translation process of hop for their own repertoire. The mounting evidence indicates that viroids and viruses exploit or interfere the plant ubiquitin/proteasome system and heat shock proteins (HSPs) for biotrophic interactions [6,71,72]. In this context, previous in-depth studies illustrated that viruses can actively manipulate the ubiquitination machinery and/or the proteasome to promote their replication and movement [73]. In addition, viruses can direct the ubiquitin proteasome system to the modification of new targets, such as argonaute to enhance or suppress gene silencing in order to suppress innate immunity [72]. In our study, we found that several genes (28 genes) associated to protein degradation via ubiquitin pathway were upregulated as compared to the down-regulated genes (13 genes) in CBCVd-infected hop. Similarly, HSP (HSP70 and HSP90) family homologs were shown to be induced several folds in many plant interaction studies with RNA viruses and viroids [74,75]. It has been suggested that HSP chaperones are targeted by RNA viruses for replication, cell-to-cell movement, and regulation of host defense directly or indirectly through interactions with DnaJ [76]. It is interesting to note that several classes of HSP (HSP33, HSP70, HSP90) transcripts are induced to higher level in CBCVd-infected hop, reinforcing the role of HSPs during the CBCVd-infection cycle. Several lines of evidence indicate that transcriptional factors (TFs) play a pivotal role in the signal transduction process in plants against abiotic and biotic response [77]. Intriguingly, the gene expression level of major families of TFs, such as bZIP, WRKY, MYBs, MYC, bHLHs, AP2/EREBP, and ERFs were significantly altered in CBCVd-infected plants (Table S3), which is consistent with the results of previous study of PSTVd-infection in potato [24].
In our study, majority of genes involved in primary and secondary metabolism pathways were fluctuated in response to the CBCVd-infection. At the link between primary and secondary metabolism, phenylalanine ammonia-lyase (PAL) serves as first enzyme that is involved in biosynthesis of secondary metabolites, such as flavonoids, anthocyanins, prenylflavonoids, and bitter acids in hop [78]. Flavonoids constitute a large group of phenolic secondary metabolites and have been shown to have diverse biological functions, such as symbiotic nitrogen fixation, pollination and floral pigmentation, defense against different biotic and abiotic stresses, etc. [79]. Nonetheless, the bitter acids (humulone or α-acid and lupulone or β-acid) and prenylflavonoids (xanthohumol and desmethylxanthohumol) secondary metabolite constituents of hop is widely used in brewing industry as a source of flavour-active secondary metabolites [78]. Notably, several genes related to general flavonoids and prenylflavonoids were influenced by CBCVd-infection in hop. The genes involved in flavonoids biosynthesis pathway were up-regulated, which indicated the stimulation of defensive substances in hop upon CBCVd-infection. Similarly, the regulatory genes encode for TFs (HlWRKY1 and HlWDR5) involved in prenylflavonoids biosynthesis pathway were up-regulated, which was in agreement with a previous study [30]. In summary, the analysis of differential expression of CBCVd-infected hop will facilitate the investigations of the detailed mechanisms of plant responses to viroid infection and contribute to a better understanding of development of strategies to combat viroid diseases in hop production.  Figure S3: Histogram presentation of clusters of orthologous groups (COGs) classification of unigenes (A); and differentially expressed genes in CBCVd-infected compared with mock inoculated hop (B). Figure S4: The analysis of protein-protein interaction network of CBCVd-responsive differentially expressed genes (DEGs) in hop by STRING (edge value cutoff of 0.70). The lines indicate network-predicted interactions. The regulatory network shows the DEGs of hop homologous to Arabidopsis thaliana. Table S1: Primers used for RT-PCR detection and qRT-PCR analyses. Table S2: Correlation analysis of sequencing libraries.