Profiling Plant circRNAs Provides Insights into the Expression of Plant Genes Involved in Viral Infection

1. Introduction

Circular RNAs (circRNAs) have recently emerged as important regulators of gene expression in various organisms, including plants. Unlike linear RNAs, circRNAs possess a covalently closed loop structure, rendering them resistant to exonucleases and exhibiting remarkable stability. While circRNAs have been implicated in diverse biological processes, including development and stress responses, their role in plant-virus interactions remains poorly understood. Viral infections pose significant threats to plant health and agricultural productivity worldwide, making it imperative to elucidate the molecular mechanisms underlying these interactions.

In this study, we investigate the dynamics of circRNA expression in the plant model system Arabidopsis thaliana (A. thaliana) during viral infections, focusing on different distinct viruses: Turnip rosette virus (TRoV) as a host virus, the small circular satellite RNA of the Lucerne Transient Streak Virus (scLTSV) and Rice yellow mottle virus (RYMV) as a non-host virus. By characterizing changes in circRNA expression profiles and analyzing their potential regulatory functions, we aim to shed light on the intricate interplay between circRNAs and viral pathogens in plants.

We hypothesize that viral RNA influences plant circRNA profiles by potentially modulating regulatory processes, including transcription. We evaluate circRNAs and their predicted impact on the dysregulation of plant cellular genes in response to viral infectious processes. We profile plant circRNAs induced by different plant viruses and circular virusoid RNA, determine which circRNAs react in response to virus infection, viral components, and circular satellite RNAs. We propose that the induced circRNAs could be potentially involved in viral replication, plant defense mechanisms against viruses, and viral disease symptoms. We catalog other virus/viral component-induced circRNAs involved in regulatory functions such as gene transcription, RNA splicing, etc. We examine the impact of the circular satellite RNA of the scLTSV virusoid on the plant endogenous circRNA profile and regulatory processes. We investigate the impact of a non-host virus (RYMV) on the profile and nature of induced plant circRNAs, and finally we determine whether the expression of an exogenous circRNAs, such as a virusoid, influences the expression of endogenous circRNAs in the host plant.

Our study offers a unique contribution to our understanding of plant circRNAs, particularly regarding their interactions in the context of viral infections.

2. Materials and Methods

Host Plant and Virus Model Systems Used in the Current Study

A. thaliana is an ideal host model system for studying circRNAs due to its well-characterized genome, short life cycle, and amenability to genetic manipulation. Furthermore, our understanding of circRNAs in plants, although limited, primarily stems from studies conducted in A. thaliana. TRoV and the scLTSV were chosen as model systems in this study due to their compatibility with A. thaliana. TRoV infects and replicates in A. thaliana, while scLTSV acts as a circular satellite RNA. Additionally, RYMV was selected as a non-host virus to elucidate host-specific responses to viral infections.

Plant Propagation, Inoculation with Virus, RNA Extraction and RT-PCR

A. thaliana seeds were pre-germinated on minimal agar plates (0.8%) for 2 days at 4°C in the dark, followed by growth in a greenhouse under a 16-hour light/8-hour dark cycle at 23°C. Seedlings at the 2-leaf stage with root length approximately 1 mm were transferred to soil. Four-week-old A. thaliana plants were inoculated with TRoV inoculum, and leaves were harvested at 5-, 10-, and 15-18-days post-inoculation to track viral systemic movement.

Total RNA was extracted from healthy and infected A. thaliana plants under RNase-free reaction conditions as specified in (Sambrook & William Russell, 2001). 30 mg of Arabidopsis leaves was used per RNA sample and each sample was subjected to phenol / chloroform extraction using TES buffer (200 mM Tris pH 7.4 and 25 mM EDTA) followed by ethanol precipitation in the presence of 1/10th volume 3M sodium acetate (pH 5.2) and 2.5 volumes of 95% ethanol. The extracted RNA samples were dissolved in nuclease-free water followed by treatment with DNase for the removal of contaminating genomic DNA as per the following reaction conditions: 15 μl of total RNA, 2.5 μl of 10X reaction buffer, 2 μl of RNase inhibitor (40 units), 1 μl of DNase (2000 units/mL), and nuclease-free water to a final volume of 25 μl. Samples were incubated at 37°C for 30 min and subsequently, the DNase was removed by phenol/chloroform extraction followed by ethanol precipitation of the remaining total RNA. The total RNA was dissolved in nuclease-free water, and the purity of the RNA extracts was estimated by measuring the ratio of A260/280 absorption with a Nanodrop-1000 spectrophotometer (NanoDrop Technologies Inc., USA). Quality and integrity of the RNA were validated by 1% (w/v) agarose gel electrophoresis using Tris/borate/EDTA (TBE) buffer. Subsequently, RNA samples were subjected to RT-PCR.

