Small RNA Profiling of Aster Yellows Phytoplasma-Infected Catharanthus roseus Plants Showing Different Symptoms

Micropropagated Catharantus roseus plants infected with ‘Candidatus Phytoplasma asteris’ showed virescence symptoms, witches’ broom symptoms, or became asymptomatic after their planting in pots. Nine plants were grouped into three categories according to these symptoms, which were then employed for investigation. The phytoplasma concentration, as determined by qPCR, correlated well with the severity of symptoms. To reveal the changes in the small RNA profiles in these plants, small RNA high-throughput sequencing (HTS) was carried out. The bioinformatics comparison of the micro (mi) RNA and small interfering (si) RNA profiles of the symptomatic and asymptomatic plants showed changes, which could be correlated to some of the observed symptoms. These results complement previous studies on phytoplasmas and serve as a starting point for small RNA-omic studies in phytoplasma research.


Introduction
Phytoplasmas are wall-less insect-transmitted bacteria. They strongly rely on the host during their entire lifetime and are not maintained as pure cultures in vitro, being provisionally classified in the 'Candidatus Phytoplasma' species [1,2]. Furthermore, they inhabit plant phloem and insect hemolymph. They were detected in hundreds of agronomically important plant species worldwide, which were associated with alterations in plant morphology and physiology due to modulating transcript and protein profiles, as well as changing hormone balances [3]. Phytoplasmas appear to have limited coding capacity in their genomes [4]. Moreover, they also appear to interfere with plant hormone metabolism, inducing witches' broom, virescence, phyllody, irregular flowering, decreased fruit size, plant yellowing, decline, and stunting [2]. Developmental processes and hormone metabolism are regulated by transcriptional factors, which can be further regulated by endogenous micro (mi) RNAs, and which are associated with RNA silencing [5]. The latter is also involved in the defense mechanisms of plants, which are directed by small interfering RNAs (siRNAs). RNA silencing is a conserved mechanism of innate immunity against viruses, fungi, oomycetes, and bacteria [6]. This basal defense mechanism can control the pathogen concentration in a sequence-specific manner. Increasing the number of pathogens has been reported to encode factors that suppress RNA silencing-based host defense mechanisms [7]. Furthermore, a recent study investigated the role of a candidate The phytoplasma quantification was performed on three DNA extracts for each plant via qPCR assay. Graphs show the relative quantification of the DnaA (stronger color) and the DnaB (lighter color); these were normalized on both CrUBQ11 and CrEF1α via LightCycler ® 96 Software (Roche, Basel, Switzerland ). Data are presented as the mean ± SD (n = 3). Statistical analysis was performed, via one-way ANOVA, only between the plants showing an infection phenotype. Multiple comparison tests were conducted using Tukey's test and were marked using a compact letter display (CLD). The samples were labeled with the same letter and were not significantly different from each other.

Phytoplasma Detection and Quantification
In order to verify the presence and identity of phytoplasmas in the selected plants, total nucleic acid extraction was carried out following the method of White and Kaper [31]. Total nucleic acid extracts were used for the small RNA HTS described below. In addition, these were employed for phytoplasma quantitative PCR analyses for the purposes of detecting two genes of AY, the DnaA gene (chromosomal replication initiator protein DnaA-AYWB_RS00005) and the DnaB gene (replicative DNA helicase DnaB-AYWB_RS00035). The analysis was conducted in triplicate for both genes using the periwinkle elongation factor 1-α (CrEF-1-α) and ubiquitin (CrUBQ11) as the internal control (Table S6.). The quantitative PCR was performed in a LightCycler 96 instrument (Roche, Basel, Switzerland). The amplification reactions were carried out using a Luminaris Color HiGreen qPCR Master Mix (ThermoFisher, Waltham, MA, USA) and 50 ng DNA as the template. A 2-step amplification protocol, as suggested by the manufacturer, was used. The results were then analyzed via the 2 ∆∆Ct method (specifically using the embedded software of the LightCycler 96 instrument (Roche, Basel, Switzerland)), using both CrEF-1-α and CrUBQ11 as the reference genes.

Small RNA HTS
A small RNA fraction was purified from the extracted total nucleic acid and used for a small RNA sequencing library preparation. In total, nine libraries were prepared using the TruSeq Small RNA Library Preparation Kit (Illumina, San Diego, CA, USA, USA). A modified protocol was used [32], whereby the sequencing was conducted using a single index on a HiScan2000 via UD GenoMed (Debrecen, Hungary), which was 50 bp long and with a single end reading. The FASTQ files of the sequenced libraries were deposited into the NCBI GEO database and can be accessed through the series accession number GSE213754. The phytoplasma quantification was performed on three DNA extracts for each plant via qPCR assay. Graphs show the relative quantification of the DnaA (stronger color) and the DnaB (lighter color); these were normalized on both CrUBQ11 and CrEF1α via LightCycler ® 96 Software (Roche, Basel, Switzerland). Data are presented as the mean ± SD (n = 3). Statistical analysis was performed, via one-way ANOVA, only between the plants showing an infection phenotype. Multiple comparison tests were conducted using Tukey's test and were marked using a compact letter display (CLD). The samples were labeled with the same letter and were not significantly different from each other.

