,
Inter-species mRNA transfer among green peach aphids, dodder parasites, and cucumber host plants
Juan Songa,b,1, Jinge Biana,b,1, Na Xuea,b, Yuxing Xua,b, Jianqiang Wua,b     
a. Department of Economic Plants and Biotechnology, Yunnan Key Laboratory for Wild Plant Resources, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China;
b. CAS Center for Excellence in Biotic Interactions, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: mRNAs are transported within a plant through phloem. Aphids are phloem feeders and dodders (Cuscuta spp.) are parasites which establish phloem connections with host plants. When aphids feed on dodders, whether there is trafficking of mRNAs among aphids, dodders, and host plants and if aphid feeding affects the mRNA transfer between dodders and hosts are unclear. We constructed a green peach aphid (GPA, Myzus persicae)-dodder (Cuscuta australis)-cucumber (Cucumis sativus) tritrophic system by infesting GPAs on C. australis, which parasitized cucumber hosts. We found that GPA feeding activated defense-related phytohormonal and transcriptomic responses in both C. australis and cucumbers and large numbers of mRNAs were found to be transferred between C. australis and cucumbers and between C. australis and GPAs; importantly, GPA feeding on C. australis greatly altered inter-species mobile mRNA profiles. Furthermore, three cucumber mRNAs and three GPA mRNAs could be respectively detected in GPAs and cucumbers. Moreover, our statistical analysis indicated that mRNAs with high abundances and long transcript lengths are likely to be mobile. This study reveals the existence of inter-species and even inter-kingdom mRNA movement among insects, parasitic plants, and parasite hosts, and suggests complex regulation of mRNA trafficking.
Keywords: Dodder    Cuscuta    Green peach aphid    Host plant    mRNA transfer    Inter-kingdom interaction    
1. Introduction

The plant vascular system not only provides mechanical support, but also interconnects all organs. Through vasculature, various compounds, including signaling molecules, are delivered throughout the entire organism, enabling vascular plants to grow in large sizes and colonize a wide range of habitats (Jung and Park, 2007). In addition to water and nutrients, recent studies have revealed that a large number of macromolecules are also transported by plant vascular system, such as small RNAs, proteins, and mRNAs (Chitwood and Timmermans, 2010; Kehr and Kragler, 2018; Lin et al., 2009).

Increasing lines of evidence have indicated extensive mRNA trafficking through plant vascular system. Using library screening, a number of mRNAs were detected in the phloem sap of pumpkin (Cucurbita maxima), among which CmNACP, encoding a protein of NAC domain family, may be involved apical meristem development (Ruiz-Medrano et al., 1999). Next Generation Sequencing of Arabidopsis heterografts composed of shoots and roots with different genotypes further revealed 2128 mobile transcripts (Thieme et al., 2015). Although the biological significance of mRNA trafficking remains largely unclear, several mobile mRNAs have been identified to be functional in distal organs: Mobile mRNAs can regulate potato tuberization (Banerjee et al., 2006; Ghate et al., 2017; Hannapel and Banerjee, 2017), tomato leaf morphology and fruit development (Haywood et al., 2005; Kim et al., 2001), and Arabidopsis floral initiation (Huang et al., 2012) and lateral root formation (Yoshida et al., 2016).

Parasitic plants, which make up about 1% of angiosperm plants, attach to their host through haustoria which provide vascular connections between parasites and hosts, allowing exchange of substances such as water, nutrients, and macromolecules (proteins, mRNAs, viruses, and metabolites) between parasites and hosts (Clarke et al., 2019; Irving and Cameron, 2009; Westwood et al., 2010; Yoshida et al., 2016). Dodders (Cuscuta spp., Convolvulaceae) are shoot holoparasites with no leaves or roots, and nearly 200 dodder species are widely distributed worldwide (Clarke et al., 2019). Like many parasitic plants, dodders also use haustoria to form xylem and phloem connections with the host's vasculature (Birschwilks et al., 2006; Yoshida et al., 2016). Large-scale inter-plant mRNA translocation between dodder and its hosts has been reported (Liu et al., 2020; Thieme et al., 2015; Westwood and Kim, 2017). In an Arabidopsis-dodder (Cuscuta pentagona revised as C. pentagona) system, about 1% of the mRNAs in dodder stems near the haustoria are from Arabidopsis, while 0.6% of the mRNAs in the Arabidopsis stems are of C. pentagona origin (Kim et al., 2014). Similarly, over 2000 mRNAs of Arabidopsis were transferred to C. reflexa (Thieme et al., 2015).

Aphids (Hemiptera: Aphididae) are phloem-feeding pests, causing growth retardation or even death of host plants (Hales et al., 1997; Smith and Boyko, 2007). Aphids are also important virus vectors and their honeydew often leads to diseases in plants. Myzus persicae (green peach aphid, GPA) is a generalist phloem-feeding herbivore and is able to feed on more than 400 plant species, including various crops (Blackman and Eastop, 2000; Tang et al., 2015). Recently, inter-kingdom mRNA transfer from GPAs to host plants was detected. When GPAs fed on Arabidopsis plants, 3186 GPA-derived transcripts were detected in the feeding sites of Arabidopsis leaves (2985 mRNAs and 201 candidate lncRNAs); importantly, among these mobile RNAs from GPAs, a Ya long noncoding RNA could move to distal leaves, where it functions to increase GPA fecundity and performance on the whole plant (Chen et al., 2020).

GPAs are known to feed on dodders. GPA feeding on Cuscuta australis increased the resistance of soybean hosts to the cotton leafworm (Spodoptera litura) and soybean aphid (Aphis glycines), indicating translocation of certain GPA herbivory-induced systemic signals from C. australis to the host soybean, where defenses were activated by the systemic signals (Zhuang et al., 2018). Here, we show that in a GPA-C. australis-cucumber tritrophic interaction system, in which GPA fed on C. australis and cucumber was the host for C. australis, GPA feeding not only activated defense-related responses in C. australis and cucumbers, but large numbers of mRNAs were found to be transferred between C. australis and cucumbers and between C. australis and GPAs; importantly, we also show that a few mRNAs from GPAs and cucumbers were also respectively detected in cucumbers and GPAs as well, indicating extensive cross-species mRNA trafficking in this tritrophic system.