Total RNA samples purified as described above (200 ng) were used as templates for cDNA synthesis. A 10 μM reverse primer specific to the 3' end of the coat protein ORF of TRoV (Table 1) was used for the reverse transcriptase reaction as follows: samples containing the template and primer mixtures were denatured at 75°C for 2 minutes and then quickly chilled on ice for 1-2 minutes. The following reagents were then added to each sample: 2.5 μl of 10x AMV Reverse Transcriptase (RT) buffer (NEB), 1 μl of (10 mM) dNTPs, 1 μl of ribonuclease inhibitor (40 U/μl), and 2.5 μl of AMV RT (200 units) in a 20 μl reaction. The mixture was then incubated for 1 hour at 42°C. Subsequently, 2 μl of the prepared cDNA generated as described above was used as template for the next PCR step wherein it was subjected to PCR amplification to detect the desired PCR product specific for the virus in question. The PCR reaction contained the following components: 12.5 μl of 2x Taq froggamix (PCR master mix), 2 picomoles each of reverse and forward primers, 2 μl of cDNA, and nuclease-free water to a final volume of 25 μl. The reaction mix was overlaid with mineral oil and subjected to PCR in a thermal cycler (Illumina Novaseq ^TM 6000) under the following conditions: initial denaturation at 94°C for 5 minutes (1 cycle), followed by 35 cycles of denaturation at 94°C for 45 seconds, annealing at the specific primer annealing temperature for 45 seconds, and extension at 72°C for 1 minute, with a final extension at 72°C for 10 minutes (1 cycle). After amplification, PCR products were analyzed by gel electrophoresis.

Transformation of Agrobacterium Tumefaciens and Generation of Transgenic A. thaliana Plants

A. tumefaciens strain GV3103 was made competent, followed by separate calcium chloride mediated transformations with the pBI121-CP-TRoV construct, pCambia 1300-RYMV construct, pCambia 1300-scLTSV construct, and empty pBI121 vector positive control, using kanamycin for selection of transformants. Specifically, 100 ul of of the competent A. tumefaciens was incubated with 50 ng of plasmid DNA and the mixture was frozen in liquid nitrogen. All subsequent steps for Agrotransformation were carried out as per (Sambrook & William Russell, 2001). Transformed Agrobacteria were resuspended in 50 mL solution (5% sucrose, 0.05% Silwet-77 solution). The culture was then sprayed on the young flowers of A. thaliana. After four weeks, seeds from these transformed plants were harvested and germinated on agar plates containing Murashige & Skoog (MS) media with Kanamycin for the selection of transgenic plants. The selected transgenic plants were screened for the presence of the CP-TRoV or RYMV or scLTSV by first total RNA extraction followed by RT-PCR as described above using primers specific for the CP-TRoV, RYMV, and scLTSV, respectively (Table 1). Seeds from one line of transgenic A. thaliana were selected and screened by germination on 20 μg/ml Hygromycin in 1/2 x MS medium. Germinating plants were considered putative transgenics and were transferred to soil. Resistance to the antibiotic Hygromycin was used as a selection marker to distinguish between transgenic and non-transgenic plants.

Total Plant Protein Extraction, SDS-PAGE and Western Analysis

0.15g of Arabidopsis leaf (fresh) samples were added to autoclaved mortars with liquid nitrogen and ground to a fine powder. The powder was transferred to Eppendorf tubes containing 200 μl protein extraction buffer (50 mM Tris-HCl [pH 7.5], 15 mM MgCl2, 18% glycerol, 0.1 mM PMSF, 7M urea, 2% Triton X-100, and 0.5% NP40). Samples were thoroughly mixed by vortexing and then heated with 5x SDS-PAGE loading dye to 90-100°C for 7 minutes. The samples were permitted to equilibrate to room temperature for a few minutes before loading. Estimation of total protein concentrations was performed using Bradford assay.

SDS-Polyacrylamide Gel Electrophoresis (PAGE) was conducted using a 4-20% glycine precast gel (BioRad), as specified in Sambrook et al., 2001. A volume of 35 μl of each prepared sample, containing between 14.35 µg and 17.15 µg of protein based on the Bradford assay, was loaded into each well. The gel was subsequently subjected to electrophoresis under a constant voltage of 200V for a duration of 50 minutes. Subsequent to SDS-PAGE, the protein samples were transferred onto nitrocellulose membrane as per the procedure stated in Sambrook et al., 2001. 1:1000 dilution of the primary antibody sample raised in rabbits was used followed by treatment with 1:3000 dilution of the goat anti-rabbit secondary antibody conjugated to alkaline phosphatase. Detection was performed by incubation with alkaline phosphatase substrate solution in the dark for 3-5 minutes.

Methodology for Observation of Morphological Symptoms A. thaliana Plants in Response to Viral/Viral Component Infection

