Large-scale MRNA Transfer Between Haloxylon Ammodendron (Chenopodiaceae) And Herbaceous Root Holoparasite Cistanche Deserticola (Orobanchaceae)

Highlights
•Ten thousand mRNAs transfer between Haloxylon ammodendron and root parasitic Cistanche deserticola

•RNA mobility and functional analysis were performed in sunflower-Orobanche cumana system

•CdNLR1 and CdNLR2 would cause root-specific HR and affect parasitic equilibrium

Cistanche Cultivation and Analysis

Summary
Exchanges of mRNA were shown between host and stem parasites but not root parasites. Cistanche deserticola (Orobanchaceae) is a holoparasitic herb which parasitizes on the roots of woody plant Haloxylon ammodendron (Chenopodiaceae). We used transcriptome sequencing and bioinformatic analyses to identify nearly ten thousand mobile mRNAs. Transcript abundance appears to be a driving force for transfer event and mRNA exchanges occur through haustorial junction. Mobility of selected mRNAs was confirmed in situ and in sunflower-Orobanche cumana heterologous parasitic system. Four C. deserticola →H. ammodendron mobile mRNAs appear to facilitate haustorium development. Of interest, two mobile mRNAs of putative resistance genes CdNLR1 and CdNLR2 cause root-specific hypersensitive response and retard parasite development, which might contribute to parasitic equilibrium. The present study provides evidence for the large-scale mRNA transfer event between a woody host and a root parasite, and demonstrates the functional relevance of six C. deserticola genes in host-parasite interactions.

Subject areas
Plant biologyInteraction of plants with organismsOmicsTranscriptomics
Introduction
Increasing evidence suggests the importance of macromolecules such as proteins and mRNAs in both short- and long-distance communication in plants.1 Plasmodesmata are considered as short-distance transport channels, whereas vascular system carries out long-distance molecular transport between distal tissues.2,3 mRNA was reported to transfer in phloem.4 Recent studies have shown that the stem-loop structure of tRNA origin endowed mRNA with long-distance transport capacity.5 AtRRP44A, a subunit of the RNA exosome, was reported to interact with plasmodesmata and to mediate the cell-to-cell trafficking of KNOTTED1 (KN1) mRNA.6

Parasitic plants account for about 1% of flowering plants. The most important feature of parasitic plants is a specialized structure called haustorium, which establishes physical and physiological connection channels between the host and parasitic plants and so dominates most of their interactions.7 Parasitic plants rely on the host to sustain their growth, absorbing water and the nutrients such as photosynthates, amino acids and other intermediate metabolites from the host through the haustorium.8 Haustoria also function as the connection bridge between the hosts to transmit signaling molecule such as herbivore signal.9 It was also shown that the biomolecules including proteins and RNAs exchange bidirectionally through host-parasite connections. The earliest example of RNA transfer between the host and parasitic plants is the transmission of RNA virus from the infected host to the uninfected plant through dodder.10 Small RNAs (sRNAs) such as small interference RNA (siRNA) and microRNA (miRNA) were shown to move between host-parasite to regulate the target gene expression in recipient organisms.11,12,13 As for mRNA transfer, mobile mRNAs from the hosts tomato, pumpkin and alfalfa to the parasite Cuscuta chinensis have been demonstrated for several genes in the initial test.14 Through microarray analysis, Roney et al. found that 474 mRNAs were transferred from tomato to dodder.14 However, large-scale identification of transfer mRNAs was achieved only on the uses of next-generation sequencing technology. Kim et al. used transcriptome sequencing technology to identify large-scale and bidirectional mRNA transfer between the parasite Cuscuta pentagona and the hosts such as Arabidopsis thaliana and tomato.15 But the functional relevance of transfer mRNAs in host-parasite interactions was not clear.

Orobanchaceae is the largest parasitic angiosperm family in which many are facultative or obligate root parasites.16 Cistanche is a worldwide genus of holoparasitic desert plants in Orobanchaceae. C. deserticola is a holoparasite which parasitizes on the roots of a woody psammophyte, Haloxylon ammodendron (Chenopodiaceae).17 It was shown that chloroplast rpoC2 gene was transferred from H. ammodendron to C. deserticola by horizontal gene transfer.18 Here, we identified nearly ten thousand mobile mRNAs between C. deserticola and H. ammodendron through next-generation transcriptome sequencing and bioinformatic analysis. The mobile mRNAs were ascertained by multiple sequence alignment, PCR validation, and especially utilization of a sunflower-Orobanche cumana parasitic system. The function of mobile mRNAs was also demonstrated for several genes by this heterologous parasitic system. To date, the reports on functional mRNAs exchange between host and parasite were very limited. Therefore, our study provides new insights into the parasitic mechanism from the view point of mobile mRNAs.