2. Materials and methods 2.1. Plant growth and aphid rearing

Dodder (Cuscuta australis) seed germination followed a previously published method (Li et al., 2015). Plants were grown in a glasshouse at ~17 ℃ (night) to ~28 ℃ (day) under a 14 h photoperiod. C. australis seedlings were infested on wild tomato (Solanum pennellii) plants to form stocks. Cucumber (Cucumis sativus cv. Chinese Long) seeds were germinated on moist filter paper, and 4-day-old seedlings were transplanted to perlite saturated with the Hoagland solution (HS) (Hoagland and Arnon, 1950) (Table S1 for recipe) for two weeks. The cucumbers were infested with excised C. australis stems from the stocks. Approximately a week later, dodder-infested cucumber plants were transferred to 1-L plastic pots filled with HS, which was renewed every 7 days until the third cucumber leaves had fully expanded (V3, 3-week-old seedlings).

The GPA (Myzus persicae) colony was originally established from a single first instar nymph collected from a field in the Kunming Institute of Botany, Chinese Academy of Sciences, in Kunming, China, in 2015. In order to minimize the risk that GPAs may carry viruses, the insect colony was reared successively on tobacco (Nicotiana tabacum), cabbage (Brassica oleracea), Chinese cabbage (Brassica rapa ssp. pekinensis), and sweet pepper (Capsicum annuum). Arabidopsis (Arabidopsis thaliana, Col-0) plants were provided as the food for GPAs for the month before the experiment.

2.2. Treatments and sample collection

GPA treatment on Cuscuta australis was performed as the following: For the treatment group, 50 adult GPAs were gently placed onto each C. australis and enclosed in a clip cage (3 cm diameter), at a position ca. 15 cm away from the host cucumber (Fig. 1A), and the plants of the control group (mock) were fixed with empty clip cages (Fig. 1A). Both control and treatment group were composed of 24 cucumber-C. australis parasitization pairs. After 24 h of GPA feeding, aphids and plant tissues (cucumber: roots and the second leaves; C. australis: the stems between the clip cage and host) were collected and each specific tissue from three individual plants were immediately pooled to form one replicate and all the samples were immediately frozen in liquid nitrogen. Samples were stored at −80 ℃ until being processed for phytohormone determination and RNA-seq analysis (five and three replicates respectively).

Fig. 1 JA, JA-Ile, and SA levels in Cuscuta australis and cucumber host plants in response to GPA feeding on C. australis (A) A schematic of the experimental setup. Empty clip cages or clip cages each containing 50 GPAs were attached to C. australis stems (15 cm from cucumber), forming the control and treatment group, respectively. After 24 h of feeding, cucumber (the second leaves and roots) and C. australis samples between the clip cages and cucumbers were harvested. The JA, JA-Ile, and SA contents in C. australis, cucumber leaves, and cucumber roots are respectively shown in (B) (C), and (D). Asterisks indicate significant differences between control and treatment group determined by Student's t-test (n = 5; ∗, P < 0.05; ∗∗, P < 0.01). Error bars are ± SE.
2.3. Phytohormone extraction and quantification

Samples were ground in liquid nitrogen, and approximately 150 mg of each sample was extracted with 1 mL of ice-cold ethyl acetate spiked with the internal standards (20 ng of D6 -JA, 5 ng of 13C6 - JA-Ile, and 5 ng of D4 -SA). After vortexing for 10 min, the samples were centrifuged at 13 000 g for 15 min at 4 ℃. The supernatants were transferred to fresh 2-mL microfuge tubes and evaporated to dryness in a vacuum concentrator (Eppendorf) at 30 ℃. The pellets were resuspended in 0.2 mL of 50% (v/v) methanol and vortexed for 10 min. The samples were centrifuged to remove particles at 13 000 g for 15 min at 4 ℃, and then immediately 0.1 mL of the supernatants were transferred to inserts in glass vials. The samples were analyzed on an UPLC-MS/MS system (LCMS-8040 system, Shimadzu), which was equipped with a Shim-pack XR-ODS III column (2.0 mm I.D. × 75 mm L., 1.6 μm, Shim-pack). The mobile phase composed of solvent A (0.05% formic acid and 5 mM ammonium formate in water) and solvent B (methanol) was used in a gradient mode for the separation. A negative electrospray ionization mode was used for detection. An ion with a specific mass-to-charge ratio generated from each endogenous phytohormone or internal standard (the parent ion) was selected and fragmented to obtain its daughter ions; a specific daughter ion was used for generating the corresponding compound’s chromatogram. The peak areas of internal standards and those of the respective target compounds were obtained, and concentration of each phytohormone was quantified by comparing its peak area with the peak area of its respective internal standard. Five biological replicates were used.

2.4. RNA-seq data acquisition and analysis

Total RNA was extracted using TRIzol reagent (ThermoFisher Scientific). For RNA-seq data acquisition, Illumina TruSeq RNA Sample Prep Kit (Illumina) was used to construct cDNA libraries. The generated cDNA libraries were sequenced on a NovaSeq 6000 platform (Illumina) to acquire sequence reads (6 G depth). The raw transcriptome data has been deposited at the Beijing Institute of Genomics under the accession number PRJCA003305. Trimmomatic (v.0.39) (Bolger et al., 2014) with default parameters was used to remove the adaptor sequences from the raw reads and low quality raw reads which contained more than 10% unknown bases (N). The Bioconductor DEseq2 (v.3.9) package (Love et al., 2014) was employed to infer differential gene expression. Transcripts with adjusted p values < 0.05 and absolute values of log2 (fold change) > 1 were selected as differentially expressed genes (DEGs).