Analyses of A. thaliana symptoms were conducted to monitor plant growth and discern morphological distinctions among three different groups: healthy, transgenic, and infected plants. All plants were simultaneously grown under the same conditions. First, all seeds were grown on 1/2 x MS medium plates for germination and placed in the growth chamber with 16 hours of light and a 22°C growth chamber program. Once the plants developed roots, they were transferred to soil. Then, the plants were examined under the magnifying stereo microscope to capture images and identify any visible symptoms. Additionally, plant height was measured across six different plants for each sample.

3. Bioinformatic Analysis

Library Construction and Sequencing by LC Sequencing Company

Total RNA was extracted from Arabidopsis leaves, and its quantity and purity were assessed using the Bioanalyzer 2100 with an RNA 6000 Nano LabChip Kit (Agilent, CA, USA). Only samples with an RNA Integrity Number (RIN) greater than 7.0 were used for further processing. Approximately 5 μg of total RNA underwent ribosomal RNA depletion following the Ribo-Zero Plant Removal Kit (Illumina, CA, USA) protocol. The remaining RNA was treated with RNase R (Epicentre Inc, Madison, WI, USA) at 37°C for 30 minutes to remove linear RNAs, thus enriching for circular RNAs (circRNAs). The enriched circRNA fractions were then fragmented using the NEBNext® Magnesium RNA Fragmentation Module (NEB, USA). The fragmented RNA was reverse transcribed into cDNA using SuperScript™ II reverse transcriptase (Invitrogen, USA). The second strand was synthesized using E. coli DNA polymerase I (NEB, USA), RNase H (NEB, USA), and a dUTP solution (Thermo Fisher, USA). An A-base was added to the blunt ends to facilitate adapter ligation. Single- or dual-index adapters with T-base overhangs were ligated to the A-tailed DNA fragments, and size selection (300-600 bp) was performed using AMPureXP beads. The ligated products were PCR amplified under the following conditions: initial denaturation at 95°C for 3 minutes, followed by 8 cycles of 98°C for 15 seconds, 60°C for 15 seconds, and 72°C for 30 seconds, with a final extension at 72°C for 5 minutes. The average insert size was 300±50 bp. The cDNA library was then sequenced using the Illumina Novaseq™ 6000 platform, generating 2×150 bp paired-end reads.

Initial sequencing produced millions of paired-end reads, which were filtered using Cutadapt (version: cutadapt-1.9) to remove adapter sequences, polyA/polyG stretches, reads with more than 5% unknown nucleotides (N), and reads with more than 20% low-quality bases (Q-value ≤ 20). Sequence quality metrics, including Q20, Q30, and GC-content, were evaluated using FastQC (version 0.11.9). Filtered reads were mapped to the Arabidopsis genome using Bowtie2 and Tophat2 (version 2.0.4). Unmapped reads were further analyzed with Tophat-fusion (version 2.1.0) to identify back-splicing junctions indicative of circRNAs. De novo assembly of mapped reads into circRNAs was performed using CIRCexplore2 (version 2.2.6) and CIRI (version 2.0.2). Differential expression analysis was conducted using the R package edgeR. CircRNA expression levels were estimated using the number of back-spliced reads, normalized as spliced reads per billion mappings (SRPBM). Twelve circRNAs were selected for qPCR validation based on their relevance in plant development and virus infection. Each 10 μl qPCR reaction contained 2 μl cDNA (50 ng), 400 nM primer mix, 2 μl nuclease-free water, and 5 μl of 2× Luna Universal qPCR Master Mix (NEB, USA) SYBER green. Real-time RT-PCR was performed to compare expression levels across different plant lines and infection conditions.

Bioinformatics Software Used.

Analysis Item	Software	Version
Quality control	FastQC	0.10.1
Adapter removal	Cutadapt	1.10
Genome mapping	Tophat	2.0.4
Back-splicing junction reads filter	Tophat-fusion	2.1
circRNA identification	CIRCexplore, CIRI	2.2.6, 2.0.2
Differential expression analysis	edgeR	NA
Interaction with miRNA	Targetmimics	NA

This comprehensive approach ensured the accurate identification and quantification of circRNAs, facilitating subsequent functional analyses.

circRNA Classification Based on Genomic Origins

CIRCexplore and CIRI, renowned circRNA detection tools, were utilized in this study. CIRCexplore employs unmapped reads from Tophat processed through Tophat-fusion (version 2.1) to identify back-splice junction reads, while CIRI2 specializes in datasets containing circRNA species. Both tools were used to detect circRNAs, leveraging their respective strengths. Tophat version 2.0.4 aligned reads to the reference genome (www.plants.ensembl.org), followed by de novo assembly of circRNA reads using CIRCexplore version 2.2.6. Back-splicing junction reads were identified among unmapped reads using TopHat-Fusion and CIRCexplore. Each sample yielded unique circRNAs. Categorization included circRNAs from exonic regions designated as circRNAs, those from intronic regions labeled as ciRNAs, and those from intergenic regions designated as intergenic circRNAs.

