Probing in vivo RNA Structure With Optimized DMS-MaPseq in Rice

RNA transcripts form various secondary and tertiary structures that have a wide range of regulatory functions. Several methods have been developed to profile in vivo RNA secondary structure in prokaryotes and eukaryotes. These methods, such as dimethyl sulfate (DMS) mutational profiling with high-throughput sequencing (DMS-MaPseq), couple small chemical-mediated RNA modifications with next-generation sequencing. DMS-MaPseq, a powerful method for genome-wide and target-specific RNA secondary structure profiling, has been applied in yeast, mammals, Drosophila, and Arabidopsis thaliana, but not in crops. Here, we used DMS-MaPseq to conduct a target-specific and genome-wide profile of in vivo RNA secondary structure in rice (Oryza sativa). The DMS treatment conditions were optimized for rice leaf and root tissues. To increase the sequencing depth and coverage of low-abundance transcripts in genome-wide DMS-MaPseq, we used streptavidin-biotin depletion to reduce the abundance of highly expressed chloroplast transcripts during library construction. The resulting target-specific and genome-wide rice DMS-MaPseq data were of high quality and reproducibility. Furthermore, we used DMS-MaPseq to profile the in vivo RNA secondary structure of an OsmiR399 target region located at 5′UTR of OsPHO2, which participates in rice phosphate homeostasis. An unfolded RNA structure downstream of miRNA target site was observed in predicted in vivo RNA secondary structure, reminiscence of the TAM (Target Adjacent nucleotide Motif) involved in mRNA structure-mediated regulation in miRNA cleavage. Our study optimized DMS-MaPseq for probing in vivo RNA secondary structure in rice, facilitating the study of RNA structure-mediated regulations in crops.

In so-called RT stop methods, 3′ ends of cDNA mapped to the transcriptome indicate the single-stranded regions of RNAs (Kwok et al., 2015b). Two early genome-wide in vivo RNA structure probing methods, Structure-seq (Ding et al., 2014) and DMS-seq (Rouskin et al., 2014), were developed based on the RT stop method and successfully applied in bacteria, yeast and plants (Ding et al., 2014;Rouskin et al., 2014;Burkhardt et al., 2017;Ritchey et al., 2017;Su et al., 2018;Zhang et al., 2018;Tack et al., 2020). However, RT stop-based methods require appropriate DMS treatment with single-hit kinetics conditions, otherwise over-reaction by DMS causes highly skewed distribution of RT stops near the primerbinding sites. In addition, 3′ end RNA structure information is difficult to obtain due to short sequencing reads at the 3′ end of RNA. Another limitation is that degraded RNAs introduce false positive signals in the RT stop methods (Wu and Bartel, 2017;Zubradt et al., 2017;Yang et al., 2020a).
Alternatively, a new strategy called mutational mapping (MaP) was developed (Siegfried et al., 2014;Zubradt et al., 2017). Instead of using RT stops, the DMS-MaPseq method transfers DMS modifications on RNAs to cDNA mutations via a special RT enzyme called thermostable group II intron reverse transcriptase (TGIRT; Mohr et al., 2013;Zubradt et al., 2017). TGIRT mostly generates mismatch mutations on cDNA, with few insertions or deletions (indels). Therefore, TGIRT is preferred over other MaP methods that use an RT with SuperScript II (SSII) plus Mn 2+ , as this generates high numbers of indels (Zubradt et al., 2017). The high-fidelity of TGIRT endows the DMS-MaPseq method with single-nucleotide resolution and a high signal-to-noise ratio (Zubradt et al., 2017). Furthermore, DMS-MaPseq can be used to specifically investigate the RNA structure of low-abundance transcripts and isoforms (Zubradt et al., 2017;Guenther et al., 2018;Wang et al., 2018). Recently, two algorithms, DREEM (Detection of RNA folding Ensembles using Expectation-Maximization; Tomezsko et al., 2020) and DRACO (Deconvolution of RNA Alternative COnformations; Morandi et al., 2021), were developed to identify coexisting alternative RNA conformations of the same transcripts based on DMS-MaPseq data.
In plants, DMS-MaPseq was first used in a study of miRNA biogenesis. Wang et al. (2018) used target-specific DMS-MaPseq to profile the secondary structure of primary miRNAs (pri-miRNAs) between the wild-type (WT) and the mutant of the SWI/SNF chromatin remodeling factor CHROMATIN REMODELING 2 (CHR2)/BRAHMA in Arabidopsis thaliana. It demonstrated that CHR2 remodels the secondary structure of pri-miRNAs to impede miRNA biogenesis (Wang et al., 2018). Later, the authors optimized genome-wide DMS-MaPseq for Arabidopsis materials (Wang et al., 2019). In addition, a recent study performed DMS-MaPseq, as well as SHAPE-seq, to query the RNA structure of chloroplast transcripts and revealed RNA structure-mediated translational regulation of psbA and other plastid genes with weak Shine-Dalgarno sequences (Gawronski et al., 2021). However, the use of DMS-MaPseq in plants has been limited to Arabidopsis.
In this study, we optimized DMS-MaPseq for profiling the in vivo RNA secondary structure in rice (Oryza sativa; Figure 1). We optimized DMS treatment conditions for rice leaves and roots. Then we assessed the quality and reproducibility of our target-specific and genome-wide DMS-MaPseq data. To improve the sequencing depth and coverage for genome-wide DMS-MaPseq, we adopted a streptavidin-biotin depletion approach to reduce the abundance of highly expressed chloroplast transcripts during library construct. Then, we validated the feasibility of our rice DMS-MaPseq data for in vivo RNA secondary structure prediction. Finally, we used our RNA structure data to model the in vivo RNA secondary structure of a key regulator of rice phosphate (Pi) homeostasis. Altogether, the optimized DMS-MaPseq for rice could facilitate the study of RNA structure-mediated biological functions in crops.

Plant Material and Growth Conditions
Rice (Oryza sativa ssp. Japonica Nipponbare) seedlings were grown in a hydroponic nutrient system at 30°C under a 12 h light-12 h dark cycle. Three-week-old seedlings were harvested for DMS treatment.

