A genome-wide analysis of the RNA-guided silencing pathway in coffee reveals insights into its regulatory mechanisms

Abstract

microRNAs (miRNAs) are derived from self-complementary hairpin structures, while small-interfering RNAs (siRNAs) are derived from double-stranded RNA (dsRNA) or hairpin precursors. The core mechanism of sRNA production involves DICER-like (DCL) in processing the smallRNAs (sRNAs) and ARGONAUTE (AGO) as effectors of silencing, and siRNA biogenesis also involves action of RNA-Dependent RNA Polymerase (RDR), Pol IV and Pol V in biogenesis. Several other proteins interact with the core proteins to guide sRNA biogenesis, action, and turnover. We aimed to unravel the components and functions of the RNA-guided silencing pathway in a non-model plant species of worldwide economic relevance. The sRNA-guided silencing complex members have been identified in the Coffea canephora genome, and they have been characterized at the structural, functional, and evolutionary levels by computational analyses. Eleven AGO proteins, nine DCL proteins (which include a DCL1-like protein that was not previously annotated), and eight RDR proteins were identified. Another 48 proteins implicated in smallRNA (sRNA) pathways were also identified. Furthermore, we identified 235 miRNA precursors and 317 mature miRNAs from 113 MIR families, and we characterized ccp-MIR156, ccp-MIR172, and ccp-MIR390. Target prediction and gene ontology analyses of 2239 putative targets showed that significant pathways in coffee are targeted by miRNAs. We provide evidence of the expansion of the loci related to sRNA pathways, insights into the activities of these proteins by domain and catalytic site analyses, and gene expression analysis. The number of MIR loci and their targeted pathways highlight the importance of miRNAs in coffee. We identified several roles of sRNAs in C. canephora, which offers substantial insight into better understanding the transcriptional and post-transcriptional regulation of this major crop.

Fig 1. Phylogenetic tree of DCL-like proteins identified in Coffea canephora.

Phylogenetic tree showing relationships between the paralogous and orthologs proteins of the DCL family. The evolutionary history was inferred using the Neighbor-Joining method [46]. The bootstrap consensus tree inferred from 2000 replicates is taken to represent the evolutionary history of the taxa analyzed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (2000 replicates) are shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the JTT matrix-based method and are in the units of the number of amino acid substitutions per site [48]. The analysis involved 33 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 286 positions in the final dataset.

Fig 2. Analysis of the catalytic residues of the CcDCL-like proteins.

The two RNase III domains (RIBOc I and II) at the glutamate (E), aspartate (D), glutamate (D), aspartate (E) (EDDE) position. The catalytic sites are highlighted.

Fig 3. Phylogenetic tree of AGO proteins identified in Coffea canephora.

Phylogenetic tree showing relationships between the paralogous and orthologs proteins of the AGO family. The evolutionary history was inferred using the Neighbor-Joining method [46]. The bootstrap consensus tree inferred from 2000 replicates is taken to represent the evolutionary history of the taxa analyzed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (2000 replicates) are shown next to the branches. The evolutionary distances were computed using the JTT matrix-based method and are in the units of the number of amino acid substitutions per site [48]. The analysis involved 55 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 333 positions in the final dataset.

Fig 4. Analysis of the DxDGD catalytic motif of the RNA-dependent RNA polymerase (RdRP) conserved domain.

Six coffee RDR possess a DLDGD motif (CcRDR1.1–1.4, CcRDR2 and CcRDR6) and two have the DFDGD motif (CcRDR3.1 and CcRDR3.2), corresponding to the RDRα clade and the RDRγ clade, respectively (Blue box). Additionally, the RDRα displays a conserved subsequences (C/A)SG(S/G) before the DLDGD motif and, all CcRDR1 and the CcRDR2 showed the CSGS sequence, while CcRDR6 showed the ASGS sequence (red box).

Fig 5. Phylogenetic tree of RDR proteins identified in C. canephora.