4. Results Obtained

Symptom Development in A. thaliana Plants Due to Viral Infection

Our investigation began with an assessment of symptom development in A. thaliana plants infected with TRoV (host), plants transgenic for TRoV capsid protein gene (CP-TRoV), and plants transgenic for the RYMV genome (non-host). Observations of plant morphology revealed notable differences between healthy and infected/transgenic plants. CP-TRoV transgenic plants exhibited reduced height compared to healthy controls, with characteristic black spots on leaf edges. TRoV-infected plants displayed stunted growth and leaf discoloration, progressing from white spots at 10 days post-inoculation (dpi) to purple spots and dry leaves by 15 dpi. Transgenic circ-LTSV and RYMV-infected plants with TRoV showed no visible symptoms but exhibited variations in plant height compared to controls. Having observed the above symptoms, we proceeded to investigate the circular RNA profiles of these plants to identify their putative molecular interactions in the context of viral infection.

Constructs Prepared in the Current Study

Detailed procedures for construct preparation, including cloning of CP-TRoV, RYMV, and circ-LTSV genes into binary vectors, are provided in the Supplementary Section. Successful cloning for each of the above constructs was confirmed through rigorous validation steps documented in the supplementary material.

RT-PCR Results

Detection of CP-TRoV, RYMV, and circ-LTSV Genes in Transgenic A. thaliana Plants

Total RNA was extracted from transgenic A. thaliana plants expressing CP-TRoV, RYMV, or circ-LTSV genes. RNA integrity was verified, and time-course experiments revealed the early presence of TRoV in the TRoV-infected plants, with RT-PCR confirming viral gene expression.

Presence of CP-TRoV in Transgenic CP-TRoV A. thaliana

To validate the presence of the CP-TRoV transgene in transgenic A. thaliana plants, an RT-PCR experiment was conducted (Figure 1). Forward primer (CP-TRoV_F) and reverse primer (CP-TRoV_R) specific to the CP gene were utilized for PCR amplification (Table 1). An 800-base pair (bp) product was detected exclusively in the transgenic plants, indicating the successful integration of the CP-TRoV transgene into the A. thaliana genome. In contrast, healthy A. thaliana plants (negative control) did not exhibit the presence of the 800 bp band. As a positive control, TRoV-infected A. thaliana plants yielded an RT-PCR product at the expected size (800 bp). These results confirm the intact expression of CP mRNA under the control of the 35S CaMV promoter in the pBI121 vector, in which it was cloned. Elongation factor was employed as a loading control to ensure the accuracy of the RT-PCR results.

Presence of Genomic RYMV Transgene in RYMV Transgenic A. thaliana

An RT-PCR experiment was conducted to verify the presence of the genomic RYMV transgene in transgenic A. thaliana plants (Figure 2). Forward primer (CP-RYMV_F) and reverse primer (CP-RYMV_R) specific to the CP-RYMV gene were employed for PCR amplification (Table 1). A 720-base pair (bp) product was detected exclusively in the transgenic plants, indicating the successful integration of the genomic RYMV transgene into the A. thaliana genome. In contrast, healthy A. thaliana plants (negative control) did not exhibit the presence of the 720 bp band. Similarly, healthy rice plants (negative control) did not demonstrate the presence of the 720 bp band. As a positive control, RYMV-infected rice plants yielded an RT-PCR product at the expected size (720 bp).

Confirmation of Circular scLTSV (circ-LTSV) in Transgenic A. thaliana Infected with TRoV

An RT-PCR experiment was conducted to verify the presence of head-to-tail circ-LTSV in transgenic A. thaliana plants infected with TRoV (Figure 3). Forward primer (circ-LTSV-F) and reverse primer (circ-LTSV-R) specific to circ-LTSV were utilized for PCR amplification (Table 1). A 322-base pair (bp) product corresponding to circ-LTSV was detected in transgenic A. thaliana plants infected with TRoV as well as in positive control samples. Additionally, bands at 522 bp and 844 bp were observed, representing concatamers resulting from the production of circ-LTSV dimer and trimer forms at 322 bp intervals due to the rolling circle replication of scLTSV. The presence of a head-to-tail dimer of sc-LTSV is crucial and a minimum requirement for expressing the circular scLTSV genomic RNA in transgenic plants, facilitated by its self-cleaving ribozyme action through rolling circle replication (Abouhaidar & Paliwal, 1988). Hence, cloning of dimer sc-LTSV was preferred over the monomer form. Elongation factor was used as a loading control to ensure the accuracy of the RT-PCR results.

Confirmation of scLTSV (circ-LTSV) in Transgenic A. thaliana