Small RNA HTS
A small RNA fraction was purified from the extracted total nucleic acid and used for a small RNA sequencing library preparation. In total, nine libraries were prepared using the TruSeq Small RNA Library Preparation Kit (Illumina, San Diego, CA, USA). A modified protocol was used [32], whereby the sequencing was conducted using a single index on a HiScan2000 via UD GenoMed (Debrecen, Hungary), which was 50 bp long and with a single end reading. The FASTQ files of the sequenced libraries were deposited into the NCBI GEO database and can be accessed through the series accession number GSE213754.

Bioinformatic Analysis of the Small RNA Sequences
For the bioinformatics analysis, a CLC Genomic Workbench (Qiagen, CLCbio, Aarhus, Denmark) was used. After quality control, the trimming and removal of the C. roseus non-coding RNA (rRNAs and tRNAs) (RNAcentral) was performed. In addition, the small RNA reads were mapped to the C. roseus genome (ASM94934v1) [33]. The reads that were not mapped on the C. roseus genome were subsequently mapped on the genome of one strain of the aster yellows phytoplasma, including its published four plasmids [34]. The percentage of the mapped reads and the coverage of the AY genome by them was calculated in each library. Then, the samples were averaged out according to the phenotype. The C. roseus-mapped small RNAs were further analyzed for the presence of published miRNAs. For this analysis, first, the C.-roseus-annotated miRNA database [35] was used. For the annotation of the miRNAs that were not annotated there, we used miRBase (release 22.1) data for the Arabidopsis, tomato, and soybean miRNAs present at miRBase.
The relative expression values of the miRNAs and siRNAs were generated using the R (4.1.0) package pheatmap (1.0.12) [36]. The heatmaps were prepared using z-scores, which were calculated by RPM (reads per million). The hierarchical clustering was performed for the rows, and the gene clusters were defined by cutting the dendrograms in the groups. The principal component analysis (PCAs) of the miRNAs and siRNAs were made using the web tool ClustVis [37]. Only the miRNAs with a ≥5 mean RPM and siRNAs with a ≥50 mean in at least one group were taken into consideration for the data in the heatmaps and PCAs.

Small RNA Northern Blot
For the small RNA Northern blot analyses, the total RNA samples (5 µg) were fractionated on denaturing 12% polyacrylamide gels, which contained 8 M urea, and were transferred to Nytran NX membrane (GE Healthcare, Chicago, IL, USA) via a capillary method using 20xSSC, and were fixed by ultraviolet cross-linking. Membranes were probed with radioactively 32P-labeled LNA oligonucleotide probes (Exiqon, Vedbaek, Hovedstaden, Denmark), which were complementary to the mature microRNAs (Table S6), using a published protocol [38]. Briefly, 5 pmol of each oligonucleotide probe was end-labeled with [γ32P] ATP by using T4 polynucleotide kinase. The prehybridization of the filters was carried out in 50% formamide, 0.5% SDS, 5xSSPE, 5x Denhardt's solution, and 20 µg/mL of sheared denatured salmon sperm DNA. The hybridization was carried out at 50 • C for 2 h in the same buffer. The washing of the membranes was performed for 10 min, two times, with a washing solution containing 0.1%SDS and 2xSSC at the hybridization temperature.

Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) Analysis of the Predicted Target Genes
RT-qPCR analysis was carried out, as described in [39], in order to verify the differential expression of the selected genes that are targeted by miRNAs with seriously altered expressions in the phytoplasma-infected plants. The M-MLV reverse transcriptase (Promega, Madison, Wisconsin, USA) was employed in order to synthesize the cDNAs by using a random hexamer primer (Thermo Fisher Scientific Baltics, Lithuania) by following the manufacturer's instructions. For all of the genes of interest,~1.5 ng of a cDNA template was used in qPCR with the expression normalized to the ubiquitin gene. The qPCR reactions were performed using the SYBR Green Master Mix (Applied Biosystems, Foster City, CA, USA). This was achieved by using gene-specific primers (Table S6) and a RealTime PCR ABIPRISM StepOne sequence detection system (Applied Biosystem, Foster City, CA, USA) in three technical replicates and three biological replicates. The relative quantifications in gene expression were determined using the 2 ∆∆Ct method, thereby obtaining the fold changes in the gene expression, which were normalized to CrUBQ11, a polyubiquitin gene. The mean data obtained were analyzed by ANOVA (P 0.05), followed by Student's t-test.