Phylogenetic tree showing relationships between the paralogous and orthologs proteins of the RDR family. The evolutionary history was inferred using the Neighbor-Joining method [46]. The bootstrap consensus tree inferred from 2000 replicates is taken to represent the evolutionary history of the taxa analyzed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (2000 replicates) are shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the JTT matrix-based method and are in the units of the number of amino acid substitutions per site [48]. The analysis involved 33 amino acid sequences. All positions containing gaps and missing data were eliminated. There were a total of 312 positions in the final dataset.

Fig 6. Validation of the main proteins of genes involved in the coffee RNA-guided silencing pathways from RNAseq libraries.

Heatmap showing the expression pattern of the C. canephora RNA-silencing genes in three leaf developmental stages—Young, Expandend (“exp” in the figure), and Old—and Stem. (Transcriptome available at https://www.ncbi.nlm.nih.gov/sra/?term=ERP003741).

Fig 7

Alignment of pre-miRNA sequences (a), comparison of secondary structures (b) and phylogenetic tree (c) of ccp-MIR156 miRNAs and their orthologues. ccp- Coffea canephora, ath–Arabidopsis thaliana, nta–Nicotiana tabacum, mtr–Medicago truncatula, gma–Glycine max, mes–Manihot esculenta, ppe–Prunus persica, mdm–Malus domestica, vvi–Vitis vinifera, tcc—Theobroma cacao, ptc–Populus trichocarpa, aly–Arabidopsis lyrata, sly–Solanum lycopersicum. The evolutionary history was inferred using the Neighbor-Joining method[46]. The bootstrap consensus tree inferred from 5000 replicates is taken to represent the evolutionary history of the taxa analyzed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (5000 replicates) are shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Kimura 2-parameter method [3] and are in the units of the number of base substitutions per site[47]. The analysis involved 23 nucleotide sequences. All positions containing gaps and missing data were eliminated. There were a total of 68 positions in the final dataset.

Fig 8

Alignment of pre-miRNA sequences (a), comparison of secondary structures (b) and phylogenetic tree (c) of ccp-MIR172 miRNAs and their orthologues. ccp- Coffea canephora, ath–Arabidopsis thaliana, cme—Cucumis melo, gma–Glycine max, lus—Linum usitatissimum, mtr–Medicago truncatula, vvi–Vitis vinifera, bra–Brassica rapa, stu–Solanum tuberosum, nta–Nicotiana tabacum, aly–Arabidopsis lyrata, mdm–Malus domestica. The evolutionary history was inferred using the Neighbor-Joining method[46]. The bootstrap consensus tree inferred from 5000 replicates is taken to represent the evolutionary history of the taxa analyzed. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (5000 replicates) are shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site[47]. The analysis involved 28 nucleotide sequences. All positions containing gaps and missing data were eliminated. There were a total of 46 positions in the final dataset.

Fig 9

Alignment of pre-miRNA sequences (a), comparison of secondary structures (b) and phylogenetic tree (c) of ccp-MIR390 miRNAs and their orthologues. ccp- Coffea canephora, aly–Arabidopsis lyrata, ath–Arabidopsis thaliana, bna—Brassica napus, gma–Glycine max, ptc–Populus trichocarpa, sly–Solanum lycopersicum, mdm–Malus domestica, tcc—Theobroma cacao, lus—Linum usitatissimum. The evolutionary history was inferred using the Neighbor-Joining method[46]. The optimal tree with the sum of branch length = 1.87754489 is shown. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (5000 replicates) are shown next to the branches. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Kimura 2-parameter method and are in the units of the number of base substitutions per site[47]. The analysis involved 22 nucleotide sequences. All positions containing gaps and missing data were eliminated. There were a total of 65 positions in the final dataset.

Fig 10. SEA (Singular Enrichment Analysis) of the GO terms of the predicted targets of the ccp-miRNAs.

(A) Biological process, (B) molecular function and (C) cellular component.