An RT-PCR experiment was conducted to verify the presence scLTSV in transgenic A. thaliana plants (Figure 4). Forward primer (LTSV_sat_F) and reverse primer (LTSV_sat_R) specific to sc-LTSV were utilized for PCR amplification (Table 1). A 322-base pair (bp) product corresponding to scLTSV was detected in transgenic A. thaliana plants as well as in positive control samples.

Confirmation of the Absence of CP-TRoV in scLTSV Transgenic A. thaliana

Since the scLTSV transgenic were generated using a construct harboring both CP-TRoV and scLTSV, an RT-PCR experiment was conducted to verify the absence of CP-TRoV in sc-LTSV transgenic A. thaliana plants (Figure 5). Forward primer (CP-TRoV_F) and reverse primer (CP-TRoV_R) specific to CP-TRoV were utilized for PCR amplification (Table 1). No 800-base pair (bp) product corresponding to CP-TRoV was detected from the transgenic plants, indicating the absence of CP-TRoV in the scLTSV transgenic A. thaliana. As expected, positive control samples (TRoV-infected Arabidopsis plants) showed a band at the expected size (800 bp). Additionally, scLTSV transgenic A. thaliana plants and healthy A. thaliana plants (used as negative controls) did not exhibit the presence of the 800 bp CP-TRoV PCR product. Elongation factor was used as a loading control to ensure the accuracy of the RT-PCR results.

Western Analysis for the Detection of CP Expression in Plants Transgenic for CP-TRoV

A Western blot analysis was conducted to evaluate the expression of the capsid protein (CP) of TRoV in transgenic A. thaliana plants (Figure 6). The results revealed the presence of a distinct 30 kDa CP-TRoV band in the virus-infected positive control sample, demonstrating the specificity and efficacy of the TRoV polyclonal antibody in detecting CP-TRoV expression. The Western blot results indicated that the protein extraction method and Western blotting conditions were optimized for accurate detection. Notably, CP-TRoV transgenic A. thaliana plants also exhibited the CP-specific 30 kDa band, confirming the expression of TRoV capsid protein under the control of the CaMV 35S promoter. In contrast, negative control samples (healthy plants) displayed no equivalent band, further validating the specificity of the CP-TRoV detection.

Absence of RYMV Capsid Protein Expression in Transgenic A. thaliana as Demonstrated by Western Analysis

A Western blot analysis was conducted to assess the expression of the Rice yellow mottle virus (RYMV) capsid protein in transgenic RYMV A. thaliana plants, wherein A. thaliana servs which serve as a non-host for RYMV (Figure 7). Protein samples extracted from RYMV transgenic Arabidopsis plants and RYMV-infected rice plants (used as a positive control) were subjected to Western blotting. The results revealed a specific CP band at 30 kDa in the protein extract from RYMV-infected rice plants, indicative of CP expression and RYMV genome replication in its natural host. In contrast, the protein extract from RYMV transgenic Arabidopsis plants did not exhibit this CP-specific band, suggesting the absence of RYMV CP-protein expression. This observation implies the lack of RYMV replication in transgenic Arabidopsis plants, consistent with their status as non-hosts for RYMV. Negative control samples, including protein extracts from wild-type Arabidopsis thaliana and healthy rice plants, displayed no detectable CP band, further validating the specificity of the Western blot analysis.

CircRNA Sequencing

The advent of high-throughput sequencing (HTS) technology has revolutionized transcriptomic analysis by enabling the detection and characterization of various RNA species, including circular RNAs (circRNAs). CircRNAs, being non-polyadenylated RNA molecules, require specialized HTS methodologies for their accurate identification. Notably, the application of linear rRNA depletion strategies, such as RiboZero, has significantly facilitated circRNA sequencing (Salzman et al., 2012). In this study, we conducted comprehensive bioinformatic analyses to explore the potential functions of circRNAs in Arabidopsis plants under different conditions. Differential gene expression analysis, functional annotation, and pathway enrichment analysis were performed to elucidate the roles of circRNAs in response to various stimuli. Additionally, we present the results of circRNA sequencing and bioinformatic analyses, followed by detailed sequence mapping.

Results Obtained from RNA-SEQ Bioinformatic Analyses: Identification and Characterization of circRNAs: Analysis of Circular RNA Abundance Across Samples and GC Content

Analysis of circRNA sequencing data revealed the generation of more than 40 million valid reads for each sample/run of Arabidopsis thaliana leaves. The GC content was approximately 47%, and the Q30 (quality score) exceeded 92%. Using CIRI2, CIRCexplore and the R package edgeR, a total of 760 circRNAs were predicted from the sequencing data. We analyzed the count of circRNAs across all samples, including healthy Arabidopsis thaliana, TRoV-infected plants, and transgenic plants expressing CP-TRoV, scLTSV, and RYMV wherein varying counts were observed.

Genomic Distribution of Identified circRNAs