Phytoplasma Concentration in Periwinkle Plants
C. roseus were initially infected with 'Ca. P. asteris'. The cuttings from this plant were maintained in micropropagation for years, then potted. In the following year, these plants showed three phenotypes, including virescence (Vir), witches' broom (WB), and no symptoms (NO) (Figure 1). Three plants (Vir1-Vir3) showed virescence, an abnormal development of green pigmentation occurring in the flowers, which is why the originally white flowers became green. The structure of the flowers was differently maintained; the symptoms of the Vir1 plant were the most severe. The flowers of this plant not only showed virescence, like the flowers of Vir 2 and Vir3, but also showed slight phyllody, which is when the flower structure starts to become leafy. Three plants (WB1-WB3) showed the witches' broom phenotype. Besides extensive branching, their internodes shortened, leading to a slightly stunted phenotype of the branched shoots. In addition, they did not flower, or their flowers turned green. The witches' broom symptom was the most severe in the case of WB1 followed by WB3; although WB2 showed severe branching, its leaves stayed bigger and their stunting was negligible. Six plants of the originally infected periwinkle descendant were asymptomatic. They had no unusual branching or stunting, and their flowers stayed white, with the original flower structure maintained. All the symptomatic plants and three plants from those in the asymptomatic group were chosen for the analysis (Figures 1a and S1).
To better verify the phytoplasma infection stage, the bacteria titer was determined. Aster yellows phytoplasma (AY) infections were detected in the Vir and WB plants showing symptoms and was only detected in incredibly low concentrations in one of the recovered plants (NO1) determined via qPCR assays calculating the delta-delta Ct, using the C. roseus elongation factor 1 and ubiquitin gene as references (Figure 1b). The phytoplasma DNA concentrations in the WB plants was higher than the concentrations in the Vir plants; however, the concentration was about the same in the Vir1 and WB2 plants. The concentration of phytoplasma DNA was very low in the NO1 plant, and no phytoplasma was detected either in the NO2 or NO3 plants.  (Table S1). The size distribution of the sequenced small RNAs showed a typical pattern. The majority of both the all-genome-mapped and the C. roseus-genome-mapped reads were 21 and 24 nt long, respectively. This demonstrates that they are products of the plant DICER enzymes, and that they also show that the RNA used for the library preparation was not degraded ( Figure S2a,b).

Small RNA HTS Profiling
A total of 11-20 million reads could be mapped to the C. roseus genome (Table S2), and the reads were analyzed for the miRNA and siRNA expression, while the remaining reads were mapped to the AY phytoplasma genome [34].

Characteristics of Phytoplasma-Derived Small RNAs
The small RNA reads that were not mapped to the C. roseus host genome were mapped to the AY genome. The analysis showed that for the Vir and WB plants, 1.7-3.4% and 5.1-12% of the total small RNA reads could be mapped to the phytoplasma genome, while for the asymptomatic plants, this value was 0.03-0.07% (Figure 2a).  Most of the small RNAs could be mapped to both of the two 16S rRNA gene copies that are encoded in the AY genome, but they were also mapped to different parts of the phytoplasma genome. The coverage of the AY genome by the small RNAs was 4.1% and 7% in the Vir and WB plants, respectively (Figure 2b). There was only a trace amount of AY-derived small RNAs in the NO plants, which were mapped to the conservative part of the 16S rRNA encoding gene. The WB plants had both a higher number of AY-mapped reads and a bigger coverage of the AY genome when compared to the Vir plants. To investigate whether these AY-derived small RNAs could have a specific size property reflecting their origin, the size distribution of the AY-mapped reads was checked ( Figure  S3). The size distribution of the AY-derived reads showed a constantly decreasing length. Most of the small RNAs could be mapped to both of the two 16S rRNA gene copies that are encoded in the AY genome, but they were also mapped to different parts of the phytoplasma genome. The coverage of the AY genome by the small RNAs was 4.1% and 7% in the Vir and WB plants, respectively (Figure 2b). There was only a trace amount of AY-derived small RNAs in the NO plants, which were mapped to the conservative part of the 16S rRNA encoding gene. The WB plants had both a higher number of AY-mapped reads and a bigger coverage of the AY genome when compared to the Vir plants. To investigate whether these AY-derived small RNAs could have a specific size property reflecting their origin, the size distribution of the AY-mapped reads was checked ( Figure S3). The size distribution of the AY-derived reads showed a constantly decreasing length.

The miRNA Profile of the Plants, Showing Different Symptoms
The miRNA profiling of the sequenced reads was performed on the reads that could be mapped to the C. roseus genome. The reads per million (RPM) were calculated for all the individual small RNA sequences. In order to be able to focus on the changes that could play a role in the symptom development, those having a high expression (at least 50 counts per million reads) and a significant expression change (higher than two times) were selected ( Table S2) Most of the small RNAs could be mapped to both of the two 16S rRNA gene copies that are encoded in the AY genome, but they were also mapped to different parts of the phytoplasma genome. The coverage of the AY genome by the small RNAs was 4.1% and 7% in the Vir and WB plants, respectively ( Figure 2b). There was only a trace amount of AY-derived small RNAs in the NO plants, which were mapped to the conservative part of the 16S rRNA encoding gene. The WB plants had both a higher number of AY-mapped reads and a bigger coverage of the AY genome when compared to the Vir plants. To investigate whether these AY-derived small RNAs could have a specific size property reflecting their origin, the size distribution of the AY-mapped reads was checked ( Figure  S3). The size distribution of the AY-derived reads showed a constantly decreasing length.