Results
Identification of mobile mRNAs between C. deserticola and H. ammodendron
Tight physical connection of host H. ammodendron root and parasite C. deserticola makes it impossible to separate them at the haustorium. Therefore, we combined high-throughput RNA-sequencing (RNA-seq) and stepwise bioinformatic classification to identify the mobile transcripts between host and parasite. Four different set of samples were collected for RNA-seq analyses on Illumina HiSeq2500 platform. These include the root tissue of parasite-free host H. ammodendron (HA), succulent stem of C. deserticola (CD), the root tissue of host H. ammodendron parasitized with C. deserticola (HC), and the haustorial interface (HI). For HC and CD samples, the respective host root and parasite stem were collected 1 cm away from the parasitic junction (Figure 1A).

Figure 1. RNA-seq analysis and confirmation of mobile RNAs

(A) Illustration of select tissues for RNA-seq analysis. CD, fresh succulent stems of C. deserticola near haustorial junction; HC, the roots of CD-infested H. ammodendron; HI, haustorial interface. Scale bars are indicated.

(B) Outline of the strategies for sequence assembly and confirmation of mobile RNAs. RNA-seq, next generation RNA sequencing. ISO-Seq, full-length transcriptome sequencing. HA, the roots of non-parasitized H. ammodendron. HA.FL_NR, full-length transcriptome sequence of H. ammodendron. CD.FL_NR, full-length transcriptome sequence of C. deserticola.

(C) ORF prediction and further confirmation of assembly result. ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/) was applied for ORF prediction, and above assembled unigenes with the expression threshold value of (FPKM≥0.3) were cross-checked. The overlapped unigenes were further filtered with dual BLASTN (E value = 1e−10) against Cis and HAC unigenes.

(D) Venn diagrams showing common and differed sets of transcripts between CD, HA and HC samples.

Reliable identification of source tissue for a given transcript from interacting organisms is a challenging issue. Therefore Ikeue et al. have developed a useful bioinformatic method to distinguish host and parasite transcripts.¹⁹ Here we used a modified bioinformatic approach which is based on their method to identify the mobile mRNAs between C. deserticola and H. ammodendron (Figure 1B). A grand total of 723,382,940 clean reads were generated from above libraries (Table S2), and these reads were assembled via different strategies (Figure 1B). As a result, all ten samples were hybrid assembled into 222,899 unigenes (termed “Combined” hereafter, Table S3). C. deserticola and H. ammodendron samples were also assembled into 107,752 and 194,720 unigenes respectively (termed “Cis” and “HAC” hereafter, Table S3). To visualize the variation as well as the similarity for all samples, we performed a principal component analysis (PCA) on the normalized FPKM values of all the detected genes in “Combined unigenes”. The PCA plot showed that the data for three biological replicates were clustered closely and were separated in different samples, especially between H. ammodendron and C. deserticola (Figure S1). Meanwhile, ORF prediction showed that 149,825 (67.23%) of the “Combined” unigenes were assumed to encode putative proteins (Figure 1C). The individually assembled unigenes with the expression threshold value of FPKM≥0.3 from four different samples were cross-checked with the ORF prediction results (Tables S4–S6). This analysis identified that 69.01–75.90% of the unigenes from the four samples CD, HC, HA and HI had putative ORF (Figure 1C). Full-length transcriptome sequencing was carried out and the data, CD_FL for C. deserticola and HA_FL for H. ammodendron, were used for assembly error correction for both host H. ammodendron and parasite C. deserticola as both species lack genomic data (Tables S8–S11). When the assigned unigenes from above four samples were filtered with dual BLASTN (E value = 1e⁻¹⁰) against Cis, HAC, CD_FL and HA_FL, 94.78–99.00% of them had reliable consensus sequences (Figure 1C and Table S7). Based on the above analysis, 17,379 (HA-exclusive 14,810 and HA-containing 2,569) common unigenes were finally retrieved from HC and CD on Venn diagram analysis of the unigenes from CD, HC and HA (Figure 1D). Among them, the former 14,810 unigenes are most probably of parasite CD origin because of their absence in HA, but a proportion of them might also represent HC→CD mobile mRNAs which are present only in HC but not in HA owing to their up-regulation on parasite attachment. The latter 2,569 unigenes most probably represent unidirectional HC→CD transfer mRNAs as they are also present in HA. The other unigenes were excluded from candidate mobile mRNAs as they were not shared between HC and CD. One should also note that 14,810 common unigenes in [HC, CD] outnumbered 166 common unigenes in [HA, CD] by 89 times, highly suggesting the reliability of our analyses as the latter was logically impossible and was caused by the inevitable errors from unigene assembly or bioinformatic analyses (Figure 1D).