circRNAs were categorized based on their start and end positions in the genome, distinguishing between exonic, intronic, and intergenic origins. Predominantly, circRNAs originated from exons, with fewer arising from introns, suggesting that mRNA-producing regions are primary sources of circRNA generation. Among the samples, empty vector samples exhibited the highest count of genes generating intron circRNAs. In contrast, TRoV-infected, circ-LTSV, CP-TRoV, and control samples predominantly produced circRNAs exclusively from exonic or intronic regions. CP-TRoV uniquely produced intergenic circRNAs among the samples.

Figure 8 shows the approximate number of circRNAs originating from the exonic, intronic and intergenic regions of the plant genome. For instance, RYMV exhibits approximately 16 SRPBM circRNAs from exonic regions and around 3 circRNAs from intronic regions.

Principal Components Analysis

All samples underwent circRNA sequencing with two replicates, and comparative analysis revealed no significant differences between the two sequencing runs for all nine samples. Principal Component Analysis (PCA) was conducted on the normalized RNA-seq data encompassing 760 differentially expressed circRNAs across all nine samples. The resulting PCA plot provides a comprehensive representation of the conditions for all samples, demonstrating consistent representation across the two sequencing runs (Figure 9). To explore the overall expression patterns and relationships within these samples, Principal Component Analysis (PCA) was conducted on the normalized RNA-seq data encompassing 760 differentially expressed circRNAs across all nine samples. The resulting PCA plot (Figure 9) is instrumental in distinguishing negative controls (healthy, EV, and mock-inoculated samples) from the remaining samples within both batches. Overall, the plot provides a comprehensive representation of the conditions for all samples, demonstrating their consistent representation across the two sequencing runs.

Chromosomal Distribution of circRNAs

Analysis revealed circRNA distribution originating from across chromosomes 1 to 5, with a limited presence in mitochondria (Figure 10). Notably, the chloroplast chromosome expressed the highest number of circRNAs.

Diversity in circRNA Splicing Signals

The splicing signals at circRNA splice sites in our study exhibited diversity. Unlike the typical GT-AG splicing signal observed in animals, our sequencing results in Arabidopsis presented a variety of splicing signals, such as TA/GT and TT/GT, expressed in terms of AGCT due to the DNA sequencing method used (Figure 11).

Patterns of circRNA counts

In this study, circRNA counts in each sample were analyzed using the R program version 3.6.3 package ggplot2. Our data revealed diverse circRNA counts across different samples. For instance, ciRNA9 exhibited 30 counts at 15 dpi, decreasing to 10 counts at 10 dpi, while the mock-inoculated sample showed 60 counts. Another circRNA, circRNA2, had 23 counts in the non-host sample (RYMV transgenic A. thaliana) and 100 counts in the host sample (CP-TRoV transgenic A. thaliana). Similarly, ciRNA11 showed varied counts across different samples, with the highest count in the CP-TRoV-scLTSV double transgenic A. thaliana. All samples exhibited zero counts of ciRNA17, while healthy controls showed a high count for this circRNA.

Gene Ontology (GO) Functional Analysis

GO enrichment analysis was conducted for host genes encoding dysregulated circRNAs. Most host genes were associated with plant development and protein binding activity, suggesting a regulatory role of circRNAs in host gene regulation.

Characteristics of Chloroplast Derived circRNAs

In our study, the chloroplast chromosome exhibited the highest number of identified circRNAs, with a total of 288 circRNAs identified. To explore the diversity of these circRNAs, a subset was selected for further analysis. Notably, six circRNAs from this subset shared the same end position but differed in their start positions on the chloroplast chromosome. For example, ciRNA101 showed 5 counts in the circ-LTSV transgenic sample, while all other samples, including the control, showed zero counts. Similarly, ciRNA103 had 3 circRNA counts compared to the control's 8 circRNA counts. Despite these variations, all circRNAs shared the same ending position and differed by a few nucleotides in the start position (Figure 12). Furthermore, selected circRNAs derived from the chloroplast were chosen for qPCR validation (Figure 13). All selected circRNAs from the chloroplast for qPCR validation exhibited distinct start and end positions without any overlap.

Internal ORF Identification in circRNAs

Variations in circRNA2 counts were observed among transgenic and virus-infected groups. This circRNA exhibited elevated counts in all transgenic samples and reduced counts in TRoV-infected samples compared to the control group. Further analysis revealed that circRNA2 contains an open reading frame (ORF) with start and stop codons, suggesting its potential translatability to a protein. Similarly, circRNA205 (ATCG00020) showed count differences between non-host (RYMV) and host-virus samples. Notably, it exhibited 20 counts in the RYMV sample, while zero counts were recorded in the control and CP-TRoV/ scLTSV transgenic samples. Furthermore, circRNA205 contained an ORF region in its sequence, indicating potential translatability to a protein with similarities to a wheat protein known for monooxygenase activity.