In order to identify mobile mRNAs, all clean reads from a specific cDNA library were aligned to the source organism's genome using Hisat2 (v2.1.0) (Kim, Landmead and Salzberg, 2015), and the mapped reads were considered to be native transcripts and thus excluded from further analysis. The unmapped reads from each sample were mapped to both of the other two organisms' reference genomes; the reads that uniquely mapped to one of the two genomes were considered to be mobile mRNAs from that specific organism, but if certain reads mapped to both of the genomes, they were discarded from the analysis. The FeatureCounts from the Subread package (v.1.6) (Liao et al., 2014) was used to count reads for all genes in the reference genome, and only the mRNAs that appeared in at least two out of three replicates for each group of samples and with an average count number≥4 were chosen for further analyses (Kim et al., 2014). The M. persicae G006 genome v.3.0 (https://bipaa.genouest.org/sp/myzus_persicae_g006/), C. sativus genome (Li et al., 2019), and C. australis genome (Sun et al., 2018) were used in this study.

To estimate the proportions of mobile mRNAs in the three species' transcriptomes, the percentage of foreign mRNAs was calculated based on the ratio of the foreign reads to the sum of foreign and native reads according to the previously described methods (Kim et al., 2014; Liu et al., 2020). Furthermore, we determined the characteristics of the mobile mRNAs as the following: The exon lengths of all genes in an organism's genome were calculated according to the GFF annotation files (GPA, C. australis, and cucumber). The counts of raw reads were normalized to FPKM (fragments per kilobase of exon per million reads, FPKM) values using exon length according to the previously described methods (Maza et al., 2013; Van Verk et al., 2013).

Gene annotations of the mobile mRNAs were retrieved from M. persicae, C. australis, and C. sativus genome databases. Gene ontology enrichment was performed using the AgriGO (v.2.0) (http://systemsbiology.cau.edu.cn/agriGOv2/). Pathway enrichment analysis for cucumber and C. australis was respectively done using the Cucurbit Genomics Databases website (http://cucurbitgenomics.org/) and Plant MetGenMAP (http://bioinfo.bti.cornell.edu/cgi-bin/MetGenMAP/help.cgi). The bubble charts of enrichment results and Venn diagrams were plotted using web-based tools (http://www.bioinformatics.com.cn and http://bioinfogp.cnb.csic.es/tools/venny/).

3. Results 3.1. Aphid infestation-induced phytohormonal and transcriptional responses in Cuscuta australis and cucumber host

First, we determined whether GPA herbivory on dodder leads to communications between dodder and host plants. Since we were also interested in the responses of different organs (leaf and root) of a host plant, cucumber was chosen as the host, as cucumber is very amenable for hydroponic culture and root can be easily harvested. In the dodder (C. australis)-cucumber parasitization system, GPAs were infested on C. australis (treatment group) or C. australis was mock treated (control group) (illustrated in Fig. 1A; photo see Fig. S1). After 24 h of feeding, GPAs, C. australis, and cucumber hosts’ leaves and roots were harvested. Given that this study mainly aimed to study signaling and mRNA movement among GPA, C. australis, and cucumber, systemic C. australis stems, namely the segments between host stems and clip cages were harvested to avoid possible contaminations from GPAs and cucumber hosts.

Phytohormones play important roles in plant response to insect feeding (Erb et al., 2012; Wu and Baldwin, 2010). Thus, jasmonic acid (JA), JA-isoleucine conjugate (JA-Ile), and salicylic acid (SA) were quantified in the samples from the control and treatment group. Compared with those in the control group, in the treatment group, JA and JA-Ile contents in the stems of C. australis between clip cages and cucumber hosts increased about 1- and 1.6-fold, respectively, while the SA content decreased 48.7% (Fig. 1B). Although cucumber plants were not directly attacked by GPAs, JA, JA-Ile, and SA levels of cucumber leaves showed a similar response as did the C. australis stems: in the treatment group, JA and JA-Ile levels in cucumber leaves elevated about 66.4% and 33.2%, while the SA levels in cucumber leaves were 77.6% lower than in the control group (Fig. 1C); by contrast, JA and JA-Ile levels in cucumber roots decreased about 61.2% and 61.3% and SA contents slightly increased (13.1%) (Fig. 1D).

Furthermore, we determined whether GPA feeding on C. australis could activate transcriptional responses in both C. australis and cucumber hosts. Compared with C. australis in the control group, only 54 (25 up- and 29 downregulated) DEGs in C. australis were detected, and 11 of them were annotated (Table S2). Notably, 258 (222 up- and 33 downregulated) and 289 (90 up- and 199 downregulated) DEGs were found in the cucumber leaves and roots, respectively. Only 9 DEGs were commonly shared by cucumber leaves and roots, and 5 out of the 9 DEGs showed the opposite patterns of regulation. Further biological process and molecular function-related gene ontology (GO) term enrichment on the 249 unique DEGs in cucumber leaves indicated that “oxidoreductase activity, acting on NAD(P)H″, “NADH dehydrogenase activity”, “NADH dehydrogenase (ubiquinone) activity”, “oxidoreductase activity, acting on NAD(P)H, quinone or similar compound as acceptor”, and “chlorophyll binding” were enriched, and pathway enrichment analysis on these 249 unique DEGs indicated that eight pathways were significantly altered in cucumber leaves after GPA feeding on C. australis: “aerobic respiration I (cytochrome C)”, “photosynthesis light reactions”, “phenylpropanoid biosynthesis, initial reactions”, “scopolin and esculin biosynthesis”, “oxygenic photosynthesis”, “NAD/NADH phosphorylation and dephosphorylation”, “adenosine ribonucleotides de novo biosynthesis”, and “adenosine nucleotides de novo biosynthesis” (Fig. 2, Table S2). The 280 unique DEGs in cucumber roots were enriched in GO terms of biological process and molecular function “endopeptidase inhibitor activity”, “peptidase inhibitor activity”, “endopeptidase regulator activity”, “cysteine-type endopeptidase inhibitor activity”, “peptidase regulator activity”, “molecular function regulator”, and “enzyme inhibitor activity”, and the significantly changed pathways were “coumarin biosynthesis”, “taxiphyllin bioactivation”, “linustatin bioactivation”, and “neolinustatin bioactivation” (Fig. 2, Table S2).

Fig. 2 Comparison between the DEGs identified in cucumber leaves and roots. Experimental setup is illustrated in Fig. 1A. DEGs in cucumber leaves and roots were obtained by comparing their transcriptomes under control and treatment conditions. Venn diagram indicates the specifically and commonly regulated DEGs in cucumber leaves and roots. The tables below are the respective GO terms (molecular function) and pathways enriched from the specifically regulated DEGs. Complete data can be found in Table S2.

Thus, after GPA feeding, C. australis and cucumber hosts responded with phytohormonal and transcriptional changes. Moreover, we inferred that certain systemic signals were transmitted from the GPA-infested region of C. australis stem to the systemic tissues of C. australis and cucumbers, activating physiological changes.

3.2. Identification of cross-species mobile mRNAs in GPA-Cuscuta australis-cucumber system

In order to determine whether there is cross-species mRNA movement in the GPA-C. australis-cucumber tritrophic system, RNA-seq was also performed for the GPAs, in addition to the above-mentioned RNA-seq done for the systemic C. australis stems, cucumber leaves, and cucumber roots. Inter-species mobile mRNAs were identified by mapping the RNA-seq reads to the sample organism's genome to obtain non-mappable reads, which were further mapped to the other two organisms' genomes; if any of these reads could be mapped to one of the other two organisms but not both, they were considered to be translocated from that specific organism.

In the control group, 3672 and 257 mRNAs of C. australis were detected in cucumber leaves and roots, respectively, while only 15 cucumber mRNAs were found in C. australis (Fig. 3A, Table S3). After 24 h of GPA feeding, 3463 and 998 C. australis mRNAs, together with 3 and 2 GPA mRNAs, were respectively found in cucumber leaf and root samples, and 3818 cucumber mRNAs and 1181 GPA mRNAs were detected in the C. australis systemic stems which were collected between the clip cages and hosts; in GPAs, 3576 C. australis mRNAs and only 3 cucumber mRNAs were identified (Fig. 3A, Table S3). These data indicate that in this tritrophic interaction system, there was transport of large numbers of mRNAs between the two partners which directly interact (i.e., cucumber and C. australis, GPA and C. australis) and small numbers of mRNAs (less than five) can even reach the third partner (i.e., from cucumber to GPA and vice versa). However, the total abundances of mobile mRNAs in the recipient organisms were relatively low: In GPAs, C. australis, cucumber leaves and roots, foreign mobile mRNAs comprised only 0.22%, 1.30%, 0.30%, and 0.04% of the transcriptome, respectively (Fig. S2; Table S4).

Fig. 3 Identification of cross-species mobile mRNAs in the GPA-Cuscuta australis-cucumber system (A) Numbers of cross-species mobile mRNAs identified in the GPA-C. australis-cucumber system. The numbers in different colors indicate the numbers of mobile mRNAs from different species (red: C. australis-originated mobile mRNAs, green: cucumber-originated mobile mRNAs, blue: GPA-originated mobile mRNAs) (B) Venn diagram and GO analysis of the mobile mRNAs under control and treatment condition. Top panel: Cucumber mobile transcripts translocated into C. australis and C. australis mobile mRNAs identified in cucumber leaves and roots under control and treatment conditions. Bottom panel: Top 10 GO terms enriched from GPA feeding-induced unique mobile mRNAs translocated from cucumber to C. australis, from C. australis to cucumber leaves, and from C. australis to cucumber roots. Complete data can be found in Table S3.

To gain further insight into the regulation of plant mobile mRNAs by GPA feeding, Venn diagrams were used to analyze the mobile mRNAs identified in a specific group of plant samples, namely C. australis and cucumber leaves and roots, under control and GPA treatment. In C. australis, the 15 mRNAs from cucumber under control conditions were also mobile after GPA herbivory, and it is striking that GPA feeding induced many more (3803) unique mRNAs (Fig. 3B, Table S3), and these unique 3803 mobile cucumber mRNAs were found to be significantly enriched in the molecular function-related GO terms “catalytic activity”, “structural molecule activity”, “structural constituent of ribosome”, and many binding-related GO terms, such as “ion binding”, “nucleotide binding”, “nucleoside phosphate binding”, and “RNA binding” (Fig. 3B, Table S3). In cucumber leaves, 2299 dodder mRNAs were always transferred from C. australis, among which 1373 were no longer mobile after GPA herbivory, whereas 1164 mRNAs became mobile (Fig. 3B, Table S3). GO enrichment of the 1164 GPA-induced mobile mRNAs indicated that “catalytic activity” and a series of binding-related GO terms, including “protein binding”, “small molecule binding”, “nucleotide binding”, “nucleoside phosphate binding”, “carbohydrate derivative binding”, “ribonucleotide binding”, and “purine ribonucleotide binding” were enriched (Fig. 3B, Table S3). In cucumber roots (Fig. 3B, Table S3), 38 C. australis mRNAs’ inter-plant mobility was abolished by GPA feeding on C. australis, and 155 C. australis mRNAs were detected in both control and GPA treatment group; 843 C. australis mRNAs were exclusively found in cucumber roots after GPA infestation on C. australis, and molecular function-related GO term analysis of these 843 mRNAs showed that the terms “oxidoreductase activity”, “catalytic activity”, and binding-related GO terms similar to what were found in cucumbers leaves were enriched (Fig. 3B, Table S3). Therefore, GPA infestation on C. australis greatly altered the profile of mRNA translocation in C. australis and cucumber.

During feeding, aphids inject saliva into plant tissues, which contain various proteins and RNAs (Louis and Shah, 2013), and recent evidence indicated that aphid RNAs may function in plants to promote aphid fecundity (Chen et al., 2020). We found that 1181 GPA mRNAs were translocated to C. australis systemic stems (Fig. 4A, Table S5), and GO analysis on these insect-to-plant mobile mRNAs showed enrichment of molecular function terms “structural molecule activity”, “oxidoreductase activity”, “structural constituent of cuticle”, “polysaccharide binding”, “chitin binding”, and several transporter activity-related GO terms (Fig. S3; Table S5).

Fig. 4 Movement of mRNAs between GAPs and cucumbers through Cuscuta australis (A) Mobile mRNAs transported from GPAs to cucumbers through C. australis. Left: Venn diagram indicating the numbers of GPA mobile mRNAs in C. australis (very soft red) and in cucumber leaves (very soft blue) and roots (very soft organe); right: FPKM values of the three GPA mobile mRNAs (g16390 g2264, and g23359) detected in C. australis, cucumber leaves, and cucumber roots (B) Mobile mRNAs transported from cucumbers to GPAs through C. australis. Left: Venn diagram indicating the numbers of cucumber mobile mRNAs in C. australis (very soft red) and GPAs (very soft blue); right: FPKM values of the three cucumber mobile mRNAs (CsaV3_UNG167240, CsaV3_UNG240520, and CsaV3_UNG208380) that were detected in C. australis and GPAs. Different small letters indicate statistical differences (P < 0.05, one-way ANOVA followed by Duncan's test, n = 3). Complete data can be found in Table S5.

Phloem is the main channel that allows transfer of mRNAs (Haywood et al., 2005; Kehr and Buhtz, 2008; Notaguchi, 2015; Ruiz-Medrano et al., 1999; Westwood and Kim, 2017) and aphids are phloem feeders. Consistently, 2071 mobile mRNAs were commonly transferred from C. australis to both cucumber and to GPA (Fig. S4; Table S6). Notably, 1510 C. australis mRNAs were specifically detected in cucumber plants, and 1505 C. australis mRNAs were specifically translocated to GPAs (Fig. S4; Table S6), implying certain selectivity during mRNA transport between different interaction partners.

Our RNA-seq data revealed that even though cucumber plants and GPAs did not interact directly, C. australis indeed facilitated the transfer of mRNAs between them (Fig. 4, Table S5). Three mRNAs (g16390, g23359, and g22642) among the 1181 GPA mRNAs in C. australis were further translocated into cucumber leaves, and g16390 and g23359 transcripts were also found in cucumber roots (Fig. 4A, Table S5). g23359 is annotated as “ATP-dependent DNA helicase PIF1-like”, but the other two mRNAs are with no annotations (Table S7). The transcript level of GPA g16390 decreased as it being transported from C. australis to cucumber leaves and roots, but g22642 mRNA could only be detected in C. australis and cucumber leaves (Fig. 4A, Table S5). Similarly, only 3 cucumber mRNAs (CsaV3_UNG167240, CsaV3_UNG240520, and CsaV3_UNG208380) were detected in GPAs, and the FPKMs of CsaV3_UNG167240 and CsaV3_UNG240520 decreased as they moved from C. australis to GPAs, but the levels of CsaV3_UNG208380 in C. australis and GPAs were similar (Fig. 4B, Table S5). They are respectively annotated as retrotransposable element Tf2 (CsaV3_UNG167240), gamma-glutamyl hydrolase 2 (CsaV3_UNG240520), and photosystem II protein D1 (CsaV3_UNG208380) (Table S7).

These analyses suggest that mRNAs are transported in the GPA-C. australis-cucumber tritrophic system, and mRNAs can not only travel between plants, but also between insects and plants.

3.3. General properties of cross-species mobile mRNAs

It has been proposed that the transcripts with high abundances tend to be mobile (Calderwood et al., 2016). We examined whether this scenario holds true in our GPA-Cuscuta australis-cucumber system. The relative abundances were compared between mobile and nonmobile mRNAs and indeed we found that mobile mRNAs generally have greater abundances than do the nonmobile ones, and even for GPAs, mobile GPA mRNAs also exhibited greater levels than did the nonmobile ones (Fig. 5A, Table S8).

Fig. 5 Association between mRNA abundance and length with mobility. The FPKM values and transcript lengths were obtained from the following groups of mobile mRNAs: those from GPAs to Cuscuta australis, from C. australis to GPAs, from cucumbers to C. australis, from C. australis to cucumber leaves, and from C. australis to cucumber roots. The FPKM values and transcript lengths of nonmobile mRNAs for each organism were also obtained (A) Comparison of the average FPKM values between mobile and non-mobile mRNAs (B) Comparison of the average transcript length values between mobile and non-mobile mRNAs. Asterisks indicate statistically significant differences between the mobile and nonmobile mRNAs (∗∗P < 0.0001, ∗P < 0.05, two-tailed Student's t-test). Complete data can be found in Table S8.

In addition to mRNA abundance, another factor that may participate in the regulation of mobility is the transcript length (Calderwood et al., 2016). Statistical analysis indicated that in C. australis, cucumbers, and even GPAs, compared with those of nonmobile mRNAs, the genes of mobile mRNAs have longer exons, suggesting that long transcripts are more likely to move (Fig. 5B, Table S8).

4. Discussion

Taken advantage of the large phylogenetic distances between dodder and its wide range of host species, dodder-host systems can be used for studying systemic signaling (Hettenhausen et al., 2017; Li et al., 2020; Qin et al., 2019; Shen et al., 2020; Zhang et al., 2021) and identification of mobile molecules, including proteins and RNAs (Kim et al., 2014; Liu et al., 2020; Shahid et al., 2018). In this study, we show that many mRNAs trafficked between Cuscuta australis and cucumber hosts; importantly, GPA feeding on C. australis greatly altered the profile of mobile mRNAs between the parasites and hosts.

Previously, it was shown that in response to GPA feeding on the dodder C. australis, systemic defense responses in host soybean were induced, resulting in enhanced resistance of the soybean plants to S. litura and A. glycine (Zhuang et al., 2018). In this C. australis-soybean system, the JA contents in C. australis were decreased by GPA feeding (Zhuang et al., 2018). In contrast, in this study the JA contents in C. australis increased after GPA infestation. It is likely the discrepancy was due to different sampling positions in C. australis: In the previous study, the local C. australis stems, which were directly fed by GPAs, were collected (Zhuang et al., 2018), while in the current study, the systemic C. australis stems were harvested, which were between the clip cages and host stems.

mRNAs are reported to transport bidirectionally between hosts and Cuscuta (Kim and Westwood, 2015; Liu et al., 2020; Thieme et al., 2015; Westwood and Kim, 2017). Recent evidence from a Cuscuta campestris-cucumber system indicated that nitrogen stress strongly affected the number and species of mobile mRNAs between cucumber and C. campestris parasite (Zhang et al., 2021). Similarly, we found that GPA feeding also strongly regulates inter-species mRNA transfer: In the control group, there were only 15 cucumber mRNAs in C. australis, while in the treatment group, in which C. australis was infested with GPAs, 3818 cucumber mRNAs were identified in C. australis (Fig. 3A); the mRNAs translocated from C. australis to cucumber leaves and roots were also largely different between the control and treatment group (Fig. 2). It remains unclear why GPA feeding, nitrogen stress, and probably other environmental factors affect mRNA mobility. Compared with GPA herbivory-induced transcriptomic changes (only 54 DEGs in C. australis and 258 and 289 DEGs in the cucumber leaves and roots, respectively), changes in the numbers and species of mobile mRNAs are much greater. Thus, regulation of mobile mRNAs seems to be at least partly independent from the regulation of transcriptome.

The physiological and ecological significance of mRNA transfer between C. australis and cucumber and between C. australis and GPA is still unclear. Recently, it has been shown that more than 3186 transcripts from GPA were identified in the feeding site of the host plant Arabidopsis and among these mobile GPA transcripts, a long noncoding RNA Ya was found to be able to move into systemic leaves of host plant to promote aphid fecundity (Chen et al., 2020). During GPA feeding, in addition to various mRNAs translocated from GPAs to C. australis, we show more than 3500 C. australis mRNAs were found in GPAs (Fig. 3A). Aphids are phloem feeders. It is very likely that mRNAs and probably other types of RNAs in the phloem sap of plants (Chitwood and Timmermans, 2010; Kehr and Kragler, 2018; Westwood and Kim, 2017) enter GPA digestion track naturally together with the flow of phloem sap. Nevertheless, if any plant RNAs play a role of defense in aphids remains to be demonstrated. In our GPA-C. australis-cucumber system, three GPA mRNAs and three cucumber mRNAs also respectively migrated into cucumbers and GPAs. Given that most mobile GPA and cucumber mRNAs identified in C. australis did not further enter cucumbers and GPAs, respectively, it is worth studying why these three GPA mRNAs and three cucumber mRNAs are “super” mobile and whether they have any functions after entering foreign organisms.

Previous studies have suggested that mRNA abundance is likely a factor that influences mRNA mobility in Arabidopsis (Calderwood et al., 2016), and RNA mobility between Arabidopsis or tomato hosts and C. pentagona seems to well correlate with mRNA abundance too (Kim et al., 2014). In line with these findings, in our GPA-C. australis-cucumber system, for all organisms, mRNAs with relatively high abundances have greater likelihood of moving into the other species than those with relatively low abundances. mRNAs need to go through plasmodesmata to move across cells, and indeed, Calderwood et al. (2016) found that larger transcripts were less mobile (Calderwood et al., 2016). However, our analysis revealed that longer mRNAs exhibited greater tendency to move across species (Fig. 5B). One reason for this discrepancy is the difference of analyses: We directly used lengths of mobile and nonmobile transcripts for mobility comparison, while Calderwood et al. (2016) did not compare the mobile and nonmobile mRNAs, instead, only mobile mRNAs were used to analyze the correlation between the abundances and mRNA lengths, and they found that small but statistically significant negative correlation between mobile mRNA abundances and transcript lengths. Another possible reason is species dependency; namely, the shorter mRNAs in Arabidopsis are more mobile, while longer transcripts in GPA-C. australis-cucumber system are likely to be mobile. More mRNA mobility studies are needed to gain better insight into the mechanisms underlying mRNA trafficking.

Recent studies revealed extensive inter-plant mRNA exchanges between dodders and hosts (Kim and Westwood, 2015; Liu et al., 2020; Thieme et al., 2015; Westwood and Kim, 2017). Here, we show that in an inter-kingdom system (GPA-C. australis-cucumber), C. australis connections enable inter-species transfer of mRNAs among GPA, C. australis, and cucumber. miRNAs from dodder C. campestris move from C. campestris haustoria to the host stem and suppress host defense against C. campestris (Shahid et al., 2018). Recent evidence has also pointed to large-scale transfer of proteins between C. australis and hosts and between C. australis-connected different hosts (Liu et al., 2020). In this GPA-C. australis-cucumber tritrophic system, taken advantage of the large phylogenetic distances among these interaction partners, which allow identification of mobile molecules, whether proteins, small RNAs, and even secondary metabolites move and even may provide ecological/physiological significance in this interaction system is interesting to explore. Importantly, even though this study used a parasitism system, it is very likely that within an individual plant, aphid infestation or other stresses also activate changes in the mobile mRNAs qualitatively and quantitatively. Grafting systems could be used to examine the hypothesis.

Author contributions

J.W. and J.S. conceived the research. J. S., J.B., N.X., and Y.X performed the experiments. All authors analyzed the data. J.S., J.B., and J.W. wrote the manuscript.

Declaration of competing interest

The authors declare no conflict of interest.

Acknowledgements

This work was supported by China Postdoctoral Science Foundation (2020M673314 (J.S.) and 2020M673315 (Y.X.)) and Postdoctoral Directional Training Foundation of Yunnan Province (J.S. and Y.X.). We are grateful to the Service Center for Experimental Biotechnology at the Kunming Institute of Botany, CAS, for supporting plant cultivation.

Appendix A Supplementary data

The following is the Supplementary data to this article:

References
Banerjee, A.K., Chatterjee, M., Yu, Y.Y., et al., 2006. Dynamics of a mobile RNA of potato involved in a long-distance signaling pathway. Plant Cell, 18: 3443-3457. DOI:10.1105/tpc.106.042473
Birschwilks, M., Haupt, S., Hofius, D., et al., 2006. Transfer of phloem-mobile substances from the host plants to the holoparasite Cuscuta sp. J. Exp. Bot., 57: 911-921. DOI:10.1093/jxb/erj076
Blackman, R., & Eastop, V. (2000). Aphids on the World’s Crops
Bolger, A.M., Lohse, M., Usadel, B., 2014. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics, 30: 2114-2120. DOI:10.1093/bioinformatics/btu170
Calderwood, A., Kopriva, S., Morris, R.J., 2016. Transcript abundance explains mRNA mobility data in Arabidopsis thaliana. Plant Cell, 28: 610-615. DOI:10.1105/tpc.15.00956
Chen, Y., Singh, A., Kaithakottil, G.G., et al., 2020. An aphid RNA transcript migrates systemically within plants and is a virulence factor. Proc. Natl. Acad. Sci. U.S.A., 117: 12763-12771. DOI:10.1073/pnas.1918410117
Chitwood, D.H., Timmermans, M.C.P., 2010. Small RNAs are on the move. Nature, 467: 415-419. DOI:10.1038/nature09351
Clarke, C.R., Timko, M.P., Yoder, J.I., et al., 2019. Molecular dialog between parasitic plants and their hosts. Annu. Rev. Phytopathol., 57: 279-299. DOI:10.1146/annurev-phyto-082718-100043
Erb, M., Meldau, S., Howe, G.A., 2012. Role of phytohormones in insect-specific plant reactions. Trends Plant Sci., 17(5): 250-259. DOI:10.1016/j.tplants.2012.01.003
Ghate, T.H., Sharma, P., Kondhare, K.R., et al., 2017. The mobile RNAs, StBEL11 and StBEL29, suppress growth of tubers in potato. Plant Mol. Biol., 93(6): 563-578. DOI:10.1007/s11103-016-0582-4
Hales, D.F., Tomiuk, J., Wohrmann, K., et al., 1997. Evolutionary and genetic aspects of aphid biology: a review. Eur. J. Entomol., 94(1): 1-55.
Hannapel, D.J., Banerjee, A.K., 2017. Multiple mobile mRNA signals regulate tuber development in potato. Plants-Basel, 6(8): 1-17. DOI:10.3390/plants6010008
Haywood, V., Yu, T.S., Huang, N.C., et al., 2005. Phloem long-distance trafficking of Gibberellic acid-insensitive RNA regulates leaf development. Plant J., 42(1): 49-68. DOI:10.1111/j.1365-313X.2005.02351.x
Hettenhausen, C., Li, J., Zhuang, H.F., et al., 2017. Stem parasitic plant Cuscuta australis (dodder) transfers herbivory-induced signals among plants. Proc. Natl. Acad. Sci. U.S.A., 114(32): E6703-E6709. DOI:10.1073/pnas.1704536114
Hoagland, D.R., Arnon, D.I., 1950. The water-culture method for growing plants without soil. Calif. Agric. Exp. St., 347(5406): 1-32.
Huang, N.C., Jane, W.N., Chen, J., et al., 2012. Arabidopsis thaliana CENTRORADIALIS homologue (ATC) acts systemically to inhibit floral initiation in Arabidopsis. Plant J., 72(2): 175-184. DOI:10.1111/j.1365-313X.2012.05076.x
Irving, L.J., Cameron, D.D., 2009. You are what you eat: interactions between root parasitic plants and their hosts. Adv. Bot. Res., 50: 87-138. DOI:10.1016/S0065-2296(08)00803-3
Jung, J.H., Park, C.M., 2007. Vascular development in plants: specification of xylem and phloem tissues. J. Plant Biol., 50(3): 301-305. DOI:10.1007/Bf03030658
Kehr, J., Buhtz, A., 2008. Long distance transport and movement of RNA through the phloem. J. Exp. Bot., 59(1): 85-92. DOI:10.1093/jxb/erm176
Kehr, J., Kragler, F., 2018. Long distance RNA movement. New Phytol., 218(1): 29-40. DOI:10.1111/nph.15025
Kim, D., Landmead, B., Salzberg, S.L., 2015. HISAT: a fast spliced aligner with low memory requirements. Nat. Methods, 12(4): 357-360. DOI:10.1038/Nmeth.3317
Kim, G., LeBlanc, M.L., Wafula, E.K., et al., 2014. Genomic-scale exchange of mRNA between a parasitic plant and its hosts. Science, 345(6198): 808-811. DOI:10.1126/science.1253122
Kim, G., Westwood, J.H., 2015. Macromolecule exchange in Cuscuta-host plant interactions. Curr. Opin. Plant Biol., 26: 20-25. DOI:10.1016/j.pbi.2015.05.012
Kim, M., Canio, W., Kessler, S., et al., 2001. Developmental changes due to long-distance movement of a homeobox fusion transcript in tomato. Science, 293(5528): 287-289. DOI:10.1126/science.1059805
Li, J., Hettenhausen, C., Sun, G.L., et al., 2015. The parasitic plant Cuscuta australis is highly insensitive to abscisic acid-induced suppression of hypocotyl elongation and seed germination. PLoS One, 10(8): 1-12. DOI:10.1371/journal.pone.0135197
Li, Q., Li, H., Huang, W., et al., 2019. A chromosome-scale genome assembly of cucumber (Cucumis sativus L.). GigaScience, 8(6): 1-10. DOI:10.1093/gigascience/giz072
Li, S., Zhang, J., Liu, H., et al., 2020. Dodder-transmitted mobile signals prime host plants for enhanced salt tolerance. J. Exp. Bot., 71(3): 1171-1184. DOI:10.1093/jxb/erz481
Liao, Y., Smyth, G.K., Shi, W., 2014. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics, 30(7): 923-930. DOI:10.1093/bioinformatics/btt656
Lin, M.K., Lee, Y.J., Lough, T.J., et al., 2009. Analysis of the pumpkin phloem proteome provides insights into angiosperm sieve tube function. Mol. Cell. Proteomics, 8(2): 343-356. DOI:10.1074/mcp.M800420-MCP200
Liu, N., Shen, G., Xu, Y., et al., 2020. Extensive inter-plant protein transfer between Cuscuta parasites and their host plants. Mol. Plant, 13(4): 573-585. DOI:10.1016/j.molp.2019.12.002
Louis, J., Shah, J., 2013. Arabidopsis thaliana-Myzus persicae interaction: shaping the understanding of plant defense against phloem-feeding aphids. Front. Plant Sci., 4: 213. DOI:10.3389/fpls.2013.00213
Love, M.I., Huber, W., Anders, S., 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol., 15(12): 550. DOI:10.1186/s13059-014-0550-8
Maza, E., Frasse, P., Senin, P., et al., 2013. Comparison of normalization methods for differential gene expression analysis in RNA-Seq experiments: a matter of relative size of studied transcriptomes. Commun. Integr. Biol., 6(6): e25849. DOI:10.4161/cib.25849
Notaguchi, M., 2015. Identification of phloem-mobile mRNA. J. Plant Res., 128(1): 27-35. DOI:10.1007/s10265-014-0675-6
Qin, Y., Zhang, J., Hettenhausen, C., et al., 2019. The host jasmonic acid pathway regulates the transcriptomic changes of dodder and host plant under the scenario of caterpillar feeding on dodder. BMC Plant Biol., 19(1): 540. DOI:10.1186/s12870-019-2161-8
Ruiz-Medrano, R., Xoconostle-Cazares, B., Lucas, W.J., 1999. Phloem long-distance transport of CmNACP mRNA: implications for supracellular regulation in plants. Development, 126(20): 4405-4419. DOI:10.1007/s004290050293
Shahid, S., Kim, G., Johnson, N.R., et al., 2018. MicroRNAs from the parasitic plant Cuscuta campestris target host messenger RNAs. Nature, 553(7686): 82-85. DOI:10.1038/nature25027
Shen, G., Liu, N., Zhang, J., et al., 2020. Cuscuta australis (dodder) parasite eavesdrops on the host plants’ FT signals to flower. Proc. Natl. Acad. Sci. U.S.A., 117(37): 23125-23130. DOI:10.1073/pnas.2009445117
Smith, C.M., Boyko, E.V., 2007. The molecular bases of plant resistance and defense responses to aphid feeding: current status. Entomol. Exp. Appl., 122(1): 1-16. DOI:10.1111/j.1570-7458.2006.00503.x
Sun, G., Xu, Y., Liu, H., et al., 2018. Large-scale gene losses underlie the genome evolution of parasitic plant Cuscuta australis. Nat. Commun., 9(1): 2683. DOI:10.1038/s41467-018-04721-8
Tang, Q., Xiang, M., Hu, H., et al., 2015. Evaluation of sublethal effects of sulfoxaflor on the green peach aphid (Hemiptera: Aphididae) using life table parameters. J. Econ. Entomol., 108(6): 2720-2728. DOI:10.1093/jee/tov221
Thieme, C.J., Rojas-Triana, M., Stecyk, E., et al., 2015. Endogenous Arabidopsis messenger RNAs transported to distant tissues. Nat. Plants, 1(4): 15025. DOI:10.1038/Nplants.2015.25
Van Verk, M.C., Hickman, R., Pieterse, C.M., et al., 2013. RNA-Seq: revelation of the messengers. Trends Plant Sci., 18(4): 175-179. DOI:10.1016/j.tplants.2013.02.001
Westwood, J.H., Kim, G., 2017. RNA mobility in parasitic plant-host interactions. RNA Biol., 14(4): 450-455. DOI:10.1080/15476286.2017.1291482
Westwood, J.H., Yoder, J.I., Timko, M.P., et al., 2010. The evolution of parasitism in plants. Trends Plant Sci., 15(4): 227-235. DOI:10.1016/j.tplants.2010.01.004
Wu, J.Q., Baldwin, I.T., 2010. New insights into plant responses to the attack from insect herbivores. Annu. Rev. Genet., 44: 1-24. DOI:10.1146/annurev-genet-102209-163500
Yoshida, S., Cui, S., Ichihashi, Y., et al., 2016. The haustorium, a specialized invasive organ in parasitic plants. Annu. Rev. Plant Biol., 67: 643-667. DOI:10.1146/annurev-arplant-043015-111702
Zhang, J., Xu, Y., Xie, J., et al., 2021. Parasite Cuscuta campestris enables transfer of bidirectional systemic nitrogen signals between host plants. Plant Physiol., 185(4): 1395-1410. DOI:10.1093/plphys/kiaa004
Zhuang, H., Li, J., Song, J., et al., 2018. Aphid (Myzus persicae) feeding on the parasitic plant dodder (Cuscuta australis) activates defense responses in both the parasite and soybean host. New Phytol., 218(4): 1586-1596. DOI:10.1111/nph.15083