The miRNA Profile of the Plants, Showing Different Symptoms
The miRNA profiling of the sequenced reads was performed on the reads that could be mapped to the C. roseus genome. The reads per million (RPM) were calculated for all the individual small RNA sequences. In order to be able to focus on the changes that could play a role in the symptom development, those having a high expression (at least 50 counts per million reads) and a significant expression change (higher than two times) were selected ( Table S2)   The Vir2 and WB2 plants had about the same percentage of the AY reads and were the closest in this analysis. Specific changes in the miRNA expression levels were visualized by preparing a heatmap of the z-scores ( Figure 4).
The z-score changes were markedly different between the three groups of plants and could be grouped into eight clusters. Some of the clusters showed similar trends (clusters 4 and 5, and 6 and 8). It was interesting to note that although there were clusters showing the same tendency in the Vir and WB groups when compared to the NO samples (clusters 2, 3, and 6), there were also clusters in which the miRNA changes were more different between the Vir and WB groups than between the Vir and NO groups, or between the WB and NO groups (clusters 1, 4, 5, and 7). Detailed data for the aforementioned are presented in Table S3.
The z-scores' heatmap shows a differential regulation, but this representation does not take into consideration the absolute expression level of the individual miRNAs. To highlight the changes in relation to the expression levels, a new list was prepared grouping the samples according to the miRNA families (Table S4). In several cases, the different forms of the same miRNA family clustered into a different group. The summarized data of the most important changes, which showed the expressions of the most abundant members of the miRNA family, are shown in Table 1. The Vir2 and WB2 plants had about the same percentage of the AY reads and were the closest in this analysis. Specific changes in the miRNA expression levels were visualized by preparing a heatmap of the z-scores ( Figure 4).  Table S3.
The z-score changes were markedly different between the three groups of plants and could be grouped into eight clusters. Some of the clusters showed similar trends (clusters 4 and 5, and 6 and 8). It was interesting to note that although there were clusters showing the same tendency in the Vir and WB groups when compared to the NO samples (clusters 2, 3, and 6), there were also clusters in which the miRNA changes were more different between the Vir and WB groups than between the Vir and NO groups, or between the WB and NO groups (clusters 1, 4, 5, and 7). Detailed data for the aforementioned are presented in Table S3.
The z-scores' heatmap shows a differential regulation, but this representation does not take into consideration the absolute expression level of the individual miRNAs. To highlight the changes in relation to the expression levels, a new list was prepared grouping the samples according to the miRNA families (Table S4). In several cases, the different forms of the same miRNA family clustered into a different group. The summarized data of the most important changes, which showed the expressions of the most abundant members of the miRNA family, are shown in Table 1.  Table S3.
In order to validate the small RNA HTS analysis, the expression of some of the conservative miRNAs was investigated by a Northern blot hybridization, using LNA probes. For this validation, the RNA samples from individual plants were used.
The expression of miR156 showed a slight increase in the Vir and a slight decrease in the WB plants, but different miR156 species exhibited a slightly different pattern. The validation of its expression in the leaves shows a slight decrease in the Vir2 plants, but an increase in the WB2 plants, when compared to the NO2 plants. The Vir2 and WB2 groups were the less typical plants in their group, and this could explain the difference between the summarized small RNA HTS and the hybridization results ( Figure S4a,b). Two types of mir157 were identified, and only the less abundant one had an altered expression, showing a down-regulation in the Vir plants and an up-regulation in the WB plants. The Northern blot validation, while using RNA that was extracted from the stems, confirmed this altered expression ( Figure S4c). Several different types of miR159 were identified by the HTS with slightly altered expression patterns and rates, showing an overall decreased expression pattern in both the Vir and WB plants. The Northern blot analysis showed a slightly different pattern, and it also showed an increased expression in the phytoplasmainfected plants ( Figure S5). The mir162 expression was down-regulated in both the Vir and WB plants. The expression of the two types of miR165 and the several miR166 family members was altered. Although several miR166 members showed an increased expression rate, these changes were masked by the highly expressed miR166 species (showing reads that were higher than 100,000 RPM when compared to 5-50) (Tables 1 and S5). The Northern blot analysis was conducted for the miR166, and its decreased expression rate was validated ( Figure S6). The miR167 expression showed a significant down-regulation in the phytoplasma-infected plants, which was also confirmed by the Northern blot analysis ( Figure S5c). The expression of miR168 showed a decrease in the phytoplasma-infected plants. This decrease was more severe in the Vir plants and was validated by the Northern blot hybridization ( Figure S7). The miR168 was annotated incorrectly on the CroFGD website, whereby the reported mature form was actually the star strand. The mature form of the miR168 was found to be identical to that of Arabidopsis thaliana; however, when using this correct sequence of miR168, no target was found in the CroFDG's CDS collection. The miR168 regulates ARGONAUTE 1 (AGO1) mRNA, but the AGO1 sequence, which is annotated in the CroFGD (accession number CRO_T016560), was only partial. When using the more recent genome of C. roseus (ASM2450571v1), it was possible to identify the entire AGO1 gene sequence and to obtain more accurate predictions of the CDS and protein sequences. A new miRNA target prediction was performed using psRNATarget [40] with miR168 and CrAGO1, whereby the analyses confirmed the possibility of their interaction. The result indicated that the target region for the miR168 was identical to the one in A. thaliana, thus confirming the miRNA-target interaction in the CrAGO1 ( Figure S7c). The expression of the miR170 was decreased (Table S4), while the expression of the miR171 was increased in the phytoplasma-infected plants, which was validated by the Northern blot analysis (in this latter case of WB2 ( Figure S8)). The miR172 expression showed an increase in both the Vir and WB plants, but this change could not be fully validated either in the Vir or in the WB plants from the RNA that was purified from the stems ( Figure S9). The miR319 showed a slightly increased expression in the Vir plants, and a slight decrease in the WB plants, but only this later alteration could be validated ( Figure S5d). The miR390 and miR391 both showed a decreased expression in the phytoplasma-infected plants. Moreover, this was validated in the case of miR390 by the Northern blot analysis ( Figure S10). Although the signal was very faint, this result was in line with the low RPM; in addition, it was visible that, in the phytoplasma-infected plants, this signal was decreased.
Although the miR395 was identified with a low RPM, its expression changed very drastically (mostly in the Vir plants). The expression of the miR396 decreased due to phytoplasma infection, and this decrease was stronger in the WB plants (Table S4). Although a pronounced signal for this RNA was expected, only a faint signal during the Northern blot validation was obtained, due to not being able to unambiguously confirm its decreased level ( Figure S11). The miR398 showed a strongly increased expression in the phytoplasmainfected plants, which was validated by the Northern blot analysis ( Figure S12).
The expression changes in several currently described C. roseus-specific miRNAs [26] were found in the phytoplasma-infected plants. Among them, cro-novel-34 was highly expressed in the asymptomatic plants (about 1000 RPM) and was significantly induced (more than twice) in both the Vir and WB plants.