Functional Insights of Identified circRNAs

This study profiled circRNAs from Arabidopsis thaliana in the presence of both host and non-host viruses for the first time. Several key findings include:

AT3G60240 (eIF4G): This gene generated circRNA expressed exclusively in the RYMV transgenic sample. The eIF4G protein, essential for translation initiation, reportedly also influences TCV multiplication (Yoshii et al., 2004)

AT4G21790 (TOM1): This gene generated circRNA expressed in the RYMV transgenics. TOM1 codes for a host factor necessary for TMV virus replication (Yamanaka et al., 2000).

AT3G05500 (SRP3): This gene showed low-level expression in the CP-TRoV transgenic sample compared to the negative control. Reportedly, overexpression of SRP3 results in improved drought stress tolerance and increased reproductive growth (Kim et al., 2016).

AT5G37475: This gene was expressed at low levels in the RYMV transgenic sample. It is involved in the formation of the eIF4G complex, which is crucial for initiating a subset of mRNAs involved in cell proliferation (TAIR, 2022).

CircRNA-miRNA Interactions

In this study, we investigated whether identified circRNAs could act as miRNA sponges. Dysregulated circRNAs were predicted as miRNA target mimics using Targetfinder and Miranda. CircRNAs were found to have one or multiple decoy sites for miRNAs (Figure 14). For example, circRNA291 (AT1G01520) might act as a decoy for miRNA ath-miR396b-5p. Similarly, circRNA236 (AT5G10550) might act as a decoy for two different miRNAs (ath-miR167a-5p and ath-miR167b). Multiple circRNAs from the same parental gene (AT1G01520) might act as decoys for one miRNA, such as circRNA ATH_circ04122, ATH_circ04121, and ATH_04124, all of which had miRNA ath-miR396b-5p binding sites.

Further examination revealed that mRNAs in the regulatory network of circRNA-miRNA-mRNA were mainly involved in plant growth and development. This suggests that circRNAs in plants respond to various stress conditions, including biotic and abiotic stresses, affecting miRNA expression profiles. While the circRNA-miRNA-mRNA network has been identified in plants, biochemical analyses are needed for further elucidation. Only one study reported the role of circRNA in plants, where exonic circRNA in A. thaliana was shown to regulate the splicing of its parental mRNA (Conn et al., 2017).

Enrichment Analysis

Enrichment analysis for all dysregulated genes was conducted using ShinyGO to extract essential biological processes, molecular functions, and KEGG pathways. Gene ontology and domain enrichment analyses were performed to elucidate the main biological functions of the identified circRNAs in our dataset. The analysis revealed expected biological functions, including plant development, RNA processing, and protein binding, with mRNA and RNA binding being the most significantly enriched function (Figure 15 and Figure 16).

Validation of Identified circRNAs by RT-qPCR

To confirm the origin of selected circRNAs, RT-PCR followed by real-time quantitative qPCR was conducted. Primers were designed using Primer 3 software, ensuring the detection of authentic circRNAs with divergent primers and linear RNAs with convergent primers. RNase R treatment was applied to remove linear RNAs, and primer efficiency was validated before qPCR. Efficient primers (90-110% efficiency with R2 > 0.98) were selected for qPCR validation. Gene expression changes observed in qPCR aligned with RNA-seq analysis, confirming the dysregulation of selected circRNAs (Figure 17). Strong correlation between qPCR and RNA-seq data validated the consistency of gene regulation patterns identified across both techniques.

5. Discussion

CircRNAs, emerging as vital gene regulators, encompass exogenous and endogenous types, yet their functions remain incompletely understood. While some roles, like miRNA "sponging," are recognized, much remains unclear, making differential expression analysis valuable. Our study sheds light on circRNA expression during virus infection and plant defense mechanisms. The escalating agricultural impact of plant viruses underscores the significance of understanding circRNA involvement. We explored circRNA profiles in virus-infected, uninfected, and transgenic plants expressing viral genes, revealing differential expression patterns across groups.

To elucidate A. thaliana's response to specific viral genes and non-host viruses, we employed transgenic plants. This novel approach unveils endogenous circRNA landscape alterations, providing insights into plant-virus interactions previously unexplored. Predictive analysis identified abundant and low-abundance circRNAs across experimental groups, revealing diverse origins in the A. thaliana genome. Dysregulated circRNAs, mainly from coding exons, underscored their potential regulatory roles in plant cells. GO enrichment analysis highlighted circRNA-associated host genes' involvement in diverse cellular processes, including stress response and mRNA binding. These findings suggest circRNAs' regulatory roles in plant defense mechanisms. Distinct circRNA expression patterns in virus-infected and transgenic plants point to potential plant genes implicated in infection processes. Notably, certain genes showed differential expression only in non-host plant samples, suggesting unique plant-virus interactions. PCA revealed distinct circRNA expression patterns among experimental groups, validating data reliability and highlighting circRNA alterations in response to viral infections. Dysregulated circRNAs identified in transgenic plants suggest exogenous circRNA (scLTSV) effects on endogenous circRNA expression, warranting further investigation into underlying mechanisms. Our data showed TA/GT and TT/GT as the highest splicing signals for the circRNAs, while Zhang et al. (2020) reported diverse splicing signals, with CT/AC and GT/AT being prominent. Another study on plant circRNAs, focusing on tomato plants under heat stress (Yang et al., 2020), found GU/AG as the predominant splice signal. Identification of dysregulated circRNAs with miRNA binding sites suggests potential miRNA sponging activity, implicating circRNAs in post-transcriptional gene regulation.