The origin of putative mobile mRNAs was further confirmed via dual BLAST analyses to assign their sequence origin (Figure 1B and Table S12). Because of the lack of genomic data for both host and parasite, we utilized available sequence data for Chenopodiaceae and Orobanchaceae family organisms to which H. ammodendron and C. deserticola belong respectively. Our hypothesis is that mobile mRNAs of H. ammodendron origin most probably have higher homologies with their orthologs from other Chenopodiaceae species than those from Orobanchaceae species, while it is vice versa for mobile mRNA of C. deserticola origin, and the conserved orthologs probably have high similarity with both Orobanchaceae and Chenopodiaceae. This analysis could help us to ascertain the origins of the mobile mRNAs as they could be assigned to family level via homology searching. To prove our hypothesis, phylogenetic tree of thirty-seven species, including six Orobanchaceae and four Chenopodiaceae species, was constructed using OrthoFinder. The far evolutionary distance between Orobanchaceae and Chenopodiaceae species ascertained the reliability of our hypothesis (Figure 2A and Table S13). Gene loss analysis showed that the transcripts for many orthogroups and photosynthesis related genes were absent in C. deserticola (Figures 2B and 2C, Table S14, Data S2 and S3). Moreover, the Venn diagram analysis between C. deserticola and three other sequenced parasitic plants also showed that 53.84% (926/1,720) of the missing orthogroups in C. deserticola were also absent in other three parasitic plants, namely 50.81% (755/1486), 48.09% (377/784) and 50.00% (260/520) in C. australis, Striga asiatica and Phtheirospermum japonicum respectively (Figure S2). These results not only indicated the parasitic property of this species but also the reliability of unigene assembly from our analyses, and also supported our hypothesis to use relative plant species for confirmation of sequence origin. Therefore, the candidate transfer unigenes from above analyses were separated with dual BLASTN (E value = 1e⁻¹⁰) against a collection of public available sequence datasets from Orobanchaceae (transcriptome of Triphysaria versicolor, Striga hermonthica and Phelipanche aegyptiaca; EST and mRNA sequences from NCBI) and Chenopodiaceae (next-generation and full-length transcriptome of H. ammodendron; EST and mRNA sequences from NCBI) (Figure 1B and Table S12). As a result, 7,496 unigenes were confirmed to originate from the parasite C. deserticola, accounting for 9.66% (7,496/77,615) and 13.70% (7,496/54,721) of the unigenes in destination HC samples and source CD, respectively; and 2,370 unigenes were assigned to host H. ammodendron, accounting for 3.38% (2,370/70,119) and 4.15% (2,370/57,091) of the unigenes in source HC and destination CD samples, respectively (Figures 2D and 2E, Tables S7, S12, and S15). The common 2,931 unigenes were too homologous to be assigned to the exact source organism. These results indicated that nearly ten thousand (7,496 + 2,370 = 9,866) unigenes could transfer between root parasitic plant C. deserticola and woody plant host H. ammodendron. Much more (7,496/2,370 = 3.16-folds) mobile RNAs were of parasite C. deserticola origin than of host H. ammodendron origin, showing the mRNA transfer bias in a parasite→host direction.

Figure 2. Gene loss analysis and further confirmation of mobile RNAs via multiple alignment

(A) Phylogeny of flowering plant genomes. OrthoFinder was used to analyze the evolution of all protein sequences of 37 species, and to find the orthogonal group including direct homologous genes and collateral homologous genes among species.

(B) Percentage of lost orthogroups. The number of orthogroup deletions in semi parasitic plants and total parasitic plants.

(C) Percentage of lost unigenes. The proportion of genes related to photosynthesis and metabolism was shown.

(D) Venn diagrams showing common and differed sets of transcripts between CD and HC samples on multiple alignment. CD_trans, CD→HC mobile RNAs. HC_trans, HC→CD mobile RNAs.

(E) Pie charts show the proportions of mobile reads mapped to CD (outer circle) and HC (inner circle) transcriptomes.