The siRNA Profiles of the Investigated Plants
Small RNAs, which were involved in the siRNA pathway, were also analyzed to identify whether phytoplasma infection could change their expression. The 24 nt long reads mapped to C. roseus were selected. To focus on the relevant changes, only the siRNAs presented with at least 50 RPM in one of the three groups of plants were analyzed. As an outcome, a total of 305 unique 24 nt long siRNAs were identified. Among these sequences, 22 belong to the ribosomal RNAs or transfer RNAs. The other siRNAs exhibited predominantly 5 adenine (about 66%), which is in accordance with the preferential recruitment of the AGO4 and AGO6 in Arabidopsis [41], i.e., the key proteins involved in the DNA methylation pathway. The PCA analysis showed that both the Vir and WB plants cluster together, where the WB plants were particularly closely related ( Figure 5). sequences, 22 belong to the ribosomal RNAs or transfer RNAs. The other siRNAs exhibited predominantly 5′adenine (about 66%), which is in accordance with the preferential recruitment of the AGO4 and AGO6 in Arabidopsis [41], i.e., the key proteins involved in the DNA methylation pathway. The PCA analysis showed that both the Vir and WB plants cluster together, where the WB plants were particularly closely related ( Figure 5). , the asymptomatic plant in which phytoplasma DNA was detected, clusters closer to the WB group than to the other two asymptomatic plants. The detailed analysis and heatmap obtained from the average RPM of the three plant groups were converted to z-scores. This showed that the siRNA clustered into 11 groups, with different expression patterns ( Figure 6). NO1, the asymptomatic plant in which phytoplasma DNA was detected, clusters closer to the WB group than to the other two asymptomatic plants. The detailed analysis and heatmap obtained from the average RPM of the three plant groups were converted to z-scores. This showed that the siRNA clustered into 11 groups, with different expression patterns ( Figure 6). circles highlights PCAs of the three groups of plants NO1, the asymptomatic plant in which phytoplasma DNA was detected, clusters closer to the WB group than to the other two asymptomatic plants. The detailed analysis and heatmap obtained from the average RPM of the three plant groups were converted to z-scores. This showed that the siRNA clustered into 11 groups, with different expression patterns ( Figure 6).   Table S5.
Among those phytoplasma-responsive siRNAs, 53 were significantly decreased (fold change ≥ 2), and only 1 was up-regulated in the Vir plants. In the WB plants, 16 siRNAs were significantly decreased, and 10 siRNAs were up-regulated (Table S5).
During the siRNA analysis, several detected siRNAs were mapped to the upstream region of the minovincinine 19-hydroxy-O-acetyltransferase (CrMAT-GenBank accession number AF253415) [24] and desacetoxyvindoline-4-hydroxylase-like (CrD4H-like-GenBank accession number GU363550) [42] genes ( Figure 7).  CrMAT and CrD4H-like are involved in the alternative pathway of vindoline synthesis, where their altered expression could affect the synthesis rate of vindoline, which is one of the main C. roseus products used in medicine. CrMAT has been reported to be active only in the root; thus, siRNAs targeting its promoter are possibly able to alter its expression and could explain its reduced expression in the leaves. In the CrMAT gene, siRNAs are mapped to the predicted promoter region that contains tandem repeats. The expression of the most abundant siRNAs decreased with 30% and 11% in the Vir and WB plants, respectively. The siRNAs mapped to a D4H-like gene also show similar behavior, although their combined expression is only slightly higher than 25 RPM in the control. The two main siRNAs are derived from the same sequence element, but with opposite orientation, showing their real siRNA origin (Figure 8).

Investigation of the TIA Pathway Key Regulator Expression during Phytoplasma Infection
Although, due to the lack of a proper annotation of the C. roseus genome, a prediction analysis of the possible miRNA target genes was not carried out, the expression pattern of certain key enzymes playing a role in the terpenoid biosynthetic pathway was investigated. ORCA2, ORCA3, and CR1 are key regulators of the TIA pathway. Meanwhile, the ORCA2 and ORCA3 promote the CR1 blocks, the transcription of terpenoid indole alkaloid (TIA) genes encoding molecules that are produced in plants under stress (Figure 7) [43]. They can be regulated by jasmonate-responsive AP2-domain transcription factors, which are regulated by the miR172. The mir172 expression was induced in the Vir and WB plants. Although there is no evidence of a direct correlation between the miR172 and the above three genes, we thought that it would be important to know how their expression was changed during phytoplasma infection. The expression of ORCA2 and ORCA3 were slightly induced in the Vir and WB plants, while the expression of CR1 was only induced in the WB plants (Figure 9).

