Identification of miRNAs and Their Targets in the Liverwort Marchantia polymorpha by Integrating RNA-Seq and Degradome Analyses

Abstract

Bryophytes (liverworts, hornworts and mosses) comprise the three earliest diverging lineages of land plants (embryophytes). Marchantia polymorpha, a complex thalloid Marchantiopsida liverwort that has been developed into a model genetic system, occupies a key phylogenetic position. Therefore, M. polymorpha is useful in studies aiming to elucidate the evolution of gene regulation mechanisms in plants. In this study, we used computational, transcriptomic, small RNA and degradome analyses to characterize microRNA (miRNA)-mediated pathways of gene regulation in M. polymorpha. The data have been integrated into the open access ContigViews-miRNA platform for further reference. In addition to core components of the miRNA pathway, 129 unique miRNA sequences, 11 of which could be classified into seven miRNA families that are conserved in embryophytes (miR166a, miR390, miR529c, miR171-3p, miR408a, miR160 and miR319a), were identified. A combination of computational and degradome analyses allowed us to identify and experimentally validate 249 targets. In some cases, the target genes are orthologous to those of other embryophytes, but in other cases, the conserved miRNAs target either paralogs or members of different gene families. In addition, the newly discovered Mpo-miR11707.1 and Mpo-miR11707.2 are generated from a common precursor and target MpARGONAUTE1 (LW1759). Two other newly discovered miRNAs, Mpo-miR11687.1 and Mpo-miR11681.1, target the MADS-box transcription factors MpMADS1 and MpMADS2, respectively. Interestingly, one of the pentatricopeptide repeat (PPR) gene family members, MpPPR_66 (LW9825), the protein products of which are generally involved in various steps of RNA metabolism, has a long stem-loop transcript that can generate Mpo-miR11692.1 to autoregulate MpPPR_66 (LW9825) mRNA. This study provides a foundation for further investigations of the RNA-mediated silencing mechanism in M. polymorpha as well as of the evolution of this gene silencing pathway in embryophytes.

Keywords: ARGONAUTE; Class III homeodomain leucine zipper; Degradome; MADS-box; Marchantia polymorpha; Transcriptome; miRNA prediction.

Fig. 1

Characterization of the Marchantia polymorpha transcriptome and the annotation of open reading frames (ORFs). (A) Four-week-old M. polymorpha plants (Tak1 accessions). Scale bar = 2 cm. (B) Pie chart showing the percentages of full-length transcripts (12.5%), partial transcripts (11.7%) and undefined transcripts (75.8%) from the 39,090 de novo assembled M. polymorpha contigs. (C) Gene length distribution comparison between M. polymorpha and Arabidopsis thaliana. Only the coding DNA sequences (CDS) of full-length transcripts were calculated. (D) Amino acid similarity of the CDS from Arabidopsis and M. polymorpha. (E) Size distribution of small RNAs in M. polymorpha.

Fig. 2

RNA-Seq workflow for miRNA–target prediction in Marchantia polymorpha. (A) Four-week-old thalli of M. polymorpha were selected for whole-transcriptome, small RNA and degradome analyses using deep sequencing. The ContigViews platform was employed to annotate open reading frames (ORFs) from a set of 39,090 de novo assembled M. polymorpha contigs. RNA secondary structure was predicted with the RNAfold program from the Vienna RNA web server. MiRNA target prediction was performed with the psRNATarget web server. (B) The algorithm for miRNA prediction. (C) Pie chart illustrating the distribution of the first 5′-end nucleotides of M. polymorpha miRNA.

Fig. 3

Detection of miRNA 3′-end 2′-O-methylation in Marchantia polymorpha. β-Elimination assay of miR166 in Col-0, transgenic Arabidopsis expressing the P1/HC-Pro gene (P1/HC), and M. polymorpha. The upper panel shows an X-ray film after a long exposure (12 h). The lower panel shows an X-ray film after a short exposure (4 h).

Fig. 4

Four conserved miRNAs in Marchantia polymorpha. (A) MiRNA precursor structures. Sequences in bold black represent miRNAs, while sequences in bold gray represent the miRNA* strands. Read counts are indicated in parentheses. (B) The detection of conserved miRNAs by Northern blot in M. polymorpha. U6 was used as a loading control. The asterisk indicates that the deep sequence of the small RNA comprises >700 read counts. (C) Conserved miRNAs and their pairwise target regions. The dashed line indicates the 10th and 11th positions of the miRNA. The arrowhead and the number indicate the position observed in the degradome and the number of read counts, respectively. (D) Degradome map of the miRNA target. The black line indicates significant (binomial test, P-value < 0.001) degradome read counts at the 10th and 11th positions of the target site. Degradome reads were normalized to the total degradome read counts.