DMS Treatment
DMS treatment was performed as described (Wang et al., 2018(Wang et al., , 2019. Three-week-old rice seedlings grown in hydroponic solution were collected and immersed in 20 ml of 1× DMS reaction buffer (40 mM HEPES pH 7.5, 100 mM KCl and 0.5 mM MgCl 2 ). Then, 200, 300, 400, 500, or 600 μl of DMS (MACKLIN, Cat#: D824267) was added to the reaction buffer to give a final concentration of 1, 1.5, 2, 2.5, and 3%, respectively. For DMS treatment control, leaves of 3-week-old Arabidopsis grown on the soils were collected and treated with 1% DMS. For the no DMS control, the same volume of DEPC-treated water was added into the reaction buffer. DMS treatment was performed for 15 min at 30°C with shaking at 250 rpm or under vacuum (approximately 12 psi) at room temperature without shaking. After the DMS treatment, 6 ml of β-mercaptoethanol (Sigma, Cat#: M6250) was added to a final concentration of 23% and incubated under vacuum for 5 min to quench the DMS reaction. Then, the samples were washed three times with DEPC-treated water, frozen in liquid nitrogen, and ground to fine powder.

RNA Extraction
The powdered rice materials (0.1 g) were mixed with 1 ml of Trizol reagent (invitrogen, Cat#: 15596018) and total RNA was extracted following the manufacturer's protocol. Denatured agarose gel electrophoresis was used to validate the integrity of total RNA.

Primer Extension Assays
The DMS treatment was validated by a primer extension assay as previously described (Wang et al., 2018(Wang et al., , 2019  Following phenol-chloroform extraction and ethanol precipitation, the cDNA was sizefractionated on a 10% urea-polyacrylamide gel. Next, the cDNA was transferred to a positively charged nylon membrane (GE Healthcare, Cat#: RPN303B) via a semi-dry blotter (Bio-rad). Immobilized cDNA was detected following the procedure of the chemiluminescent nucleic acid detection module kit (Thermo Fisher, Cat#: 89880) and the signal was collected with iBright1500 (Thermo Fisher).

Target-Specific DMS-MaPseq
Target-specific DMS-MaPseq was performed as described (Zubradt et al., 2017;Wang et al., 2018) with some modifications. Following DNase treatment, 3 μg of non-DMS-treated sample or 6 μg of DMS-treated sample was mixed with 0.5 μl of 10 μM genespecific RT primers mixture (Supplementary Table 1). The mixture was precipitated and re-suspended in a 10 μl of Tris-KCl solution (50 mM KCl and 10 mM Tris-HCl, pH 7.5). The solution was heated at 75°C for 3 min, followed by incubation at 57°C for 15 min. Then, 4 μl of 5× First-Strand buffer, 1 μl of 0.1 M DTT, 1 μl of SUPERase-In RNase inhibitor (Thermo Fisher), 1 μl of RNase-free H 2 O and 1 μl of TGIRT-III (InGex, Cat#: TGIRT50) were added to the solution. After incubation at room temperature for 30 min, 2 μl of 10 mM dNTPs was added and reverse transcription was conducted at 60°C for 2.5 h. Then, 2 μl of 2.5 M NaOH was added to stop the reaction and decay the RNA. The mixture was incubated at 95°C for 3 min and neutralized by adding HCl. Next, the cDNA was purified with RNAClean XP beads (Beckman). Then, the targets were amplified with KOD-FX hot-start DNA polymerase (Toyobo) using genespecific primers (Supplementary Table 1). PCR products were gel purified and normalized according to band intensity. The library was constructed with NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB). The libraries were quantified using Agilent TapeStation before sequencing by 2 × 250 bp paired-end reads on the Illumina Novaseq 6000 at Novogene. In vitro Transcription of Biotinylated Anti-chloroplast RNA Probes Chloroplast genes were amplified using listed primers containing T7 promoter sequence (Supplementary Table 1) and cloned into TA/Blunt-Zero vector (Vazyme). The resulting plasmids were used as the PCR template to amplify DNA templates for in vitro transcription. Then, PCR products were gel purified with a PCR purification kit (Vazyme). The reaction mixture including 250 ng of purified DNA templates, 1 μl of 100 mM ATP, 1 μl of 100 mM GTP, 1 μl of 100 mM CTP, 0.4 μl of 100 mM UTP, 0.8 μl of 50 mM biotin-16-UTP (Lucigen), 1.5 μl of T7 RNA Polymerase Mix (NEB) and 1 μl of SUPERase-In RNase inhibitor (Thermo Fisher) was incubated at 37°C overnight. Then, TURBO DNase (Thermo Fisher) was added to digest the DNA template at 37°C for 1 h. Finally, RNA probes were purified by RNAClean XP beads (Beckman) and mixed with a final concentration of 160 ng/μl anti-PSBA probes, 70 ng/μl anti-RBCL probes, 15 ng/ μl anti-PSAB probes, 15 ng/μl anti-PSAA probes, 10 ng/μl anti-PSBC probes, 10 ng/μl anti-PSBB probes, and 10 ng/μl anti-PSBD probes.

Chloroplast RNA Depletion
Total RNA was extracted from DMS-treated rice leaves and was treated with TURBO DNase (Thermo Fisher). One microgram of DNase-treated total RNA was mixed with 0, 0.25 or 0.5 μl anti-chloroplast RNA Probe mixture in a 20 μl hybridization reaction (50 mM Tris-HCl pH 7.5, 100 mM NaCl). The mixture was put in a thermocycler at 68°C for 5 min, then ramped down by −0.1°C/s to 22°C, and finally held at 22°C for 5 min. Then, 100 μl of Dynabeads MyOne Streptavidin C1 (Thermo Fisher, Cat#: 65001) was washed and re-suspended in 40 μl of 2× binding and washing buffer according to the manufacturer's manual. Biotinylated probe-target hybrids were immobilized by streptavidin beads twice. The depleted RNA solution was precipitated by ethanol and re-suspended in 10 μl of RNase-free H 2 O. Then, 0.75 μl of random primer and 0.75 μl of 10 mM dNTP were added to 10 μl of RNA. The mixture was heated to 65°C for 5 min and incubated on ice for at least 1 min. cDNA was synthesized by adding 4 μl of 5× First-strand buffer, 0.75 μl of 0.1 M DTT, 0.75 μl of SUPERase-In RNase inhibitor (Thermo Fisher) and 0.5 μl of TGIRT (InGex), and 2.5 μl of RNase-free H 2 O. The reverse transcription was performed at 25°C for 10 min, then 42°C for 30 min, then 60°C for 1.5 h. Next, 2 μl of 2.5 M NaOH was added to stop the reaction and decaying RNA. The expression of chloroplast RNAs were measured by RT-PCR. The intensity of PCR products was quantified by ImageJ. The primers used for RT-PCR were listed in Supplementary Table 1.

Sequencing Alignment and Analysis
Data analysis for target-specific or genome-wide DMS-MaPseq was performed as described (Wang et al., 2019;Tomezsko et al., 2020). Briefly, after read quality filtering with TrimGalore, 1 clean reads were mapped to the reference genome using TopHat 2 with parameter settings: --library-type fr-firststrand --no-novel-juncs -N 15 --read-gap-length 10 --read-edit-dist 15 --max-insertion-length 5 --max-deletionlength 5 -g 3. In that, the 10% mismatch tolerance setting (-N 15 for 150 nt sequencing reads) is based on a previous study (Zubradt et al., 2017). Next, uniquely mapped reads were extracted from the bam file using the Linux command grep with NH:I:1 tag. After discarding mismatches located within 3 nt of an indel, mutations and sequencing depth were counted from the uniquely mapped bam file. The DMS mutation signal was calculated for each adenine (A) and cytosine (C) nucleotide as mismatch/sequencing depth. The used Python scripts are described in the previous study (Wang et al., 2019).

In vivo RNA Structure Prediction
Based on the DMS mutation signal, the secondary structures were modeled by RNAStructure (Reuter and Mathews, 2010). 3 DMS signals were color coded on structure models using VARNA. 4

Graph Drawing
Graphs with dot plots (individual data points) were drawn using GraphPad Prism 8 5 , R, 6 or Adobe Illustrator CC.

Accession Codes
The GEO accession number of the DMS-MaPseq data in this study is GSE197245.

Optimization of DMS Treatment Conditions for Rice
DMS-MaPseq requires that DMS penetrates into cells to modify RNA in vivo. However, rice absorbs a lot of silicon from the soil and deposits it in the leaves, stem, and husks to form silica bodies . These silica bodies serve as a physiological barrier, which not only resists pathogen infection and lodging but also hinders chemical penetration (Tamai and Ma, 2003;. To optimize the DMS treatment conditions for rice tissues, we treated 3-week-old rice leaves with varying DMS concentrations (1, 1.5, 2, 2.5, and 3%, v/v) and different incubation conditions (30°C with shaking at 250 rpm for 15 min, or with vacuum for 15 min). We also treated 3-week-old Arabidopsis leaves as a reference.
As seen in the Arabidopsis sample, total RNAs of DMS-treated rice samples were slightly degraded while untreated RNA was intact (Figure 2A), consistent with previous finding that high DMS concentrations cause RNA degradation (Wang et al., 2019). Moreover, the extent of decay was correlated with the DMS concentration, in which RNA treated with higher concentrations of DMS appeared more severely decayed. This indicated DMS modifications on the RNAs. We also found that shaking-and vacuum-treated samples exhibited similar RNA degradation patterns in denaturing agarose gel electrophoresis (Figure 2A). Unlike Arabidopsis that is a taproot system dominated by the primary root, rice has a fibrous root system dominated by vast crown roots and lateral roots. Rice roots contain multiple layers of cortex cells from several to more than 10 layers, impeding the penetration of DMS into inner cells (Henry et al., 2017). To assess the efficacy of the DMS treatment conditions on rice roots, we treated roots of 3-week-old hydroponically cultured rice plants with different DMS concentrations (from 1 to 3%, v/v) and incubation conditions (shaking or vacuum for 15 min). Consistent with the leaf results, total RNAs of root samples were modified under our DMS treatment conditions ( Figure 2B). Of note, the treatments with greater than 2% DMS caused very severe RNA decay (Figure 2B), which may compromise subsequent profiling of in vivo RNA structure. This suggested the optimal DMS treatment varies for different tissues. The tested DMS treatments for rice roots are much harsher than the optimized condition for Arabidopsis roots (0.75% DMS and 1 min incubation; Tack et al., 2020), indicating DMS is more difficult to penetrate rice roots than Arabidopsis roots.
To further measure the extent of DMS modification of RNA, we performed primer extension assays for 18S ribosome RNA (rRNA) in untreated and DMS-treated leaf samples. The primer extension assay was modified from a previously published protocol (Wang et al., 2018), with some modifications. Briefly, we used a 5′ end biotinylated 18S rRNA-specific primer and chemiluminescent detection, instead of 32 P-labeled primer and autoradiography. Compared to the untreated samples, DMS treatments led to less full-length cDNA and more truncations (Figure 2C), implying DMS modifications of RNA.
Consistent with this, image quantification showed that higher DMS concentrations caused an increase in the ratios of truncations-to-full-length cDNA ( Figure 2C). Of note, the vacuum-treated samples had a moderately higher truncations-to-full-length cDNA ratio than the corresponding shaking-treated samples, suggesting that DMS penetrates into plant cells more efficiently under vacuum than shaking ( Figure 2C). Altogether, these results suggested that our DMS treatment conditions efficiently modified in vivo RNA. Higher DMS concentration treatment results in more DMS modifications on RNA, but also causes greater RNA decay that deteriorates the quality of DMS-MaPseq library (Figure 2). We recommended using 1-2% DMS with shaking or vacuum for rice leaves and 1% DMS with shaking for rice roots.

Assessment of Target-Specific DMS-MaPseq Data
DMS alkylates A and C located in single-stranded regions. Then, DMS lesions on RNA are decoded as mismatches on cDNA through TGIRT (Zubradt et al., 2017). To assess the quality of DMS-MaPseq data with different DMS treatment conditions, we amplified a region of 18S rRNA (93-445 nt) from untreated and treated samples, followed by target-specific DMS-MaPseq and bioinformatic analysis. We examined the enrichment of DMS-induced mismatches on nucleotides. Compared to the untreated sample, the percentage of mismatches located at A and C, but not guanosine (G) and thymidine (T), was dramatically increased in DMS-treated samples ( Figure 3A). This indicates high signal-to-noise ratios in our DMS-MaPseq data. Furthermore, higher DMS concentrations led to more mismatches in A and C, indicating a dosage-dependent effect of DMS on RNA modification ( Figure 3A). In addition, vacuum treatment generated more mismatches than shaking, consistent with the results of the primer extension assay (Figures 2C, 3A). We also found root sample with 1% DMS and 15 min shaking treatment exhibited a similar mismatch percentage on each nucleotide with leaf samples of Arabidopsis and rice, suggesting that 1% DMS with 15 min shaking is optimized to rice root tissue ( Figure 3A).
Then, to assess the fidelity of DMS-MaPseq data from the different DMS treatments, we conducted a correlation analysis among samples treated with different DMS concentrations. The results revealed excellent correlation in the DMS mutation signals (the ratios of mismatches to total reads) among the 1, 2, and 2.5% DMS-treated samples ( Figure 3B). We also compared DMS mutation signals between shaking and vacuum treatment and observed high correlation ( Figure 3B). In addition, the different DMS treatments resulted in highly similar DMS mutation signals distribution patterns along the tested region of 18S rRNA ( Figure 3C). These results showed high fidelity of the in vivo RNA structure data generated by various treatments, suggesting that the DMS-MaPseq method can tolerate high DMS concentrations.
To test whether our DMS-MaPseq data accurately profiled in vivo RNA secondary structure, we mapped our DMS mutation signals of 1% DMS treated leaf sample to the evolutionally conserved 18S rRNA secondary structure (Ding et al., 2014). The results revealed high consistency between our DMS mutation signals and the well-known 18S RNA secondary structure, indicating our DMS-MaPseq data are consistent with the in vivo RNA structure (Figures 3D,E).
In summary, these results demonstrated the high quality of our rice DMS-MaPseq data.

Optimization of Genome-Wide DMS-MaPseq for Rice Materials
Genome-wide DMS-MaPseq requires extensive sequencing depth and coverage to generate reliable global RNA structure information (Zubradt et al., 2017;Wang et al., 2019). However, for species with large genome sizes, such as crops, in-depth sequencing would be costly and a computational burden. One way to address this issue is to reduce the amount of highly expressed transcripts (Wang et al., 2019).
Since chloroplast transcripts are highly expressed and account for almost half of the total coding RNA in plant leaves, we adopted a streptavidin-biotin depletion approach to decrease the abundance of chloroplast transcripts from total RNA ( Figure 4A). Several highly expressed chloroplast transcripts were selected, which account for approximately 55% of the total amount of chloroplast RNAs. Next, we made biotinylated anti-chloroplast RNA probes through in vitro transcription with biotin-UTP. Denatured RNA gel images showed the high purity of these home-made biotinylated RNA probes ( Figure 4B). We annealed the biotinylated probes on DMS-treated total RNA and depleted the targeted chloroplast transcripts with streptavidin magnetic beads ( Figure 4A, see section "Materials and Methods"). RT-PCR results showed that, compared with the no-probe control, the amount of targeted chloroplast transcripts was decreased in probe-treated samples, suggesting that our depletion method was successful ( Figure 4C). The result also showed a dosage-dependent effect of antisense probes on RNA depletion ( Figure 4C). In addition, we observed the depletion efficiency of PSBC and PSAB is lower than other targeted chloroplast transcripts, indicating insufficient annealing between probes and targets ( Figure 4C). It might contribute from DMS-induced modifications at the Watson-Crick face of A and C that compromise the probe-target annealing, or from the strong intramolecular RNA structure in targets that impedes probe interaction.
Notably, this streptavidin-biotin depletion approach could be applied to reduce tissue-specific highly-expressed transcripts in different tissue samples, increasing the sequencing depth and coverage of genome-wide DMS-MaPseq.

Application of Genome-Wide DMS-MaPseq on Rice Materials
To apply genome-wide DMS-MaPseq on rice materials, we chose total RNAs from untreated and 1% DMS-treated leaf samples to construct genome-wide DMS-MaPseq libraries. rRNAs were depleted with Ribo-Zero Kit and the amount of chloroplast transcripts were reduced by the streptavidin-biotin depletion approach discussed above. Then, following a previously published protocol for Arabidopsis materials (Wang et al., 2019), we prepared the libraries for rice using a commercial RNA-seq library kit with some modifications in the RT step to be compatible with TGIRT (Figure 1). We generated one biological repeat for untreated and three biological repeats for DMS-treated samples. After sequencing and quality filtering, we mapped clean reads to the Oryza sativa L. ssp. Nipponbare reference genome (MSU Rice Genome Annotation Project Release 7) by TopHat with 10% mismatch tolerance (Zubradt et al., 2017). For both untreated and treated samples, most reads were mapped to the reference genome. Next, we assessed the quality of our genome-wide DMS-MaPseq data. Compared with the untreated sample, the increased mismatches were specific to A and C in DMS-treated samples ( Figure 5A). This is consistent with the mode of DMS modification. Moreover, the mismatch ratio of A was slightly higher than that of C, similar to the published results in human, yeast, and Arabidopsis (Zubradt et al., 2017;Wang et al., 2019). Next, we used Pearson's r value and Gini index to measure the reproducibility among three DMS-treated biological replicates and obtained a high r value and a small Gini index difference among repeats, indicating the strong reproducibility of our data ( Figure 5B). Therefore, we merged the DMS-MaPseq data of the three DMS-treated biological replicates for further analysis.
To validate the feasibility of our genome-wide DMS-MaPseq data for in vivo RNA secondary structure prediction, we used our DMS mutation signals as constraints to predict the in vivo RNA secondary structure of U1 small nuclear RNA (snRNA) with known RNA structure. The predicted RNA secondary structure exhibited a four-way junction, highly consistent with the reference structure (Krummel et al., 2011 ; Figures 5C,D). Altogether, these results demonstrated that genome-wide DMS-MaPseq could be applied to rice materials and produce high-quality in vivo RNA structure information.

In vivo RNA Secondary Structure Modeling
We applied our optimized DMS-MaPseq to investigate the regulatory functions of in vivo RNA secondary structure in rice. FIGURE 3 | The ratiometric DMS signal per position normalized to the highest mismatch/total in the displayed region, which was set arbitrarily to 1.0. In tested 18S rRNA region (from 98 to 390 nt), 78.57% (true positive) of As and Cs that showed high DMS mutation signal (defined as normalized DMS activity ≥ 0.1) in our DMS-MaPseq data corresponded to single-stranded regions in the phylogenetic structure, whereas 63.27% (true negative) of As and Cs that showed low DMS mutation signal (defined as normalized DMS activity ≤ 0.05) in our DMS-MaPseq data corresponded to base-paired regions in the phylogenetic structure. Of the 21.43% (false positive) nucleotides (defined as normalized DMS activity ≥ 0.01) that were annotated as base-paired in phylogenetic structure, 66.67% nucleotides were positioned either at the end of a helix or proximal to a bulge or loop, which were known be flexible. Corrected for these positions, the values in parentheses showed higher true positive and lower false positive percentages. miRNA-mediated cleavage participates in various aspects of developmental and stress responses by suppressing gene expression or translation in plants (Rogers and Chen, 2013;Li et al., 2017). miR399 was the first identified miRNA involved in stress responses in plants (Fujii et al., 2005;Bari et al., 2006;Chiou et al., 2006). It is a key regulator of inorganic phosphate (Pi) homeostasis and the phosphate-starvation response pathway (Fujii et al., 2005;Bari et al., 2006;Chiou et al., 2006). miR399 binds to the 5′UTR of PHO2 and reduces its expression through posttranscriptional miRNA-mediated cleavage of the PHO2 mRNA (Lin et al., 2008;Pant et al., 2008). PHO2 encodes a ubiquitinconjugating E2 enzyme, that participates in protein degradation of the Pi exporter PHO1 (Liu et al., 2012). The effect of mRNA structure on miRNA-mediated cleavage has been studied in human, C. elegant, Drosophila, and Arabidopsis (Ameres et al., 2007;Long et al., 2007;Yang et al., 2020a), but it is unclear how mRNA structure regulates miR399-mediated cleavage.
Detecting the effect of mRNA structure on miR399-mediated cleavage requires probing in vivo RNA structure of PHO2 before cleavage. Regular genome-wide RNA structure probing methods only provide population-average RNA structure information, resulting in that RNA structures of pre-cleaved and cleaved miRNA-target mRNAs are indistinguishable. However, target-specific DMS-MaPseq can probe isoform-specific RNA secondary structures (Zubradt et al., 2017). Therefore, we used target-specific DMS-MaPseq to profile the in vivo RNA secondary structure of the OsmiR399 target site and its flanking region of pre-cleaved OsPHO2 transcripts in rice. The DMS-MaPseq data showed no obvious difference in DMS mutation signals between the miRNA target site and its flanking region, consistent with a previous finding that miRNA target sites are not structurally accessible for binding the miRNAinduced silencing complex (Yang et al., 2020a; Figure 6A).
Next, we modeled the in vivo RNA secondary structure of this region based on our DMS-MaPseq data and observed a single-stranded structure immediately downstream of the miRNA target site (Figure 6B). Our predicted in vivo RNA secondary structure supported a published conclusion that the singlestranded structure downstream of miRNA target sites, named Target Adjacent nucleotide Motif (TAM), facilitates miRNA cleavage (Yang et al., 2020a). Our results suggested a regulatory function of mRNA structure on miR399-mediated Pi homeostasis, and also validated that our optimized DMS-MaPseq method could be used to investigate the biological function of in vivo RNA secondary structure in rice.

DISCUSSION
RNA structure is considered to be another layer of gene expression regulation, participating in various aspects of RNA metabolism (Bevilacqua et al., 2016;Vandivier et al., 2016;Yang et al., 2018;Zhu et al., 2021). However, compared to model species, studies of in vivo RNA structure-dependent biological functions in crops are rare (Deng et al., 2018;Yang et al., 2021).
There are some hurdles for RNA structure study in crops. First, the complexity of the cell wall and presence of multiple Frontiers in Plant Science | www.frontiersin.org cell layers hinder penetration of small chemicals into cells to react with RNA. For instance, rice deposits silica bodies on the leaf and stem surface to resist pathogen infection and abiotic stress , which also block penetration of small chemicals, thereby hindering chemical-modificationbased in vivo RNA structure probing. Furthermore, rice roots contain multiple layers of cortex cells (from several to more than 10 layers), while Arabidopsis roots contain only one cortical-cell layer (Henry et al., 2017). Multiple cell layers hinder chemical uptake in inner cells. To deal with this issue, we tested several DMS treatment conditions, including varying DMS concentrations and incubation conditions. We balanced the DMS-induced mutation ratio and RNA decay, and suggested that 1-2% DMS is suitable for rice leaf samples and 1% DMS is suitable for rice root samples (Figure 2). Our DMS treatment optimization for rice could serve as a reference for other crops, such as wheat (Triticum aestivum) and maize (Zea mays).
Secondly, due to their large genome size, crops require substantially more sequencing reads than Arabidopsis to achieve sufficient sequence depth and coverage for reliable genomewide RNA structure information. The genome size of rice is relatively small compared to other major cereal crops, but it is still approximately threefold larger than the Arabidopsis genome (Yu et al., 2002). The level of coverage needed for the DMS-MaPseq method, such as 20× mismatch coverage, dramatically increases sequencing costs. In addition, low-abundance transcripts are difficult to detect. To solve this problem, some studies used etiolated plants and mRNA enrichment to reduce the amount of highly-expressed chloroplast RNAs (Deng et al., 2018;Su et al., 2018). The downsides of this are that etiolated plants are under stress condition, and mRNA enrichment could miss RNAs lacking a poly(A) tail. Alternatively, we used a streptavidin-biotin depletion approach to specifically reduce the abundance of highly-expressed chloroplast transcripts from total RNA, allowing for greater sequencing depth and coverage of low-expressed transcripts (Figure 4). Another strategy is using the target-specific DMS-MaPseq method to specifically profile the RNA secondary structure of low-expressed genes. It is worth noting that targetspecific DMS-MaPseq can also detect in vivo RNA secondary structure from different isoforms ( Figure  6; Zubradt et al., 2017).
The first in vivo RNA structurome of rice was profiled by the Structure-seq (Deng et al., 2018). Both Structure-seq and DMS-MaPseq use DMS modification for RNA structure probing. DMS modifications on RNA are decoded through RT mutation in our optimized DMS-MaPseq method, instead of RT stop in Structure-seq. DMS-MaPseq increases RNA structure information content in sequencing data and reduces false-positive signals from unwanted RNA decay (Wang et al., 2021). However, RT mutation-based DMS-MaPseq method requires a greater sequencing depth to generate accurate RNA structurome (Zubradt et al., 2017;Wang et al., 2019). Together, our optimized DMS-MaPseq is complementary to the Structure-seq.
Studies have shown that RNA structure plays important roles in abiotic stress responses in plants (Anderson et al., 2018;Su et al., 2018;Chung et al., 2020;Kramer et al., 2020;Tack et al., 2020;Reis et al., 2021). However, the biological functions of RNA structure in plant nutrientdeficiency stress are still elusive. Pi is an essential nutrient for crop growth and production (Oldroyd and Leyser, 2020). Due to the low solubility and slow diffusion of Pi in soil, approximately 70% of global cultivated land suffers from Pi deficiency (Raghothama, 1999;Lopez-Arredondo et al., 2014;Paz-Ares et al., 2022). To sustain modern agriculture and global crop yield, it is of great importance to understand A B FIGURE 6 | Profiling RNA secondary structure of a OsmiR399 target site located at 5′UTR of OsPHO2. (A) Ratiometric DMS signals of each A and C were plotted along the tested 5′UTR of OsPHO2. The OsmiR399 target site was shown by a black line. (B) The RNA secondary structure of an OsmiR399 target site located at 5′UTR of OsPHO2 was predicted based on our DMS-MaPseq data. The ratiometric DMS signal per position normalized to the highest mismatch/total in the displayed region, which was set arbitrarily to 1.0.
plant Pi-starvation responses (PSR) and improve Pi utilization efficiency of crops. A recent report showed that the Pi-starvation induced long non-coding RNA cis-NAT PHO1;2 enhances the translation of the Pi exporter gene PHO1;2 via an internal RNA-RNA interaction and RNA structure change (Reis et al., 2021). Here, we used target-specific DMS-MaPseq to profile the in vivo RNA secondary structure of an OsmiR399-target region of OsPHO2, which is a key PSR gene involved in Pi transport (Lin et al., 2008;Pant et al., 2008;Liu et al., 2012). We found a single-stranded region downstream of the miR399 target site (Figure 6), which may facilitate miRNA-mediated cleavage (Figure 6). These findings shed light on the regulatory function of RNA structure in plant nutrient metabolism. As OsPHO2 contains five miR399 target sites, it would be interesting to investigate whether single-stranded RNA structure are exhibited in other miR399 target regions. Moreover, the miR399-PHO2 regulatory mechanism is conserved across angiosperms (Bari et al., 2006), whether such RNA structure exists beyond rice also be an attractive topic for future study.
In this study, we presented an optimized and powerful DMS-MaPseq method for studying the biological functions of RNA structure in rice. We hope that this method, together with other advanced RNA structure probing approaches, will promote RNA structure-guided molecular breeding and crop improvement.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI) BioProject database under accession number GSE197245.

AUTHOR CONTRIBUTIONS
ZW conceived and designed research. QJ, SH, and GW conducted experiments. QJ and LZ performed data analysis. ZW and QJ wrote the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank colleagues and all Wang lab members for suggestions and discussion during the course of this work. We thank Changhao Li for bioinformatics advice and Jiaying Zhu for helpful comments on the manuscript.


Exosomes facilitate intercellular communication between uterine perivascular adipose tissue and vascular smooth muscle cells in pregnant rats

Perivascular adipose tissue (PVAT) is distinct from other adipose depots as it has differential gene and protein profiles and vasoactive functions. We have shown that pregnancy affects the morphology of PVAT surrounding the uterine arteries (utPVAT) differentially than the morphology of non-perivascular reproductive adipose depots (i.e., periovarian adipose tissue, OVAT). Here, we hypothesized that pregnancy modifies the profile (size and molecular weight) of exosome-like extracellular vesicles released by utPVAT (Exo-utPVAT) compared to OVAT (Exo-OVAT) and that primary uterine vascular smooth muscle cells (utVSMCs) can internalize Exo-utPVAT. Our findings indicate that utPVAT from pregnant and non-pregnant rats secrete exosome-like vesicles. Exo-utPVAT from pregnant rats were smaller (i.e., molecular size) and heavier (i.e., molecular weight) than those from non-pregnant rats, while pregnancy did not affect the size of Exo-OVAT. Immunocytochemistry and confocal microscopy showed that primary utVSMCs internalized Exo-utPVAT (both tissues from the same pregnant rat) labeled by the lipophilic tracer DiO. Treatment of isolated uterine arteries with Exo-utPVAT did not affect relaxation responses to acetylcholine (ACh) in pregnant or non-pregnant rats. Collectively, these findings demonstrate a novel type of intercellular communication between Exo-utPVAT and utVSMCs and indicate pregnancy modulates the morphology and cargo of Exo-utPVAT.

INTRODUCTION
Perivascular adipose tissue (PVAT) contributes to vascular tone regulation via its ability to release products with paracrine vasoactive properties (1)(2)(3)(4). We recently showed that perivascular adipose tissue surrounding uterine arteries (utPVAT) contributed to increases in uterine artery blood flow in pregnant but not in nonpregnant rats (5). In addition, utPVAT-conditioned media reduced relaxation responses to acetylcholine (ACh) in isolated uterine arteries (5,6). These effects were explained in part by pregnancy-induced changes in responsiveness of the uterine vascular wall to utPVAT-derived factors (5).
Exosomes have a critical role in the paracrine actions of nonperivascular adipose tissues because they carry cargo that modifies epigenetic and signaling pathways in target cells (7)(8)(9). Recent studies reported that similar to nonperivascular adipose tissues (10)(11)(12), PVAT also secretes extracellular vesicles (13)(14)(15). Our previous studies have shown that pregnancy, which is a physiological state of dynamic changes and adaptations, affects PVAT morphology, genotypic profile, and vasodilatory properties (5,6). Interestingly, pregnancy differentially affected utPVAT morphology compared with periovarian adipose tissue (OVAT), which is a nonperivascular adipose tissue surrounding reproductive organs and thus in proximity with utPVAT.
Based on previous findings, we sought to investigate whether utPVAT secretes exosome-like extracellular vesicles, if pregnancy modulates the morphology of these vesicles and if such vesicles are internalized by adjacent uterine vascular smooth muscle cells (utVSMCs). Furthermore, we sought to determine if exosome-like extracellular vesicles released by utPVAT (Exo-utPVAT) could explain in part the effects of utPVAT-derived factors on uterine artery dilation that we had seen in our previous work and if the effects of pregnancy on Exo-utPVAT differ from those on vesicles released from nonperivascular adipose depots. We selected OVAT as a reference tissue to address the latter because of its localization (reproductive adipose depot). We hypothesized that pregnancy modifies the profile of exosome-like extracellular vesicles released by utPVAT compared with OVAT, and utVSMCs internalize exosomes released from its adjacent utPVAT.

Animals and Experimental Design
Every protocol was conducted in agreement with the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health and the ARRIVE guideline (16). These protocols were approved by the Institutional Animal Care and Use Committee of the University of North Texas Health Science Center (IACUC No. 2017-0042).
Male and virgin female Sprague-Dawley rats were purchased from ENVIGO (Indianapolis, IN). Upon arrival, the animals were acclimatized for a week before experiments began and housed in standard conditions as previously described (5). Female rats were double housed, whereas male rats were single housed. All animals were fed with standard chow and tap water ad libitum. Females were placed with males overnight for mating. Gestational day 1 (GD1; term = 22-23 days) was defined as the day we identified spermatozoa in vaginal smears (5,17). Nonpregnant rats were tested at diestrus.
Final experiments were performed in pregnant (n = 24) and nonpregnant rats (n = 16) between the ages of 18-24wk old. Male rats (15-26 wk of age) were used only for mating. Euthanasia and tissue harvest of all animals were performed on GD16 (5, 6) between 8:00 AM and 10:00 AM. PVAT tissue cultures were immediately established after tissue harvest, and extracellular vesicles were isolated the next day.

Tissue Harvest
Pregnant rats (GD16) and nonpregnant rats were euthanized using isoflurane overdose followed by cutting the diaphragm and removing their hearts. Uterine artery remodeling in rodents occurs between gestational days 15-18 (18, 19). Furthermore, our previous studies on the vasoactive effects of PVAT and the effects of pregnancy on PVAT protein/gene profile have been performed on GD16 (5,6). Thus, in the present study, all pregnant rats were tested on GD16. Fetuses were euthanized via decapitation. Upon removal of the surrounding mesometrium, ovarian adipose tissue (OVAT), and utPVAT (Fig. 1A), segmented uterine arteries were isolated and used immediately for ex vivo vascular reactivity experiments or processed for isolation of utVSMCs.

Isolation of Exosome-Like Extracellular Vesicles from OVAT and utPVAT Tissue Culture Media
We isolated and purified exosome-like extracellular vesicles from OVAT and utPVAT tissue culture media using ultrafiltration and ultracentrifugation (20). After OVAT and utPVAT were weighed, tissues were minced in a tissue culture hood and placed in six-well plates (Cat. No. 3516, Corning, NY) with 5 mL sterile exosome-free Dulbecco's modified Eagles medium (DMEM, Cat. No. 10567-014, Gibco, Thermo Fisher Scientific, Waltham, MA) supplemented with 100 U/mL penicillin and 100 mg/mL streptomycin (Pen-Strep; Cat. No. 15140-122, Gibco, Thermo Fisher Scientific) and 10% exosome-free fetal bovine serum (FBS; Cat. No. 26140095, Gibco, Thermo Fisher Scientific). DMEM and FBS were depleted of exosomes via ultracentrifugation. Minced adipose tissues were then cultured at 37 C in a humidified incubator (HeraCell 150i, Model No. 1335, Thermo Fisher Scientific) containing 5% CO 2 for 24 h. At 24 h, tissue culture media was collected and filtered through a 0.2 mm-sized 50-mL tube top vacuum filter system (Cat. No. 430320, Corning, NY) at room temperature to remove any contaminating microvesicles and cell debris (15). Filtrates were then centrifuged in a fixed angle rotor at 100,000 g for at least 2 h at 4 C (Sorvall WX þ Ultra Series with FiberLite; F50L-8x39 rotor, Thermo Fisher Scientific) to pellet exosomes. The supernatant (Exo-free supernatant) was removed and later used as a negative control for Western blot experiments. The exosome-like vesicle pellets were resuspended in 110 mL of 0.2 mm-filtered 1Â phosphate-buffered saline (PBS; 10Â; Cat. No. BP399, Thermo Fisher Scientific, Waltham, MA). Samples were stored at 4 C until needed for further experimentation. Exosome-like extracellular vesicles from OVAT or utPVAT were designated as exosome-like extracellular vesicles released by OVAT (Exo-OVAT) and Exo-utPVAT, respectively.

Characterization of Exo-utPVAT and Exo-OVAT
Size of adipose tissue-derived extracellular vesicles was detected using a particle size analyzer (Zetasizer Nano ZS system, Malvern Panalytical, Malvern, UK) according to the manufacturer's instructions. Data were captured and analyzed with Zetasizer Nano software (Malvern Panalytical).
In addition to vesicle size, Western blot analysis was used to confirm the characterization of exosome-like extracellular vesicles by measuring content of known extracellular vesicle-associated protein markers [tumor susceptibility gene 101, TGS101, protein related to exosome biogenesis; antiprogrammed cell death 6-interacting protein, Alix, protein associated with macrovesicle biogenesis (21)]. Exosome-free supernatants were used as negative controls. Purified exosomal samples were reconstituted in ice-cold T-PER Tissue Protein Extraction reagent, and the concentration of extracted proteins was determined with a modified bicinchoninic acid protein assay as in other studies (6,17,22). Proteins (30 mg) were resolved by electrophoresis on precast gradient 4%-20% SDS-PAGE gels (Cat. No. 4561094, Bio-Rad, Hercules, CA).
Molecular masses of Exo-utPVAT and Exo-OVAT were determined using AKTA purifier (GE Healthcare) fast protein liquid chromatography (FPLC) system by monitoring the 280nm wavelength. Purified Exo-utPVAT and Exo-OVAT were suspended in 1Â PBS and then injected and passed through a Superose 6 (increase 10/300) column (fractionation range of 5,000-5,000,000 Da molecular mass), and 1-mL fractions were collected using an autofractionator. Known molecular mass standards (Cat. No. 28403842, GE Healthcare, Chicago, IL) were run on the FPLC systems to estimate the molecular mass of exosomes based on elution profile. Two independent samples per group were analyzed for each tissue. To ensure that the FPLC outcomes were specific to exosomes, we used as controls exosome-free supernatants generated during the exosome isolation process (ultracentrifugation). Two controls (control 1 and control 2) were used and represent two independent biological samples.

Isolation of Uterine Artery Vascular Smooth Muscle Cells
UtVSMCs from pregnant rats (GD16) were isolated from the main uterine arteries using an enzymatic digestion method (24). Cell pellets were resuspended in DMEM (with 10% FBS and 100 U/mL penicillin, 100 mg/mL streptomycin), and cell solutions were distributed in culture dishes and kept in a humidified incubator at 37 C (5% CO 2 ). Cell culture media was replaced every 48 h. Purity of utVSMCs was tested by immunocytochemistry for positive staining of VSMC marker smooth muscle actin (anti-a-smooth muscle actin, 1:400, Sigma Cat. No. A2547; secondary: anti-mouse Alexa 488, 1:500, Abcam, Cat. No. ab150113). Anti platelet endothelial cell adhesion molecule 1 (PECAM 1; CD31, 1:3,200, Cell Signaling, Cat. No. 3528S) was used as a negative control. Cultured utVSMCs were used at passage 2 to avoid genotypic and phenotypic drifts.

Coculture of Exo-utPVAT and utVSMCs
utVSMCs were plated into 24-well plates (Cat. No. 3526, Corning) containing cover glasses (Cat. No. 12-545-80, Thermo Fisher Scientific) at a seeding density of 40,000 viable cells/mL (20,000 cells/well, 0.5 mL growth media/ well). utVSMCs were allowed to incubate overnight for adherence and were then treated for 3 h (37 C, 5% CO 2 , and 100% relative humidity) with Exo-utPVAT (170 mg/well) stained with lipophilic tracer DiO cell-labeling solution (Invitrogen, V22886) as described in Horibe et al. (25) or Exofree supernatant (170 mg/well) stained with DiO. Exo-utPVAT and Exo-free supernatant concentrations were determined as the concentration of protein present in the samples, measured with NanoDrop (Thermo Scientific NanoDrop 1000, Thermo Fisher Scientific) and analyzed with its corresponding NanoDrop software v 3.8.1. Exo-utPVAT and Exo-free supernatant were incubated with 5 mL of 1 mM DiO for 20 min at 37 C. Stained solutions were separated from the remaining DiO by filtering through an exosome spin column (MW 3,000, Cat. No. 4484449, Invitrogen) at 750 g for 2 min at room temperature (25). Imaging controls were as follows: 1) untreated utVSMCs to determine autofluorescence of isolated cells, 2) 170 mg/well of unlabeled Exo-utPVAT to determine autofluorescence, and 3) DiO dye alone (5 mL of 1 mM DiO into 55 mL PBS) to ensure effective functioning of spin column removal of excess DiO. After treatment, utVSMCs were washed with PBS, fixed with 4% paraformaldehyde (PFA) for 10 min, permeabilized with 0.1% Triton X-100 for 5 min, stained with (1:1,000) Alexa Fluor 555 phalloidin (Cat. No. A34055, Invitrogen) for 30 min, and mounted with ProLong Diamond Antifade Mountant with DAPI (Cat. No. P36962, Thermo Fisher Scientific) and allowed to cure for 24 h at 4 C as described in Horibe et al. (25). All reagents and procedures were in room temperature until slides were placed at 4 C. Cells were washed with PBS (37 C) between steps. For these studies, uterine Exo-utPVAT and utVSMCs were isolated from the same timed-pregnant rat for each experiment. Two independent experiments (2 biological replicates) were performed.
Internalization of Exo-utPVAT into utVSMCs was visualized using a confocal microscope (LSM 880, Zeiss, Germany). Images were captured at Â20. Laser energy required to image DiO-labeled Exo-utPVAT-treated utVSMCs was the energy used for acquiring images for the rest of the samples from the same animal on the same experimental day. Confocal images were processed using ZEN imaging software (Zeiss, Germany). To confirm the uptake of exosomes into utVSMCs, a super-resolution Z-stack confocal image was generated from DiO-labeled utVSMCs treated with Exo-utPVAT.  100% relative humidity. After 3 h, uterine arteries were mounted onto wire myograph chambers (Danish Myo technology, Aarhus, Denmark), baseline tension was determined for individual segments, and vascular smooth muscle and endothelial integrity were assessed as in Ref. 5

and Ref. 17.
To determine whether Exo-utPVAT contributes to the effects of utPVAT on ACh-induced relaxation in uterine arteries as in Osikoya et al. (5, 6), we performed concentration-response curves to ACh (10 À9 -3 Â 10 À5 M) in uterine arteries contracted with phenylephrine (1 μM). For these experiments, Exo-utPVAT and uterine arteries were from different animals, because uterine artery reactivity requires immediate use of the arteries, whereas isolation of Exo-utPVAT requires at least 24 h. Uterine arteries from pregnant rats (GD16) were incubated with Exo-utPVAT from gestational age-matched pregnant rats, and uterine arteries from nonpregnant rats at diestrus were incubated with Exo-utPVAT from age-matched nonpregnant rats at diestrus.

Data Analysis
The number of animals necessary to complete the characterization of PVAT-derived exosomes and reactivity experiments was estimated using power analysis based on preliminary data and previous studies (5,6,27). These numbers were selected with a goal to achieve a power of 0.80-0.85 with a probability of a type I error of 0.05. Concentrationresponse curves were fit and analyzed using robust nonlinear regression analysis as we previously described (5). To quantify vascular responses to ACh, we calculated EC 50 (expressed as pEC50: negative logarithm of EC 50 ) for each sigmoidal curve, and group means were compared using an unpaired t test. Unpaired t test with Welch's correction and Mann-Whitney test were used to compare the size of Exo-utPVAT and Exo-OVAT, respectively, between pregnant and nonpregnant rats. All statistical analyses were conducted with GraphPrism software (v. 9.4, San Diego, CA). Exact P values are reported.

RESULTS
Particle size determination revealed that Exo-utPVAT had a median diameter ranging between 47 and 99 nm (Fig. 1, B and C), in line with the size range previously reported for exosomes (28). Exo-utPVAT from pregnant rats were smaller compared with Exo-utPVAT from nonpregnant rats (Fig. 1B), but this difference was not observed in Exo-OVAT (Fig. 1C). Figure 1, D and E, shows representative graphical displays of intensity particle size distribution (PSD) with logarithmic xaxis and the linear y-axis from measurements of the size of Exo-utPVAT (Fig. 1D) and Exo-OVAT (Fig. 1E) derived from pregnant rats.
Vesicles isolated and purified from OVAT and utPVAT expressed exosome-associated proteins TSG101 and Alix (Fig. 2, A and B). With FPLC, two major protein peaks were recognized in chromatograms from all exosome samples Exo-utPVAT, exosome-like extracellular vesicles derived from uterine perivascular adipose tissue; utPVAT, uterine perivascular adipose tissue; utVSMCs, uterine vascular smooth muscle cells. (Fig. 2, C and D). Control samples (Fig. 2C) did not show these peaks and most of the peaks for the control samples were concentrated at late elution times beyond 40 min suggesting small molecular mass proteins/complexes. The two major peaks noted at exosome samples were with elution times of 26.3 and 34.4 min in the chromatogram. These elution times were used to estimate the molecular mass of exosomes based on the known molecular mass of standard proteins that were ran on the FPLC system. Based on this extrapolation, molecular mass of Exo-utPVAT from pregnant rats was $95-100 and $2,000 kDa, whereas molecular mass of Exo-utPVAT samples from nonpregnant rats was $97, $250-310, and $2,000 kDa. The molecular masses of Exo-OVAT were $80-100 and $2,000 kDa in samples from pregnant rats and $80 kDa, $270 kDa, and $2,000 kDa in samples from nonpregnant rats. In Exo-utPVAT samples, the protein peak corresponding to $2,000 kDa was $9.2Â greater in samples from pregnant compared with nonpregnant rats, whereas in Exo-OVAT, the same peak was 13.8Â greater in samples from pregnant compared with nonpregnant rats.

DISCUSSION
Perivascular adipose tissue is a metabolically active adipose depot that secretes adipocytokines with paracrine actions on the adjacent vascular wall (5,(29)(30)(31). A recent study by Li et al. (15) demonstrated that mouse mesenteric PVAT secretes extracellular vesicles that carry their cargo into VSMCs. This investigation added extracellular vesicles in the list of utPVAT-derived products and opened new horizons in the field of PVAT biology (15). The present study provides new knowledge about intercellular communication between PVAT and maternal vessels in the uterine circulation, which may be relevant for the establishment of blood flow to the growing fetus during pregnancy.
Previous studies have shown increases in maternal nonperivascular adipose tissue mass (27,32), adipocyte cell number and size (27,33), fat accretion, and lipid metabolism (34) in response to pregnancy in humans and in rodent experimental models. With regard to PVAT, we reported that pregnancy differentially modified the cellularity of utPVAT and nonperivascular adipose tissue (i.e., OVAT) such that the adipocytes of OVAT were larger in pregnant compared with nonpregnant rats on gestational day 16 (third trimester of rat pregnancy; 5), whereas there were no differences in the size of utPVAT adipocytes between groups (5). Furthermore, we found that pregnancy downregulated gene expression of metabolic adipokines in utPVAT (5) but not in aortic PVAT (5) or OVAT (unpublished observations). These results suggest that utPVAT has a different phenotypic and gene expression response to pregnancy compared with nonperivascular adipose depots.
The differential effects of pregnancy on types of adipose tissue depots were further demonstrated in the present study because pregnancy reduced the size of exosomes from utPVAT but had no effect on the size of exosomes from OVAT. Size separation of purified exosomes from utPVAT and OVAT revealed that the relative abundance of vesicles with molecular mass of $2,000 kDa was greater in the samples from pregnant animals. In addition, Exo-utPVAT and Exo-OVAT samples from pregnant rats did not have vesicles of $250-310 kDa. These data indicate that pregnancy induces differential effects on phenotypic and gene expression of maternal utPVAT and non-PVAT tissues and also modulates the morphology and possibly the content of utPVAT-derived exosome-like extracellular vesicles. To the best of our knowledge, this is the first report of a physiological condition that reduces the size of exosomes but increases their molecular mass.
When we investigated the communication between uterine Exo-utPVAT and utVSMCs, we discovered that primary utVSMCs internalized Exo-utPVAT. In our previous work, we reported that utPVAT-derived factors suppressed uterine artery relaxation responses to ACh and this effect was specific to pregnancy (5). Here, we sought to investigate whether these effects were attributed to Exo-utPVAT but found that the presence of Exo-utPVAT did not have an effect on uterine artery relaxation. Since we have no evidence to demonstrate uptake of Exo-utPVAT from the isolated uterine arteries, we should consider the possibility that Exo-utPVAT was degraded in the myograph chamber during the incubation period. However, our cell culture experiments suggest that utVSMCs have internalized exosomes after 3 h of coculture with them. This time period is the same as the one in the vascular reactivity studies. Others have shown internalization of Exo-utPVAT by primary VSMCs even after 18 h of coculture (15). Our findings are interpreted with caution in the context of pregnancy because although primary utVSMCs were used near isolation (passage 2), potential genotypic or phenotypic drifts should not be disregarded.
During pregnancy, the uterine vascular bed becomes a lowresistance shunt that allows increases in uterine blood flow to the growing fetus, the newly developed placenta, and the myometrium (35). These changes are brought by an increase in vascular dilatory responses and extensive regional remodeling (36). Li et al. (15) demonstrated that extracellular vesicles from mesenteric adipose tissue and isolated mesenteric adipocytes carried miRNAs that affected VSMC proliferation and migration in vitro. Thus, the involvement of utPVAT-derived exosome-like vesicles in uterine vascular remodeling during pregnancy is an exciting possibility.
In conclusion, utPVAT secretes exosome-like vesicles, which are internalized by adjacent utVSMCs. To the best of our knowledge, this is the first study to demonstrate exosomal interaction between tunica media and tunica adiposa in uterine arteries. Our findings are timely and highly relevant to recent advancements in the field of uterine vascular adaptations to pregnancy. Allerkamp et al. (37) used mathematical modeling to demonstrate distinct functional and structural profiles of uterine arteries at the pregnant compared with the nonpregnant states. Our findings suggest that the interaction between utPVAT and utVSMCs via exosomes should be considered in pregnancy-specific multiscale models that aim to improve understanding of mechanisms underlying uterine artery adaptive responses to pregnancy.


Emergence of Asian endemic begomoviruses as a pandemic threat

Plant viruses are responsible for the most devastating and commercially significant plant diseases, especially in tropical and subtropical regions. The genus begomovirus is the largest one in the family Geminiviridae, with a single-stranded DNA genome, either monopartite or bipartite. Begomoviruses are transmitted by insect vectors, such as Bemisia tabaci. Begomoviruses are the major causative agents of diseases in agriculture globally. Because of their diversity and mode of evolution, they are thought to be geographic specific. The emerging begomoviruses are of serious concern due to their increasing host range and geographical expansion. Several begomoviruses of Asiatic origin have been reported in Europe, causing massive economic losses; insect-borne transmission of viruses is a critical factor in virus outbreaks in new geographical regions. This review highlights crucial information regarding Asia’s four emerging and highly destructive begomoviruses. We also provided information regarding several less common but still potentially important pathogens of different crops. This information will aid possible direction of future studies in adopting preventive measures to combat these emerging viruses.

Plant viruses are responsible for the most devastating and commercially significant plant diseases, especially in tropical and subtropical regions. The genus begomovirus is the largest one in the family Geminiviridae, with a single-stranded DNA genome, either monopartite or bipartite. Begomoviruses are transmitted by insect vectors, such as Bemisia tabaci. Begomoviruses are the major causative agents of diseases in agriculture globally. Because of their diversity and mode of evolution, they are thought to be geographic specific. The emerging begomoviruses are of serious concern due to their increasing host range and geographical expansion. Several begomoviruses of Asiatic origin have been reported in Europe, causing massive economic losses; insect-borne transmission of viruses is a critical factor in virus outbreaks in new geographical regions. This review highlights crucial information regarding Asia's four emerging and highly destructive begomoviruses. We also provided information regarding several less common but still potentially important pathogens of different crops. This information will aid possible direction of future studies in adopting preventive measures to combat these emerging viruses.

Introduction
Geminiviruses are the most destructive group of viruses infecting agricultural crops in the tropical and subtropical regions of the world (Laufs et al., 1995;Mansoor et al., 2003;Marwal et al., 2012;Inoue-Nagata et al., 2016). They are also among the largest group of viruses transmitted by whiteflies (Bemisia tabaci; Bedford et al., 1994;Inoue-Nagata et al., 2016). The members of family Geminiviridae have circular single-stranded DNA (ssDNA) genome encapsulated in a twinned icosahedral coat protein.
In the nineteenth century, symptoms of geminivirus infection were first observed in plants in the tropical and subtropical regions (Moffat, 1999;Inoue-Nagata et al., 2016;Orfanidou et al., 2019). They seriously threatened crop production and agricultural fields (Moffat, 1999). However, in this century, multiple newly emerging geminiviruses have destroyed crops, such as legumes, cotton, and tomatoes (Cohen and Antignus, 1994;Varma and Malathi, 2003). For example, in Africa, the geminiviruses on cassava plants caused an economic loss of approximately US$ 1,300-2,300 million (Hillocks et al., 2002), and in the United States, the tomato yellow leaf curl disease (TYLCD) resulted in an average loss of US$ 140 million every year (Morales and Anderson, 2001).

Geminivirus strains and species
The International Committee on Taxonomy of Viruses group working on the family Geminiviridae proposed a new standard for species identification (Francki et al., 2012). Based on this classification, more than 500 geminiviruses have been distinguished and characterized into 14 distinct genera (Chiumenti et al., 2021;Lal et al., 2021). The order of various genera depends on the genome, host range, and viral vectors (Fauquet et al., 2003;Varsani et al., 2014). Currently recognized genera of Geminiviridae, their properties, and type members are described in Table 1. Figure 1 depicts the species distribution in each genus.

General characteristics of begomoviruses
The largest and most important genus of Geminiviridae is Begomovirus, transmitted by the whitefly (B. tabaci): (Khan et al., 2012;Yaqoob et al., 2020;Vivek et al., 2021). It is proposed that single stranded DNA viruses developed from prokaryotic episomal DNA and recombine with the target plant's genome to acquire new characteristics (Rojas et al., 2005). The size of the begomovirus genome is between 2,600 and 2,800 nucleotides (Ho et al., 2014). They replicate by rolling circle amplification (RCA) using a double-stranded DNA (dsDNA) intermediate (Yadava et al., 2010). Begomoviruses are classified into four categories based on their genomic structure and evolutionary relationships: New World (NW), Old World (OW), Sweepoviruses, and Legumoviruses (Fiallo-Olivé et al., 2021).
In OW and NW, pepper, tomato, and cucurbit plants are common hosts for begomoviruses (Seal et al., 2006). Most NW begomoviruses have a bipartite genome made up of DNA-A and DNA-B components. DNA-A found both in monopartite and bipartite begomoviruses, they are similar from genomic structure and arrangement point of view and have similar proteins encoding genes like C1/AC1, C2/AC2, C3/AC3, C4/ AC4, C5/AC5 and V1/AV1 (Gutierrez, 1999;Harrison and Robinson, 1999;Zhao et al., 2022). They also have V2 protein which is absent in OW bipartite begomoviruses. But NW bipartite has AV2 genes analogue to this V2 gene of monopartite (Sudarshana et al., 1998;Rojas et al., 2005). V3 protein is also present, essential for complete viral infection and works as an RNA silencing suppressor (Gong et al., 2021). DNA-B is found in bipartite genome. It interprets two proteins for cell-to-cell movement . Two proteins are BV1 and BC1. It also has Intergenic region (IR). IR has a region named as conserved region in which about 200 nucleotides are there which carry almost sequence similarity of more than 85%. This has nonanucleotide sequence which is found in all geminiviruses, and it is their speciality (Roshan et al., 2017). Begomoviruses in the OW can be both monopartite or bipartite and are associated with betasatellites, alphasatellites and deltasatellites (Fiallo-Olivé et al., 2021). Alphasatellites, which are mostly found in monopartite OW begomoviruses, have a genome that encodes a replication-associated protein that is required for replication. Betasatellites, which are found in association with many monopartite OW begomoviruses, are required for the onset of common disease symptoms (Zhou, 2013;Gnanasekaran et al., 2019). The betasatellite genome encoded C1 protein plays vital functions in symptom induction and the inhibition of transcriptional and post-transcriptional gene silencing (Li et al., 2018). Deltasatellites, unlike betasatellites and alphasatellites, do not encode any gene. All deltasatellites share several genomic features, including a stemloop containing the conserved begomovirus nonanucleotide TAATATTAC, a putative secondary stem-loop structure located near begomovirus iteron-like sequences, a short region with the sequence identity of the betasatellite conserved region, and an A-rich region (Fiallo-Olivé et al., 2012). NW bipartite begomoviruses have also potential to interact with betasatellite. Cotton leaf curl Multan betasatellite (CLCuMuB) associated with a serious disease of cotton is capable to interact with NW begomovirus. NW Cabbage leaf curl virus of interact with CLCuMuB to enhanced symptoms development in Nicotiana Benthamiana (Nawaz-ul-Rehman and Fauquet, 2009). Figure 2 depicts the genomic structure of the begomovirus.

Host range and symptoms
Begomovirus infects a wide range of dicot and monocot plants and causes severe loss, exhibiting symptoms such as stunting, deformed growth, and decreased seed production in infected Frontiers in Plant Science 03 frontiersin.org plants (Dasgupta et al., 2003); leaf crumpling, curling, distortion, golden-light green-yellow mosaic/mottle, interveinal yellowing, yellow spots, vein swelling, purpling, and yellowing in dicotyledonous plants (Rishi, 2009). Table 2 describes the host range and symptoms of some important begomoviruses.

Geminiviruses are unequally distributed across the continents
The Geminiviruses are present on all the continents except Antarctica. Maximum diversity of geminiviruses exists in Southeast Aisa. In contrast to several tropical and subtropical regions, the genetic variability of Geminiviridae in Europe is low (Bendahmane et al., 1995;Lindsten and Lindsten, 1999). Wheat dwarf virus, a solitary Mastrevirus localized to northern and central Europe, causes scattered yield loss in wheat (Bendahmane et al., 1995;Lindsten and Lindsten, 1999). All other geminiviruses found in Europe belong to the genus Begomovirus and are distinguished based on their genetic characteristics and region of origin (Briddon, 2002). Table 3 lists some crucial begomoviruses on the European and Mediterranean Plant Protection Organization (EPPO) alert list. Tomato yellow leaf curl virus (TYLCV) and Tomato leaf curl New Delhi virus (ToLCNDV) are prominent begomoviruses in Europe, which were previously limited to Mediterranean and Indian subcontinent respectively.

Tomato yellow leaf curl virus
Solanum lycopersicum (Tomato) is the host of a large number of viruses. These viruses cause significant losses in fruit production and quality (Hanssen et al., 2010). TYLCV is one of the important begomoviruses and it currently ranks third on the list of crucial plant viruses worldwide (Scholthof et al., 2011;Rybicki, 2015). The first case of TYLCD was reported in the late 1930s in Jordan Valley, Israel, and TYLCV was officially declared the virus of this disease after the 1960s . After that, TYLCV expanded uncontrolled across the Mediterranean basin and most tropical and subtropical areas of the globe, and it is now recognized as one of the world's most damaging virus for tomato . The disease continues to spread to new areas, with significant outbreaks in Trinidad and Tobago (Chinnaraja et al., 2016) and Costa Rica (Barboza et al., 2014). Species distribution of geminivirus. Number of species in each genus. It can be clearly seen that the genus Begomovirus has the highest number of species among the family Geminiviridae.
TYLCV and 12 TYLCV-like viruses are members of a virus complex that causes TYLCD https://talk.ictvonline.org/ taxonomy/. The significant symptoms of TYLCD in tomatoes are leaf discoloring, curling, and plant stunting. Furthermore, during extreme infection, flowers, and fruit were abscised, followed by a full reduction in plant growth. (Yan et al., 2021).
Although TYLCD spread worldwide, however, just two strains namely Mild (TYLCV-Mild) and Israel (TYLCV-IL) are truly international TYLCD-causing entities (Navas-Castillo et al., 2011). Numerous begomoviruses associated with TYLCD were found only in specific geographic areas, such as Tomato yellow leaf curl Sardinia virus and Tomato yellow leaf curl China virus, which was observed only in the Mediterranean and China (Navas-Castillo et al., 2011;Pan et al., 2012). The trading of plant materials is a significant contributor to the global spread of TYLCD (Seal et al., 2006). This global distribution is also associated with the worldwide increase in the population of insect vectors and the rapid evolution of virus variants .
TYLCV is potentially distinctive among begomoviruses in that it can replicate within B. tabaci (Czosnek et al., 2017;He et al., 2020) and seeds (Kil et al., 2016;Pérez-Padilla et al., 2020) Such feature may have contributed significantly to its global distribution and is spreading to new areas in the Indian and Pacific regions, including New Caledonia, Australia, and Mauritius. A study was conducted in which they performed a temporal-scaled, phylogeographic analysis of all publicly released TYLCV complete genome sequences and 70 new genomic sequences from Australia, Iran, and Mauritius. This indicated that epidemics in Australia and China were probably the product of multiple individual viral introductions from the East Asian region surrounding Japan and Korea. The New Caledonian epidemic was due to a variant from the Western Mediterranean region, and the Mauritian outbreak by a variant from the adjacent Island of Réunion. This study also revealed that the movement of TYLCV to East Asia has stopped temporarily, while that to America and Australia continues . The TYLCV spread in new countries after 2010 can be seen in Figure 3.

Tomato leaf curl New Delhi virus
One of the most noticeable disease on tomatoes and other vegetables in Asia and now in Europe is the ToLCNDV (Briddon, 2002;Hussain et al., 2004). Unlike TYLCV, which is monopartite in nature, the ToLCNDV is a bipartite begomovirus that causes damage to cultivated plant species of the Solanaceae family, including tomato, Solanum tuberosum (potato), Capsicum frutescens (chili), Piper nigrum (pepper), and Solanum melongena (eggplant). In its early days, its outbreak was restricted to Asian countries . Recently, ToLCNDV has increased its host plant range, including Euphorbiaceae, Cucurbitaceae, Fabaceae, and Malvaceae .
ToLCNDV was first reported in the western Mediterranean region in 2012, infecting Cucurbita pepo L (zucchini squash), Cucumis melo L (melon), and Cucumis sativus L (cucumber) crops in Southern Spain (Juárez et al., 2014;San Ambrosio and Fernández, 2015). Recent ToLCNDV epidemics in the Mediterranean basin have been associated with the emergence of a new strain, ToLCNDV-ES, which primarily attacks cucurbits, such as cucumber, zucchini, and melon. The molecular characterization of this virus revealed that it had 98% similarity to partial Coat protein gene sequences from isolates of ToLCNDV, which infects cucumber plants in India (Desbiez et al., 2020). The emergence of ToLCNDV in the Mediterranean Basin poses a new threat to commercially vital cucurbit crops and tomato production, as Spain is a significant producer of cucurbits worldwide and Europe's first exporting country (Juárez et al., 2014;Ruiz et al., 2017).
ToLCNDV was then found to infect zucchini crops in Greece in 2018 (Orfanidou et al., 2019), adding to its geographical diversity. Later, ToLCNDV expanded to Tunisia, Italy, Morocco, and Greece (Panno et al., 2016;EFSA Panel on Plant Health et al., 2020). Furthermore, ToLCNDV has recently been identified in cucurbit plants in Portugal and Estonia and in members of the Solanaceae family in Italy, indicating that it is rapidly spreading throughout Europe (Parrella et al., 2018(Parrella et al., , 2020. The expansion of ToLCNDV can be seen in the world map in Figure 4.

Cotton leaf curl disease (CLCuD) in the Indian sub-continent is a major limiting factor on cotton crop
Begomoviruses pose a serious threat to all major crops in Asia, e.g., CLCuD has been observed to be caused by various begomoviruses, i.e., Cotton leaf curl Multan virus (CLCuMuV), Cotton leaf curl Kokhran virus (CLCuKoV), Cotton leaf curl Gezira virus (CLCuGV), Cotton leaf curl Burewala virus etc. (Briddon and Markham, 2000;Farooq et al., 2011). These cottoninfecting viruses may work in association with betasatellites to cause infection, i.e., CLCuMuB in Asia and cotton leaf curl Gezira betasatellite (CLCuGB) in Africa (Tahir et al., 2011).
Asia and the Mediterranean regions suffer from leaf curling disease due to the begomoviruses and satellite DNA that make up a complex system (Nawaz-ul-Rehman and Fauquet, 2009).

Cassava mosaic disease can cause famine in Africa
Cassava mosaic disease (CMD)-related viruses are another example of essential begomoviruses in Africa, Sri Lanka and India.
Currently, the Manihot esculenta (cassava) crop is highly influenced by CMD due to the activity of at least nine bipartite local African begomoviruses, called cassava mosaic geminiviruses (CMGs;De Bruyn et al., 2016). This high level of diversification indicated several local virus introduction events in this crop. Madagascar, a restricted geographical area, can serve as an example of various introduction events and interactions of this group of CMGs (De Bruyn et al., 2016). Since the late 1980s, CMGs have accelerated the spread of CMD outbreaks, damaging cassava crops in almost 12 African countries and the southwest Indian Ocean islands (Legg et al., 2011).
African cassava mosaic virus (ACMV) is the most common cassava mosaic begomovirus in Nigeria. Which spread by contaminated stem cuttings and whitefly vectors (Eni et al., 2021). Although CMD has been prevalent in Africa since the nineteenth century, it was not known in Southeast Asia until May 2015, when the first CMD epidemic was discovered in Cambodia . The viral species causing CMD in Southeast Asia is Sri Lankan cassava mosaic virus (SLCMV), a bipartite begomovirus that is common in Cambodia, Vietnam, Thailand, and south China (Siriwan et al., 2020).

A global movement of begomoviruses that pose a serious threat
Emerging begomoviruses pose a serious threat to agricultural production worldwide. These viruses can be new (i.e., previously unknown) or already known; however, they share the common characteristic of occupying and spreading within new niches (Elena et al., 2014;Ertunc, 2020). The factors driving the emergence of plant viruses include genetic variability and the global movement of plant materials, i.e., nursery plantations (Rojas and Gilbertson, 2008;Elena et al., 2014). By applying modern science, improved cultivated crop varieties ("cultivars") can significantly enhance crop yields and improve agricultural productivity. Improved cultivars can lead to greater production combined with other modern inputs and good crop management practices. This has increased the demand for improved seeds worldwide (Baker and Smith, 1966;Johansen et al., 1994;Evenson and Gollin, 2003). Therefore, the study of seed transmission of begomoviruses is crucial to limit the spread of viruses.
Several studies on the localization and movement of begomoviruses in host plant system have concluded that the virus is limited to cambium and phloem parenchyma cells of plants. Sometimes, they exit the phloem and enter the mesophyll parenchymatous tissue (Rojas et al., 2005;Rojas and Gilbertson, 2008). One advantage of horizontal virus transmission is the limited transmission of viruses among nearby plants. Compared with horizontal transmission, viruses transmitted by infected seeds can play a crucial role in viral stability from parents to offspring and between seasons (Simmons et al., 2011).
It is crucial to identify viral transmission methods, as it helps to understand the epidemiology and outbreak of viruses  (Kim et al., 2015). The transmission of begomoviruses by infected seeds has been limited for many years, and viruses can only be transmitted through insects, sap, and inoculation mechanisms (Stanley et al., 2001). However, multiple begomoviruses, such as Mung bean yellow mosaic virus, TYLCV, ToLCNDV (Sangeetha et al., 2018;Kil et al., 2020), and Dolichos yellow mosaic virus, are seed-transmissible viruses (Kil et al., 2016). It was recently discovered that TYLCV could be transmitted through infected tomato seeds (Kil et al., 2016;Pérez-Padilla et al., 2020). Transovarial transmission of viruses from vector parents to offspring is important for their epidemiology. The majority of begomoviruses are limited to transmit only by insect vector but TYLCV have been found to show transovarial transmission. It was observed that TYLCV access into the reproductive organ of its vector was mostly determined by the developmental stage of the whitefly ovary, and that TYLCV transmission to offspring increased with whitefly adult age. The precise interaction between viral coat protein with whitefly vitellogenin was required for virus entry into the whitefly ovary (Wei et al., 2017).

Role of insect vectors
The global increase in the population of insect vectors and the rapid evolution of virus variants are also associated with the worldwide distribution of begomoviruses . The role of insect vectors, which is the key factor in the spread of begomoviruses, also contributes to their evolution (Seal et al., 2006;Materatski et al., 2021). The adaptive response of co-evolution among viruses and insects is supported by the continual and circulative dispersal of the virus in the insect vector and the direct connection between begomoviruses and insects (Nawaz-ul-Rehman and Fauquet, 2009;Gupta et al., 2021).
Bemisia tabaci, is the key player in the global spread of begomoviruses. Its increased adaptability resulted in the introduction of many begomoviruses in previously unreported areas (Islam and Wu, 2017). Squash leaf curl virus, for example, was discovered in Israel in 2002 and rapidly spread to the Middle East. Some other latest examples of the global distribution of begomoviruses, which followed the invasion path of B. tabaci, include ToLCNDV from India to southern Spain (Hagen et al., 2008) and cucurbit leaf curl virus from the southwestern United States to Florida (Juárez et al., 2014).

The connection between virus complexes and betasatellites
Another aspect involved in begomovirus outbreaks is the connection of complex viruses with betasatellites. In TYLCD, the first betasatellite, i.e., CLCuGB associated with either TYLCV-IL or TYLCV-Mild in tomato plants was reported in very recent Frontiers in Plant Science 08 frontiersin.org times in Israel (Gelbart et al., 2020), as the Mediterranean basin and the Middle East are suggested to be hubs of both TYLCV complex root and diversification. This is a major issue for farmers globally, particularly native to the Mediterranean area. Currently, 61 betasatellite species have been reported. More than 90% of betasatellite species are found in China and the Indian subcontinent. The coordinated efforts to prevent the further spread of betasatellites into the complex's genetic pool can restrict the emergence of begomovirus-betasatellite disease complexes.

Recent introduction of Cotton leaf curl Gezira virus to the United States
Two complexes of CLCuD have been identified in the OW: the African and Asian complexes. Although CLCuD was originally identified in Africa in 1912 (Farquharson, 1912), the causal agent was not discovered until much later (Sattar et al., 2013). Only CLCuGV has been found in African cotton; however, few cotton samples have been examined, and the exact diversity of begomoviruses may be considerably greater than we currently recognize (Idris and Brown, 2002). CLCuGV is a geographically widespread virus from central Africa to Jordan that infects various plant species, including cotton, Abelmoschus esculenta (okra), Alcea (hollyhock), and Salix alba (white willow; Idris and Brown, 2002;Tahir et al., 2011).
CLCuGV-satellite complex members have lately been found in the African Sahel, Arabian Peninsula and the Middle East, Pakistan, and recently in okra in southern Texas, USA (Tiendrébéogo et al., 2010;Tahir et al., 2011;Idris et al., 2014;Villegas et al., 2019). The origin of this OW CLCuGV-satellite complex in NW okra plants is unknown; however, it was possibly introduced by diseased plant materials or virulent whiteflies brought into Southern Texas, which is known for its extensive cotton and vegetable production. The spread of this disease poses a significant threat because cotton and vegetable cultivars grown in the US have no resistance to CLCuGV (Villegas et al., 2019).

Geographical distribution of some important begomoviruses
The growing phenomenon of begomoviruses and the attack of native viruses on crops, such as tomato, cotton, or cassava, indicates that existing outbreaks are the outcome of current begomoviruses colliding with recently launched hosts or vectors instead of dynamic changes in geminiviruses (Seal et al., 2006;Nawaz-ul-Rehman and Fauquet, 2009). Begomoviruses exploit the same whitefly vector, are extremely recombinogenic, and are widespread in tropical and subtropical areas, posing a severe threat to global food security (Prasanna et al., 2010). Table 4 lists some examples of begomoviruses re-emerged in Asia in new hosts. Figure 5 shows the geographical distribution of some important begomoviruses in Asia.

Asia: a geographically diverse region for emerging begomoviruses
Several begomoviruses have emerged in Asian countries, including TYLCKaV in Thailand, Tomato yellow leaf curl Indonesia virus, Pepper yellow leaf curl Indonesia virus (PepYLCIV), and Papaya leaf curl virus (PaLCuV). These begomoviruses cause substantial economic losses in Asian countries and increase their host range through mutation and recombination. Brief descriptions of some emerging begomoviruses that are not widely known but very destructive are given below. Figures 6, 7 show the symptoms and phylogenetic analysis, respectively, of these emerging viruses.

Tomato yellow leaf curl Kanchanaburi virus
TYLCKaV is a bipartite begomovirus with a large host range identified in Thailand's pepper plants. Resistance and pathogenic determinants in host plants of TYLCKaV have been studied less extensively than those of TYLCV . TYLCKaV predominantly infects Capsicum species, Solanum melongena, and Solanum lycopersicum. The spread of TYLCKaV can be reduced using simple procedures, such as isolating greenhouses for transplant production from outdoor sources of begomoviruses. The destruction of whitefly  Geographical location: The map shows the distribution of some important begomoviruses in Asia.

FIGURE 6
Symptoms of begomoviruses: populations is another key management method that might be beneficial for controlling the spread of TYLCKaV (Díaz-Pendón et al., 2010). Figure 6A depicts the symptoms exhibited by TYLCKaV infection.

Pepper yellow leaf curl Indonesia virus
PepYLCIV is a bipartite begomovirus in which DNA-A and DNA-B play significant roles in infected plants (Fondong, 2013). PepYLCIV was first detected in Central Java in 2003 and was found to infect Capsicum annuum (chili peppers) on Sumatra Island in 2005 (De Barro et al., 2008). Chili peppers are largely grown in Indonesia to meet the strong market demand. This crop has recently been decimated in Indonesia due to PepYLCIV infection (Fondong, 2013). This pathogen can infect hot peppers, sweet peppers, tomatoes, and weeds, causing symptoms such as leaf curling, yellowing in young leaves, and stunting . Figure 6B depicts the symptoms of PepYLCIV in chili fields and plants.

Papaya leaf curl virus
PaLCuV is a monopartite begomovirus and is transmissible through B. tabaci on a wide range of hosts. It can infect Carica papaya L (papaya), Cestrum nocturnum, Capsicum species, Cyamopsis tetragonoloba, Nicotiana tabacum, and S. lycopersicum . It causes one of the most severe viral diseases affecting papaya crops in tropical areas. Infected papaya plants exhibit symptoms, such as yellowing, downward, or upward curling, and vein thickening, whereas severely infected plants have curled petioles and wilting (Varun et al., 2017). Symptoms can be seen in Figure 6C. Currently, the most common treatments used to Phylogeny of DNA-A: The DNA-A sequences of different endemic begomoviruses were downloaded and assembled by Muscle Alignment module in MEGA-11. The phylogeny was constructed by maximum likelihood method with 500 bootstrap replications. The resulting tree was exported and visualized in iTOL (Letunic and Bork, 2021

Tomato leaf curl Joydebpur virus
ToLCJoV is a monopartite virus that poses a significant threat to tomato production in eastern India. ToLCJoV was first identified in Bangladesh by Maruthi et al. (2005). It infects chili and Hibiscus cannabinus (kenaf) worldwide (Paul et al., 2009). The virus is transmitted only by whiteflies (Paul et al., 2009 Figure 6D (Hamim et al., 2020).

Tomato yellow leaf curl Thailand virus
TYLCTHV (Genus: Begomovirus, Family: Geminiviridae) has a monopartite genome that contains a single molecule of circular ssDNA (DNA-A). Bemisia tabaci is responsible for its transmission. The major symptoms of TYLCTHV infection are leaf curling and yellowing (Khan et al., 2012). The virus was first reported in Thailand in 1994 (Blawid et al., 2008) and has now dispersed to Myanmar (Green et al., 2001), South China (Guo et al., 2009), and was recently identified in Taiwan in 2005 (Jan et al., 2007). As a prominent tomato-infecting begomovirus, TYLCTHV alters the dynamics of pre-tomato begomoviruses and adversely affects Taiwan's tomato production. TYLCTHV is probably more pathogenic and invasive, with a broad host range, and it outlasts and substitutes Tomato yellow leaf curl Taiwan virus, possibly due to more efficient transmission by indigenous or introduced whitefly biotypes (Tsai et al., 2011). Figure 6E shows symptoms produced by TYLCTHV in tomato.

Discussion
Begomovirus disease complexes are rapidly growing in terms of geographical distribution and host range. For example, TYLCV is spreading to new areas in the Indian and Pacific regions, including New Caledonia, Australia, and Mauritius . Before 2010, ToLCNDV was only limited to Asia but now it has been spread throughout Europe and has increased its host plant range, including Euphorbiaceae, Cucurbitaceae, Fabaceae, and Malvaceae . Similarly, CMD has been widespread in Africa until 2015, when SLCMV was spread in Southeast Asia .
The presence of multiple hosts and expansion in the whitefly population around the globe has contributed to the spread of the begomoviruses. With the existence of such a wide population of begomoviruses along with their ability to exchange genetic material through recombination, new begomoviruses have also been emerging has the potential to spread as a pandemic. In addition to the evolution, human beings also played a significant role in the spread of begomovirus diseases in different ways, such as moving crops from their original location to another environment or through germplasm transfer to a different environment (Seal et al., 2006). This review highlighted the introduction of emerging begomoviruses in Asia that have potential to spread as pandemic. We also highlighted currently less important endemic viruses: TYLCKaV, PepYLCIV, PaLCuV, ToLCJoV, and TYLCTHV. The factors driving the emergence of plant viruses include global movement of plant materials and increase in the population of insect vectors. Another aspect involved in begomovirus outbreaks is the connection of complex viruses with betasatellites. Given the alarming rate at which novel begomoviruses are appearing in tomato, peppers, papaya, eggplant, chili and other crops, there is always an urgent need to create economic and ecologically appropriate strategies for successful virus control.

Future perspectives
Recent developments highlight the diversity in host ranges of begomoviruses with strong evidence indicating that an increasing number of host species reported with time in different regions in the world. It is envisioned to properly investigate the Intra-and inter-regional diversification of begomoviruses with applying preventive measures to control their transmission. Moreover, certain aspects i.e., associated satellite molecules with unique properties particularly those related to the specific morphological changes caused by infection, begomovirus interaction with insect vectors and transmission methodology adopted by begomoviruses would give an additional approach to investigating begomoviruses. It is the collective responsibility of the scientific community to develop a thorough plan and policy to counteract this before a devastating effect on food security and the global economy is realized.

Author contributions
MQ, AL, E-JK, and SL outlined and conceptualized the review theme. MQ wrote the first draft of the manuscript. BN, TV, GS, PH, E-JK, SJ, K-YL, C-WT, HD, TH, T-TA, NW, JL, S-MK, MN-u-R, and SL contributed to the manuscript preparation and revision, and also read and approved the submitted version. All authors contributed to the article and approved the submitted version.

Conflict of interest
JL was employed by company NongWoo Bio. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


Complete chloroplast genome sequence of a bulbous flowering plant, Allium triquetrum Linnaeus, 1753 (Amaryllidaceae)

Abstract Allium triquetrum (Linnaeus, 1753) is a bulbous flowering plant of the genus Allium (Amaryllidaceae), native to the Mediterranean basin, and is now widespread and invasive in different parts of the world via ornamental horticultural trade. However, to date, the genomic study of A. triquetrum has lagged, which impedes the development of appropriate utilization and management practices for this species. Here, we report the complete chloroplast genome sequence of A. triquetrum. The chloroplast genome size of A. triquetrum was 153,298 bp, consisting of a pair of inverted repeat regions (26,547 bp), separated by a large single-copy (82,875 bp) region and a small single-copy (17,329 bp) region. Genome annotation predicted 133 genes, including 87 protein-coding genes, 38 tRNA genes, and eight rRNA genes. Phylogenetic analysis based on 60 whole chloroplast genome sequences of Allium species suggested that A. triquetrum and A. moly are sister to each other along with the clade of A. fasciculatum, A. hookeri, and A. macranthum.

Allium triquetrum Linnaeus, 1753 (Amaryllidaceae; also known as three-cornered leek) is a bulbiferous monocotyledonous herb native to the central and western Mediterranean regions of Eurasia (Tehranchian et al. 2014). Not only is this species valued as ornamental for its flowers, but it can also be traditionally used in folk medicine, especially for wound healing (Rabah et al. 2020). Due to anthropogenic propagule dispersal via the ornamental horticulture trade, followed by 'escape' into natural environments, this species has been invasive in different parts of the world (Lesieur et al. 2018(Lesieur et al. , 2022. However, to date, no genomic studies have reported on A. triquetrum, which impedes the development of its utilization and management practices (Lesieur et al. 2022). Here, we report the chloroplast genome sequence of A. triquetrum and reconstruct its phylogenetic relationship with other Allium species.
The chloroplast genome sequence of A. triquetrum was 153,298 bp in length and exhibited a typical quadripartite structure, consisting of a pair of inverted repeat regions (IRs) of 26,547 bp, separated by a large single-copy (LSC) region of 82,875 bp and a small single-copy (SSC) region of 17,329 bp. The chloroplast genome encoded 133 genes, of which 111 (77 protein-coding genes, 30 tRNA genes, and four rRNA genes) were unique and 22 (10 protein-coding genes, eight tRNA genes, and four rRNA genes) were duplicated in the IRs. The gene rps12 was trans-spliced with three exons: C-terminal exons 2 and 3 were located in the IRs, while exon 1 was $28.1 kilobase pairs downstream of the nearest copy of exons 2 and 3, and $69.8 kilobase pairs away from the distal copy of exons 2 and 3. The overall GC content was 36.8%, whereas the GC contents in IR, LSC, and SSC were 42.5%, 34.6%, and 29.9%, respectively.
The phylogenetic relationship of Allium was reconstructed using the maximum-likelihood (ML) method based on the multiple alignment of A. triquetrum and 59 previously reported chloroplast genomes of Allium, with Agapanthus coddii (KX790363) and Narcissus poeticus (MH706763) as outgroups. ML analysis was conducted using RAxML-HPC v.8.2.8 (Stamatakis 2014) with 1000 bootstrap replicates on the CIPRES Science Gateway (https://www.phylo.org/). The phylogenetic topology strongly supported that A. triquetrum and A. moly were sister to each other, and jointly sister to the clade of A. fasciculatum, A. hookeri, and A. macranthum (Figure 1).

Author contributions
ZQ and GS conceived the project. MX and LP collected samples and analyzed data. MX wrote the manuscript. ZQ and LP revised the manuscript. All authors read and approved the manuscript.

Disclosure statement
The authors are grateful for the open raw genome data from the public database. The authors reported no potential conflict of interest.

Data availability statement
The genome sequence data that support the findings of this study are openly available in GenBank of NCBI at https://www.ncbi.nlm.nih.gov under the accession no. ON310895. The associated BioProject, SRA, and Bio-Sample numbers are PRJNA825297, SRR18706756, and SAMN27511544, respectively.