Investigation of the TIA Pathway Key Regulator Expression during Phytoplasma Infection
Although, due to the lack of a proper annotation of the C. roseus genome, a prediction analysis of the possible miRNA target genes was not carried out, the expression pattern of certain key enzymes playing a role in the terpenoid biosynthetic pathway was investigated. ORCA2, ORCA3, and CR1 are key regulators of the TIA pathway. Meanwhile, the ORCA2 and ORCA3 promote the CR1 blocks, the transcription of terpenoid indole alkaloid (TIA) genes encoding molecules that are produced in plants under stress (Figure 7) [43]. They can be regulated by jasmonate-responsive AP2-domain transcription factors, which are regulated by the miR172. The mir172 expression was induced in the Vir and WB plants.
Although there is no evidence of a direct correlation between the miR172 and the above three genes, we thought that it would be important to know how their expression was changed during phytoplasma infection. The expression of ORCA2 and ORCA3 were slightly induced in the Vir and WB plants, while the expression of CR1 was only induced in the WB plants (Figure 9). [43]. They can be regulated by jasmonate-responsive AP2-domain transcription factors, which are regulated by the miR172. The mir172 expression was induced in the Vir and WB plants. Although there is no evidence of a direct correlation between the miR172 and the above three genes, we thought that it would be important to know how their expression was changed during phytoplasma infection. The expression of ORCA2 and ORCA3 were slightly induced in the Vir and WB plants, while the expression of CR1 was only induced in the WB plants ( Figure 9). MYB33, as a member of the MYB superfamily, plays a role in plant hormone responses during the defense against pathogens [44]. In addition, it exists under the miR159 regulation, and its expression was slightly induced in the WB plants ( Figure 10). The same MYB33, as a member of the MYB superfamily, plays a role in plant hormone responses during the defense against pathogens [44]. In addition, it exists under the miR159 regulation, and its expression was slightly induced in the WB plants ( Figure 10). The same trend of changes occurred, where the highest expression in the WB plants was observed for the superoxide dismutase, CSD2, which is the key enzyme in ROS detoxification and is regulated by the miR398 [45] (see Figure 10). trend of changes occurred, where the highest expression in the WB plants was observed for the superoxide dismutase, CSD2, which is the key enzyme in ROS detoxification and is regulated by the miR398 [45] (see Figure 10).

Discussion
In this work, the small RNA profile of nine originally phytoplasma-infected C. roseus plants was studied. While six plants showed phytoplasma symptoms of virescence and witches' broom, three plants were asymptomatic. The phytoplasma symptoms can show remission, especially after long-term micropropagation; this status is usually temporary, and symptoms usually re-emerge (A. Bertaccini, unpublished data). One of the asymptomatic plants was PCR positive but had significantly lower phytoplasma concentration when compared to the symptomatic plants. Quantitative PCR showed significantly high phytoplasma concentrations in the symptomatic plants, which was typically higher in the WB plants than in the Vir plants. The phytoplasma concentration in the Vir1 plant, which was the most symptomatic plant in the Vir group, and the WB2 plant, i.e., the least symptomatic plant in the WB group, were about the same, thus suggesting that the phytoplasma concentration in the leaves correlates with the symptom intensity, which is in agreement with previous reports [46]. For the study, instead of uninfected plants, the recovered plants were chosen as a control because they originated from the same plants, so their genetic expression was identical. The miRNA pattern of the plants can be altered by slight changes in the environment, and the aim was also to minimize the uncontrolled changes. The recovered plants that were the descendants of the same plants and underwent the same propagation, and were maintained at the same place, consequently differed from the symptomatic plants only in the presence of the symptoms, which turned out to be correlated very well with the AY concentration. In this case, the difference in the miRNA and siRNA profile compared to the symptomatic plants could be undoubtably