Fig. 5

Four novel miRNAs that regulate AGO1 and MADS genes in Marchantia polymorpha. (A) Four novel miRNA precursor structures. Sequences in bold black represent miRNAs, while sequences in bold gray represent miRNA* strands. Read counts are indicated in parentheses. (B) Four novel miRNAs and their pairwise target regions. The dashed line indicates the 10th and 11th positions of the miRNA. The arrowhead and the number indicate the position observed in the degradome and the number of read counts, respectively. (C) The detection of novel miRNAs by Northern blot in M. polymorpha. U6 was used as a loading control. The asterisk indicates that the miRNA sequence comprises >700 read counts. (D) Degradome map of the miRNA target. The black line indicates the significant (binomial test, P-value < 0.001) degradome read counts at the 10th and 11th positions of the target site.

Fig. 6

Reporter assay conducted to monitor the miRNA-mediated down-regulation of predicted targets by transient expression in Nicotiana benthamiana. (A) Schematic of miRNAs pairing with their target sites on the YFP fusion constructs. The red nucleotides are mutated target site nucleotides. All of the genes were cloned into a binary vector under the control of the 35S promoter. A YFP-only construct (i) was used as a negative control. The MpC3HDZ1-YFP (ii), MpC3HDZ1^mut-YFP (iii), MpAGO1-YFP (iv), MpMADS1-YFP (v), MpMADS1^mut-YFP (vi) and MpMADS2-YFP (vii) contructs contain 120 nt of the 5′ end of the gene; thus, they contain the Mpo-miR166a, Mpo-miR11707.1, Mpo-miR11687.1 and Mpo-miR11681.1 target sites, respectively, fused with the YFP gene. (B) YFP reporter assays were assessed using confocal microscopy. Scale bar = 250 µm. (C) The relative expression of the YFP reporter with/without miRNAs. The miRNAs were detected by Northern blot. 5S rRNA and tRNAs were used as loading controls. (D) The relative expression of MpMADS-YFP genes with/without miRNAs. The miRNAs were detected by stem–loop real-time RT–PCR, and miRNA levels were normalized to endogenous miR159 expression in N. benthamiana. Relative expression of the targets was normalized to the expression of NbActin as assessed by real-time RT–PCR. Bars represent the SEs (n = 3). The relative expression of the target in the presence of a miRNA precursor is significantly different from that in the presence of the target only (without miRNA treatment) for each RNA sample, according to the results of Student’s t-test; *P < 0.05; **P < 0.01. (E) Relative expression of wild-type MpMADS1 and the mutagenized miRNA-resistant MpMADS1^mut in proMpEF1α:MpMADS1 and proMpEF1α:MpMADS1^mut M. polymorpha plants, respectively. MpMADS1 expression in M. polymorpha transformed with an empty vector was set to ‘1’ and served as the wild-type (WT) reference. The numbers below indicate three independent transgenic M. polymorpha lines. Bars indicate the mean ± SE (n = 3).

Fig. 7

The MpPPR_66 (LW9825) transcript generates Mpo-miR11692.1 to target its mRNA and other targets. (A) The secondary structure of the MpPPR_66 transcript, and the hairpin structure for the generation of Mpx-miR15.1. (B) Paired regions of Mpo-miR11692.1/MpPPR_66 (i) and Mpo-miR11692.1/LW573 (ii). The dashed line indicates the 10th and 11th positions of the miRNA. The arrowhead and number indicate the position observed in the degradome and the number of read counts, respectively. (C) Degradome map of the miRNA target. The black line indicates significant (binomial test P-value < 0.001) degradome read counts at the 10th and 11th position of the target site. (D) Reporter assay to monitor Mpo-miR11692.1/MpPPR_66 (i) and Mpo-miR11692.1/LW573 (ii) interactions following transient expression in Nicotiana benthamiana. The relative expression of the targets was normalized to the expression of NbActin determined by real-time RT–PCR. Bars represent the SEs (n = 3). The relative expression of the targets in the presence of a miRNA precursor was significantly different from the expression in the presence of the target only (without miRNA treatment) for each RNA sample, according to the results of Student’s t-test; *P < 0.05; **P < 0.01. The miRNAs were detected by Northern blot. 5S rRNA and tRNAs were used as loading controls.