The complete organellar genomes of the entheogenic plant Psychotria viridis (Rubiaceae), a main component of the ayahuasca brew

Plant sampling, DNA extraction, and sequencing

Samples of P. viridis were collected in the Nucleo Menino Galante located at Serra da Cantareira (São Paulo) (23.380562 S, 46.58857 W) with the collaboration of União do Vegetal (approval number #23092021) and deposited at the Herbarium JABU from the Universidade Estadual Paulista (Unesp), campus Jaboticabal (voucher #3390). For the P. viridis genome sequencing, fresh leaves were collected after 24 h of dark treatment and submitted to flash freezing in liquid nitrogen, and stored at −80 °C in a super freezer. The plant samples were shipped to the Arizona Genomics Institute (Tucson, AZ, USA) where a high molecular weight (HMW) DNA extraction was performed using a modified CTAB protocol. The HMW DNA quality was checked for concentration on a Qubit dsDNA high-sensitivity assay, and absorbance ratios (260/280 and 260/230 nm) on a NanoDrop NP-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The DNA size was validated by a Femto Pulse and pulse-field gel electrophoresis, and the purity was validated using restriction enzyme digestion test. The HMW DNA was sheared to the appropriate size range (10–30 kb) using a Covaris g-TUBE, followed by bead purification with PB Beads (PacBio, Menlo Park, CA, USA) and libraries construction according to the manufacturer’s protocol using a SMRTbell Express Template Prep Kit 2.0. Sequencing was performed on the PacBio Sequel IIe (CCS mode – HiFi reads) platform at Arizona Genomics Institute (Tucson, AZ, USA).

Organelle genome assembly

We adopted a straightforward bioinformatics protocol to extract organellar-related reads from the sequenced reads. In the first step, we downloaded all ptDNA and mtDNA available sequences from GenBank Organelle Genome Resources (https://www.ncbi.nlm.nih.gov/genome/organelle/). For the ptDNA, the RefSeq fasta sequences were retrieved from the https://ftp.ncbi.nlm.nih.gov/refseq/release/plastid/ website, whereas the mtDNA RefSeq fasta sequences were obtained from the https://ftp.ncbi.nlm.nih.gov/refseq/release/mitochondrion/. In the second step, the fasta sequences of each RefSeq were concatenated and used to construct a custom Kraken2 (Wood, Lu & Langmead, 2019) database, according to the instructions provided at kraken2 GitHub and wiki (https://github.com/DerrickWood/kraken2/wiki/Manual#custom-databases). To facilitate the separation of the organellar HiFi reads from the sequencing reads, the ptDNA and mtDNA fasta sequences retrieved from the RefSeq had their “kraken:taxid” arbitrary set to “3702” for all ptDNA fasta sequences and “4530” for all mtDNA fasta sequences. In the third step, the Kraken v2.0 software was used to classify all sequenced reads using the custom database. Next, the “extract_kraken_reads.py” script from the KrakenTools (https://ccb.jhu.edu/software/krakentools/) was used to extract ptDNA and mtDNA related reads using the options:—include-children, —fastq-output, —taxid 3702 for ptDNA reads, and —taxid 4530 for mtDNA reads. Only organelle-related reads showing length >13 kb (N50 > 17 kb) and PHRED score >90 were used for the assembly. In the fourth step, each organelle was independently assembled with Flye v2.9 (Kolmogorov et al., 2019) using default parameters. In the final step, the “get_organelle_from_assembly.py” script from the GetOrganelle tool (Jin et al., 2020) was used to bait the ptDNA and mtDNA for each assembly graph generated (gfa file) using default parameters. The above-mentioned steps guaranteed a high-quality assembly showing average coverage of >500× for both organelles. The ptDNA and mtDNA candidate circular molecules were manually visualized and validated using the Bandage tool (Wick et al., 2015).

Organellar genome annotation, and comparative analyses

The ptDNA and mtDNA were annotated using the GeSeq platform (Tillich et al., 2017). For ptDNA, the Rubia cordifolia (NC_047470.1) and Coffea arabica (NC_008535.1) were used as references. For the mtDNA, the Apocynaceae mitogenomes Asclepias syriaca (NC_022796), Rhazya stricta (NC_024293.1), Cynanchum auriculatum (NC_041494.1) and the Magnoliaceae Liriodendron tulipifera (NC_021152.1) were used as references. All sequences were obtained from the GenBank Organelle Genome Resources. For the evaluation of the potential transfer of Transposable Elements (TE) to the mitogenomes, the The Extensive de novo TE Annotator (EDTA) package was employed for de novo TE detection (Ou et al., 2019). Relics of TE-related genes were predicted with TEsorter using the genome mode (Zhang et al., 2022). The final annotation was manually verified and refined according to the known plant mtDNA genes previously reported (Mower, Sloan & Alverson, 2012), and the organellar genome maps were made with the combination of the OGDraw (Greiner, Lehwark & Bock, 2019) and DNAPlotter (Carver et al., 2009) tools. The OGDraw and DNAPlotter figures were combined and manually edited with Inkscape (https://inkscape.org).

ptDNA comparative analyses were carried out with IRscope (Amiryousefi, Hyvönen & Poczai, 2018). Phylogenomics analyses were estimated using the maximum likelihood (ML) approach with a concatenated matrix generated by the “catfasta2phyml.pl” script (https://github.com/nylander/catfasta2phyml) and composed of 70 ptDNA genes (accD, atpA, atpB, atpE, atpF, atpH, atpI, ccsA, cemA, infA, matK, ndhA, ndhB, ndhC, ndhD, ndhE, ndhF, ndhG, ndhH, ndhI, ndhJ, ndhK, petA, petD, petG, petL, petN, psaA, psaB, psaC, psaI, psaJ, psbA, psbB, psbC, psbD, psbE, psbF, psbH, psbI, psbJ, psbK, psbL, psbM, psbN, psbT, psbZ, rbcL, rpl14, rpl16, rpl20, rpl22, rpl23, rpl32, rpl33, rpl36, rpoA, rpoB, rpoC1, rpoC2, rps2, rps3, rps4, rps7, rps8, rps11, rps14, rps15, rps16, rps18) shared across 39 Rubiaceae species (MK203877, MK203879, MN883829, MT674513, MT674517, MT674520, MT674521, MT674522, MW548283, MZ425928, MZ958829, NC_008535, NC_028009, NC_028614, NC_030053, NC_034698, NC_036300, NC_036970, NC_041149, NC_044102, NC_047302, NC_047470, NC_049078, NC_049155, NC_049160, NC_050962, NC_053701, NC_053762, NC_053818, NC_054151, NC_057496, NC_057533, NC_057593, NC_057983, NC_057984, NC_057985, NC_058252, OK236359, OK236360). The ML tree was calculated using IQ-Tree 2 (Minh et al., 2020) with the best-of-fit model GTR+F+R3, according to AIC criteria (Akaike, 1974), with the software ModelFinder (Kalyaanamoorthy et al., 2017). The clade support was estimated with the ultrafast bootstrap (UFBoot) and SH-aLRT algorithms (Hoang et al., 2018) with 1,000 replicates. Tree rooting was based on the Passiflora edulis (NC_034285), Vitis vinifera (NC_007957) and Arabidopsis thaliana (NC_000932) species as outgroup. mtDNA comparative analyses were carried out with Sibelia (Minkin et al., 2013). The long reads have been aligned to the assembled mitogenomes with pbmm2 (https://github.com/PacificBiosciences/pbmm2), and the coverage was calculated with the script “wgscoverageplotter.jar” from the Java utilities for Bioinformatics (JVARKIT) (https://github.com/lindenb/jvarkit). All figures were manually edited with Inkscape (https://inkscape.org).

Data availability

The P. viridis samples (GenBank BioSample SAMN26930853) were included in the Brazilian National System of Management of Genetic Heritage and Associated Traditional Knowledge (SisGen) under the access number #A3D39A2. The organellar genomes were deposited into GenBank, BioProject PRJNA819420, under the accession ON064099 for the ptDNA and ON064100 (mtDNA haplotype I) and ON968433 (mtDNA haplotype II) for the mtDNA.

Psychotria viridis ptDNA structure and main features

The P. viridis ptDNA genome span 154,106 bp, exhibiting a GC content of 37.55% and a typical quadripartite structure encoding all known ptDNA genetic repertoire found in photosynthetic organisms (Fig. 1 and Table 1). Phylogenomic analyses based on ptDNA, clustered all tribes as monophyletic groups and placed P. viridis in the Psychotrieae tribe, sister to Spermacoceae and Moriendae tribes (Fig. 2). The close relationship among these tribes was already suggested based on nuclear and plastidial markers (Löfstrand, Razafimandimbison & Rydin, 2019; Razafimandimbison et al., 2020, 2017). The main differences are related to Spermacoceae clustered with Psychotrieae, while in the previous studies, the Moriendae is more closely related to Psychotrieae than Spermacoceae. Most of the Moriendae species shown in the Fig. 2 are known as medicinal plants in Southeast Asia, particularly in traditional Chinese medicine. However, among the three close-related tribes, only P. viridis is capable to produce alkaloids and DMT, making this species particular in terms of genetic content.

Comparative ptDNA analysis evidenced that Psychotrieae, Spermacoceae, and Moriendae tribes are exhibiting slight junction site variation between the plastid regions (IR-SSC-IR and IR-LSC-IR) (Fig. 3). For instance, a conserved junction site pattern is observed between the two available Psychotria ptDNA, where P. rubra and P. viridis exhibit the ndhF and ycf1 genes delimiting the SSC borders (Fig. 3). The other junction site variation shown in Fig. 3 is common in other angiosperms and can be found in other plant clades (Amiryousefi, Hyvönen & Poczai, 2018; Chen et al., 2021; Li et al., 2021; Zhang et al., 2019).

Psychotria viridis mtDNA presents a complex structure, composed of two main circular chromosomes that are arranged into several isoforms and sub-circles

The plant mitochondrial genomes exhibits a complicated mixture of physical configurations and sequences, thereby showing variable and dynamic connections, such as branched linear structures containing multigenomic concatemers together with rare circular forms (Kozik et al., 2019). Recently, it was proposed that analyses based on the assembly graphs approaches are valid for visualizing the actual complexity of the plant mtDNA’s variable forms within the mitochondria organelles (Fischer et al., 2022). Graph-based analyses can also provide sufficient evidence to avoid the naïve and erroneous conclusion that a single circular mtDNA molecule exists (Fischer et al., 2022; Kozik et al., 2019).

The P. viridis mtDNA assembly graph presents two alternative and circular structures showing 44.37% of GC content, suggesting that heteroplasmy may be common in plants (Woloszynska, 2010). The mtDNA haplotype I spans 615,370 bp and contains one long (24 kb) repeat region, while the mtDNA haplotype II spans 570,344 bp and contains a smaller repeat region (740 bp) (Figs. 4A and 4B). Repeated regions can be the substrate for different homologous recombination or other recombination events that may symmetrically isomerize and rearrange the mitogenomes (Sloan, 2013). Therefore, these repeated regions may support the existence of a dynamic spectrum of up to 16 isoforms and up to eight different sub-circles arrangements for each mitogenome haplotype. All those variants probably coexist in vivo (Figs. 4C and 4D). Moreover, the 24 kb repeated region from the mtDNA haplotype I encodes for several mitochondrial genes, resulting in a redundant set of gene copies. This kind of complex structure is commonly found in other plant mitochondrial genomes and it is now considered more a rule than an exception that can result from repeat-mediated or other recombination events (Fischer et al., 2022; Kozik et al., 2019).

We further evaluated the macrosynteny patterns, revealing that six large collinear blocks are shared among the mtDNA haplotype I and haplotype II (Figs. 5A–5E), in addition to the internal and shared repeats (Figs. 5F and 5G). The average nucleotide identity between both mitogenomes haplotypes is 100% between the shared collinear blocks. Moreover, alignment analyses indicated that 86% of haplotype I is covered by haplotype II, while 92% of haplotype I is covered by haplotype II. Therefore, there are 85,631 and 45,965 bp of exclusive regions for haplotype I and haplotype II, respectively.

Psychotria viridis mtDNA gene content and main features

In general, plant mtDNA genomes carry approximately 40 protein-coding genes, a variable set of tRNAs, and the insertion of ptDNA-related sequences as a product of lateral endosymbiotic gene transfer (Mower, Sloan & Alverson, 2012). The P. viridis mtDNA follows these trends, encoding most genes typically found in mitochondrial genomes (Fig. 6, Tables 2 and 3). The main finding is related to the absence of the rps2, rps11, rpl2, rpl6, rpl10, and rpl16 genes, but the ribosomal genes are highly variable across the angiosperms (Mower, Sloan & Alverson, 2012). Additionally, the sdh3 (Uniprot: P0DKI0 and GO:0000104), which encodes the succinate dehydrogenase cytochrome b subunit, is absent in both versions of the P. viridis mtDNAs, but also other plants mitogenomes (i.e., Oryza sativa). On some plants like Arabidopsis thaliana, this gene has been integrated into the nuclear genome (Mower, Sloan & Alverson, 2012), supporting a potential function redundancy mechanism using nuclear homologous genes. Moreover, the nad5 gene, which is formed by five exons, presents a duplicated copy of exons 1 and 2. Since nad5 is processed by trans-splicing, this finding may suggest the existence of additional transcriptional paths responsible to join the spliceosome’s primary RNA transcripts into the nad5 mRNA assembly.

Other remarkable features are related to the differences between the two mtDNA haplotypes, including: the presence of a LAGLIDADG homing endonuclease located between only the cox1 intronic region of the mtDNA haplotype I, and an additional and chimeric version of nad4L (326 bp) identified in each of the mtDNA haplotypes. The LAGLIDADG homing endonuclease are cleaving enzymes that are generally encoded by introns, conferring intron mobility, and generally found in fungal mitochondrial genomes (Megarioti & Kouvelis, 2020). The presence of the LAGLIDADG in only mtDNA haplotype I suggests horizontal transfer from a fungal donor species, but still maintaining the cox1 loci intact in the mtDNA haplotype II. Both chimeric nad4L share only the first 16 amino acids, including the start codon, located on the N-terminus of the protein. Interestingly, a closer alignment inspection on the nad4L chimeric using the UniProt Blast (Pundir et al., 2015), reveals that, while the N terminus of the protein (including the 16 amino acids above mentioned) matched to different plant mitochondrial genome (i.e., Capsicum chinense, Nicotiana tabacum, Zea mays, Vitis vinifera, Glycine max, and many others), the rest of the protein exhibits matches to ortholog proteins from Fungi, Bacteria and Nematoda (UniProt IDs: A0A2G9TNM8, A0A1Y1ZXS8, R6LZY0, A0A238BKX1, and A0A1Y2DMI6); thereby, supporting a chimeric structure possibly acquired by lateral transference.

Additionally, the atp1 (1,776 bp) from haplotype II is 405 bp longer than its respective homolog located on haplotype I, indicating considerable differences at the protein C-terminal. Using the Arabidopsis thaliana atp1 (Uniprot: P92549 – 1,521 bp) as a reference, we observed 95–96% of similarity across almost the entire proteins (from the start codon to the position 1,350), supporting that P. viridis exhibits indeed two different C-terminal domains of Atp1 protein. Therefore, it is possible to speculate that these differences may impact the electron transport complexes of the respiratory chain and ATP production in a unique way for this species. This hypothesis will need further investigation.

Slight differences in the nad1 gene were also observed. For instance, the nad1 from haplotype II lacks three nucleotides (AAT); however, it does not impact the protein translation, making both genes 99% identical, showing only two amino acids difference.

Several ptDNA sequences representing complete and fragmented genes, and pseudogenes were found dispersed across the mitogenomes of P. viridis, accounting for up to 7, 9 kb in each mtDNA haplotype (Table 4). Noteworthy, we have found a chimeric gene formed by fusion of the chloroplast rpl2 and rpl23 genes. Furthermore, the group II intron splicing reverse transcriptase gene (ltrA), which is commonly observed in the ptDNA of some Angiosperms, was found as a fragment in P. viridis mtDNA (Table 4), and no mitovirus-derived sequences were identified in P. viridis mtDNA. Interestingly we have found a gene encoding a potential plant aconitate/2-methylaconitate hydratase, which is commonly involved with the tricarboxylic acid cycle, the second stage of cellular respiration, and the glyoxylate cycle. In addition, a plant glyceraldehyde-3-phosphate dehydrogenase involved in the glycolysis, and an elongation factor related to the protein synthesis were identified in both mitogenomes haplotypes, thus, suggesting acquisition from the nuclear genome.

Finally, genes encoding a ribonuclease H and an unknown polyprotein were identified in both mtDNA haplotypes. The presence of these genes is a clear indicator of a previous event of TEs migration from the nuclear to the mitochondrial genome. Indeed, several TEs scars were identified in P. viridis mtDNA, as revealed in Fig. 6, thus evidencing that these mobile elements helped to shape the structure of chacrona’s mitogenomes.