In our study, we found that circRNA2, involved in oxidoreductase and cytochrome complex assembly, showed varied expression patterns across different experimental conditions. Interestingly, it was upregulated in samples from RYMV, circ-LTSV, and CP-TRoV transgenic plants but downregulated in TRoV-infected samples compared to controls. The induction of circRNA2 seems to be triggered by various stressors on the plant, including mechanical stress induced by water inoculation and expression of scLTSV through the 35S promoter. This induction mechanism may involve RNA-RNA or RNA-DNA interactions, as scLTSV lacks coding capacity but can upregulate circRNA2, warranting experimental validation. Our findings suggest that transgenic circ-LTSV A. thaliana may induce resistance to helper virus TRoV infection, aligning with previous studies on viral satellite RNAs eliciting resistance to viral infections. Additionally, several cellular circRNAs, such as ATCG00030 (ciRNA106) and ATCG00130 (ciRNA11), are induced solely by the non-coding exogenous circLTSV.

Enriched pathways associated with photosynthesis and chloroplast-related processes in dysregulated circRNAs indicate intricate virus-chloroplast interactions influencing plant-virus dynamics. The exploration of chloroplast-derived circRNAs in our study revealed intriguing characteristics, particularly regarding their differential regulation in circ-LTSV transgenic plants. This selective disruption suggests a potential regulatory role for circ-LTSV in chloroplast-derived circRNAs. Additionally, the targeted approach for qPCR validation of circRNAs provided valuable insights into their structural diversity and potential functional roles, indicating nuanced regulation possibly associated with chloroplast-related processes or responses. The dysregulation of circRNAs across various samples, including TRoV-infected plants, RYMV (non-host) plants, and plants expressing viral genes, suggests complex interactions between viral components and plant circRNA pathways. While the exact mechanisms remain to be elucidated, hypotheses such as RNA-RNA interactions, miRNA sponging activity, or sequence similarities with plant miRNA or circRNA may contribute to circRNA dysregulation.

Our study's findings on circRNA expression patterns, especially in the context of viral infections and transgenic influences, shed light on the intricate regulatory landscape of circRNAs in plants. The genomic distribution of circRNAs and their differential expression patterns across samples provide valuable insights into the potential roles of circRNAs in plant biology and response to viral infections. The observations on plant morphology further underscore the impact of different viral infections on Arabidopsis plants, with notable differences in growth and symptom progression among TRoV-infected, RYMV-transgenic, and CP-TRoV-infected plants. These findings contribute to our understanding of the interactions between viruses and their plant hosts, potentially informing strategies for managing viral diseases in crops.

Unlike other studies that analyze mRNA to detect changes in circRNAs, our research focused exclusively on circRNAs. This approach allowed us to explore specific circRNAs that lack complementary sequences to miRNA and mRNA, suggesting a novel regulatory mechanism potentially involving RNA-RNA or RNA-DNA interactions, especially in circRNAs dysregulated in transgenic plants expressing the non-coding virusoid, circ-LTSV. Remarkably, our findings reveal that the majority of these circRNAs are localized to the chloroplast chromosome and intricately linked to the photosynthesis pathway, as evidenced by KEGG analysis. This association suggests a significant role for circRNAs in regulating critical photosynthetic processes, including light capture and energy conversion. This research not only deepens our comprehension of the functional dynamics of circRNAs within plant virology but also lays the groundwork for novel agricultural biotechnological applications aimed at enhancing crop resilience to viral threats, thus supporting food security.

Overall, this study provides extensive insights into the role of circRNAs in plant-virus interactions, highlighting their potential as targets for further research in plant biology, virology, and agricultural biotechnology.

6. Conclusion

In conclusion, this study provides novel insights into the regulatory role of circRNAs in plant-virus interactions, particularly in A. thaliana. Our findings highlight the diverse functions of circRNAs in modulating host gene expression and cellular responses to viral infections. We propose that further investigation into circRNA-mediated mechanisms could lead to the development of innovative strategies for enhancing plant resistance to viral diseases, thereby contributing to sustainable agriculture and food security.