Discussion
In this work, the small RNA profile of nine originally phytoplasma-infected C. roseus plants was studied. While six plants showed phytoplasma symptoms of virescence and witches' broom, three plants were asymptomatic. The phytoplasma symptoms can show remission, especially after long-term micropropagation; this status is usually temporary, and symptoms usually re-emerge (A. Bertaccini, unpublished data). One of the asymptomatic plants was PCR positive but had significantly lower phytoplasma concentration when compared to the symptomatic plants. Quantitative PCR showed significantly high phytoplasma concentrations in the symptomatic plants, which was typically higher in the WB plants than in the Vir plants. The phytoplasma concentration in the Vir1 plant, which was the most symptomatic plant in the Vir group, and the WB2 plant, i.e., the least symptomatic plant in the WB group, were about the same, thus suggesting that the phytoplasma concentration in the leaves correlates with the symptom intensity, which is in agreement with previous reports [46]. For the study, instead of uninfected plants, the recovered plants were chosen as a control because they originated from the same plants, so their genetic expression was identical. The miRNA pattern of the plants can be altered by slight changes in the environment, and the aim was also to minimize the uncontrolled changes. The recovered plants that were the descendants of the same plants and underwent the same propagation, and were maintained at the same place, consequently differed from the symptomatic plants only in the presence of the symptoms, which turned out to be correlated very well with the AY concentration. In this case, the difference in the miRNA and siRNA profile compared to the symptomatic plants could be undoubtably coupled to the higher concentration of the phytoplasma. The phytoplasma-derived small RNAs were present in the phytoplasma-infected plants, raising the possibility of their use in pathogen diagnostics. However, as their size distribution did not show an abundance of 21-24 nt long small RNAs, they are very likely not products of specific DICER enzyme activity, but are more likely products of other RNase enzymes, which are active in the C. roseus host.
The PCA analysis shows that the miRNA pattern in the phytoplasma-infected symptomatic plants was different from those of the asymptomatic plants. Consequently, these different expression patterns can also play a role in the different observed symptoms.
The z-score display of the miRNA expressions showed eight distinguishable clusters. Their analysis suggests that there are differentially regulated miRNAs comparing not only the NO and phytoplasma-infected groups, but also comparing the WB and Vir groups, since their different expression is present in plants with different symptoms.
The changes in the miRNA expression were detected by the small RNA HTS, which were further investigated via the Northern blot analysis in some cases. While the small RNA HTS results were analyzed by the mean of the three different groups of plants showing different phytoplasma concentrations and symptoms, the Northern blot analysis was carried out with the use of the RNA extracts of individual plants via LNA probes that were used to detect the conserved members of the various miRNA families. The two methods showed the same trend in most but not in all cases. There are only a limited number of reports investigating the changes in the miRNA profile of the phytoplasmainfected plants. In these reports, a small RNA profiling of Mexican lime, which was infected with 'Ca. P. aurantifolia'; mulberry, which was infected by aster yellows phytoplasmas; paulownia, which was infected with paulownia witches' broom; and a witches' broomdiseased Chinese jujube was carried out and based on the data from only the one-one libraries, which were prepared from the infected and healthy plants [15][16][17][18][19]. Although in the infected libraries, the samples from the different individuals were pooled, these limited sources of data could explain the differences in their results. Additionally, small RNA profiling was conducted on the grapevine var. Chardonnay, which was infected by either aster yellows phytoplasma [20] or 'flavescence dorée' [21]. The different observations of these two works suggest that the effects of different phytoplasma strains in the same host species can be different. As the pathogen presence is limited to the phloem, its effect could be also different in the stems and in the distinct parts of the plant [47]. Keeping in mind the limits of these techniques, certain conclusions can be drawn from the results of the work presented here.
Several studies in annual and perennial plants have identified the conserved roles of miR156 and its targets, i.e., SBP/SPL genes, in guiding the switch of plant growth from juvenile to adult phases. The miR157 is thought, as is the case in the miR156 family, to target mRNA coding proteins containing the Squamosa promoter binding protein (SBP) box [48]. However, neither the miR157 nor the miR156 have targets that are predicted in C. roseus; instead, they are thought to target ten mRNAs in other plant species that are fine-tuning the plant developmental processes. Their altered expression in the phytoplasma-infected plants suggests that they could play a role in the specific symptom development, as it was found in most of the cases when miRNA profiling was investigated in phytoplasmainfected plants [15][16][17][18][19]21].
The mir159 and miR319 share a common target, a MYB transcription factor, which is annotated in C. roseus. Moreover, the miR319 targets are members of the TCP gene family, which participate extensively in plant development and in responses to environmental stresses. In sepals, petals, and anthers, it was found that the GAMYB and TCP proteins were expressed and that they directly interacted with regulating miR167, which creates a miR159-miR319-miR167 network [49]. The expression of the miRNAs in this network was found to change during phytoplasma infection, but its induction and repression were also both reported. In this network, the concentration of the miR167 was kept low, thereby resulting in a strong ARF6/8 expression, which, in turn, regulates many of the genes that are required for flower development, including auxin signaling regulation. In the experiments performed here, both the miR159 and miR319 were down-regulated in the phytoplasma-infected plants, thus leading to a down-regulation of the miR167. This results in a disturbance in the ARF-regulated flower development, which could result in the virescence phenotype. The expression of the MYB33 gene was shown to be increased in the WB plants, whereas the expression of the miR159 was down-regulated, thereby suggesting a direct miRNA-target correlation.
The mir165/166 family represents the most abundant miRNAs in plants, having nine members in A. thaliana and in their targeting of five members of the HD-ZIP III transcription factor family, as well as in controlling developmental processes, such as apical shoot development, through the expression of several auxin biosynthesis, transports, and response genes [50]. As they show differential patterns in the phytoplasma-infected plants, it is highly possible that their altered expression can shape the developmental changes that lead to specific symptom development. In the data presented here, in contrast to the other species [15,17,18], the expression of all the miR166 species, in total, was down-regulated. However, it is possible that this overall expression change was masked because of the high number of the reads that can be different in distinct tissues of the plant, thus resulting in an altered regulation mechanism.
In its canonical form, the miR162 targets DCL1; meanwhile, the miR168 targets AGO1, which are the main executor molecules of the miRNA pathway. For the miR162, no target could be predicted in the C. roseus, but its over-expression or blocking its effect on the target induced resistance in Oryza sativa. In addition, in its pathogen-susceptible variants, a decrease in the miR162 could be detected [51]. Its level was decreased in both the Vir and WB plants, which suggests an efficient pathogen counter-defense strategy, but this phenomenon should be further investigated for better clarity. The AGO1 level was finely tuned by a precise miR168 concentration, which also showed an altered expression during phytoplasma infection. In previous works, it was found that during virus infection, the increased miR168 concentration can translationally inhibit the AGO1 expression and can also decrease the activity of antiviral silencing [52]. During phytoplasma infection, the concentration of miR168 decreases, which can result in an increased AGO1 level and activity; however, this hypothesis should be further addressed.
The miR170 and mir171 target SCARECROW-like genes, which are a family of transcription factors that play a role in gibberellin-controlled radial patterning in roots. Their altered pattern in the phytoplasma-infected plants were detected and, in the case of the miR171, it was validated by the Northern blot analyses. The increased expression of the miR171 could be detected in the WB2 plants; meanwhile, in the Vir2, the less typical plant showing virescence, this increase could not be seen. The increased expression of the miR171 was reported previously from Ziziphus jujube and mulberry [16,19], but its down-regulation was found in paulownia [17].
The mir172 regulates the Apetala2 gene, and its overexpression could result in the early flowering, development of abnormal flowers, and altered leaf morphology of transgenic plants [53]. In phytoplasma-infected plants, its increased expression was observed but could not be validated with the Northern blot analysis. The reason for the failure could be that only a stem-specific RNA sample was probed, which could show a different expression pattern. It is interesting to note that in Z. jujube, a decrease was also noted in the miR171 level, which was found during phytoplasma infection [19].
The miR390 and miR393 regulates auxin-based lateral root development, which can be altered by different types of stresses [54]. It appears that during phytoplasma infection, their levels decrease, but the expression changes. This would be better investigated in the roots to find out the precise mechanism.
It is quite interesting that the miR395 level was highly altered in the phytoplasma-infected plants and that its increase was also reported in an aster yellows-infected grapevine [20]. Furthermore, it plays a role in sulphate assimilation and redox signaling [55]. As its level was markedly changed, its role could be important during the phytoplasma infection, as well as in the development of phytoplasma-specific symptoms.
The miR396 regulates targets that belong to the growth-regulating factor (GRF) family of transcription factors, thereby controlling the cell proliferation in Arabidopsis leaves [56]. Its decreased level could be found in the phytoplasma-infected leaves, suggesting its possible role in the development of the symptoms.
The miR398 targets Cu/Zn superoxide dismutase, which is responsible for scavenging reactive oxygen species and for being produced during various biotic stresses, including during bacterial, fungal, and viral infections. Interestingly, the level of the miR398 is usually up-regulated during viral infection; meanwhile, it could be either up-or down-regulated during bacterial infection [57]. Its increased expression was reported in Z. jujube [19], but its decrease was shown in the phytoplasma-infected mulberry [15]. Here, the up-regulation of the miR398 and in one of the SOD (CSD2) during phytoplasma infection were found, suggesting that this pathway can also play a role in the defense reaction of the plant.
The expression patterns of novel C. roseus miRNA were also changed. This change was remarkable in the case of the cro-novel-34. However, due to the lack of detailed annotations of the genome, it was only possible to speculate whether its induced expression could play a role in the symptom development during phytoplasma infection.
In the phytoplasma-infected plants, the induced level of ORCA2, ORCA3, and CR1 were found, which could result in an increased TIA level. ORCA2, ORCA3, and CR1 can be jasmonate-induced and, as they contain AP2-like motifs in their promoter, it is possible that they are further controlled by the miR172. In parallel, it was found that a decreased level of siRNAs-mapped to the MAT and D4H-like promoter, which could induce the expression of these genes-led to a similarly increased TIA level; however, this point should be further investigated.
Phytoplasma infection highly interferes with the hormone balance and is associated with developmental changes in plants. MiRNAs and siRNAs can play a role during this altered regulation. The pathogen itself is confined to the phloem in an uneven concentration. This spatially different expression of the pathogen can further change temporally during the infection cycle, thus making it particularly difficult to precisely monitor the molecular mechanisms lying behind the symptom expression. An miRNA-based hormone regulation can be one of the main components of this complex process. C. roseus is not only the model plant for phytoplasma research, but it is also widely used as a medicinal plant; this is why it is very interesting to know whether the production of its important alkaloids can be affected during phytoplasma infection. The alkaloid content of the phytoplasma-infected plants would be important to further investigate in the future not only to better understand this important process, but also to potentially increase the yield of the beneficial compounds.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/genes14051114/s1, Figure S1: Photos about the groups of periwinkle plants used for small RNA NGS, showing different symptoms; Figure S2: A column diagram of the size distribution of the sequenced small RNA reads in the nine sequenced small RNA libraries; Figure S3: A column diagram of the size distribution of the AY-phytoplasma-mapped reads in the sequenced small RNA libraries; Figure S4: Summarized data (a) and Northern blot-based validation of the miR156 (b) and miR157 (c) expressions; Figure S5: Summarized data (a) and Northern blot-based validation of the miR159 (b), miR167 (c), and miR319 (d) expressions; Figure S6 Table S1: Initial statistics of the sequenced small RNA libraries; Table S2: The RPM values of the miRNAs in the different libraries and their average RPM in the three different groups of plants, ordered according to the miRNA families; Table S3: The z-scores of the miRNAs in the three different groups of plants, showing the data for the detected miRNA and ordered according to the clusters shown in Figure 4; Table S4: Summary of the main changes in the miRNA expression of the phytoplasma-infected plants; Table S5: The z-scores of the siRNAs in the three different groups of plants, showing the data for the detected siRNA and ordered according to the clusters shown in Figure 6; Table S6: The sequences of the primers used for the RT-qPCR or Northern blot hybridization.