Metagenomics of African Empogona and Tricalysia (Rubiaceae) reveals the presence of leaf endophytes

Plant material, DNA isolation, and sequencing

This study investigates four species of the Coffeeae tribe: Empogona congesta (Oliv.) J.E.Burrows, Tricalysia hensii De Wild., T. lasiodelphys (K.Schum. & K.Krause) A.Chev., and T. semidecidua Bridson (Table 1). These species were included in a previous study about chloroplast genome evolution in Rubiaceae (Ly et al., 2020). The plant material was obtained from the collection of Meise Botanic Garden, Belgium. Total DNA isolation from the leaves and DNA sequencing was done as described in Ly et al. (2020). Raw Illumina reads (BGI-seq 500 platform, 2 × 100 bp paired-end) are available under the BioProject PRJNA880288 (Table 1). While processing the raw sequencing reads, Ly et al. (2020) encountered “contamination” (i.e., non-chloroplast reads) and removed those reads to be able to reconstruct the chloroplast genomes of these four species. However, in this study, we are interested in this “contamination” in the raw reads because it contains information on possible leaf endophytes.

Read filtering and assembly

The contaminants in the Illumina reads were explored using metagenomic analysis tools. First, Kaiju v.1.9 (Menzel, Ng & Krogh, 2016) was used to identify and classify raw reads not belonging to the plant genome with the NCBI non-redundant RefSeq protein database (NCBI nr_euk). When a significant number of contaminants was present, it was examined in more detail in the Kaiju output. Second, the presence of contaminants was confirmed by exploring the distribution of GC count per sample using FASTQC: several GC count peaks in a single sample might suggest the presence of different organisms in the reads. Finally, reads from contaminants showing different GC count were filtered out using KAT v.2.3.4 (k-mer Analysis Tool; Mapleson et al., 2017). The reads were analysed with KAT gcp to create a matrix of the number of k-mers found, given k-mer frequency (27-mer) with GC count for each distinct k-mer to explore GC bias. The matrix was displayed via a density plot of the k-mer coverage versus GC count. KAT filter tools were used to filter out reads according to the GC bias (k-mer coverage of 100 to 500X and GC count of 10 to 22%). The reads left after filtering (without trimming, cleaning, or error correction) were subsequently assembled using MaSuRCA v.3.2.6 (Zimin et al., 2013) into scaffolds (Table 2) with the default parameters. The draft genome assemblies of the four endophytes are available on GenBank (Table 3) and on Zenodo (Verstraete et al., 2022a). Additionally, the reads were assembled using metaSPAdes v.3.15.5 (Nurk et al., 2017) and those draft assemblies are also available on Zenodo (Verstraete et al., 2023).

Analysis of the assembled bacterial draft genomes

The scaffolds obtained after assembly were compared to a Burkholderia s.l. sequence database of 2,288 genomes (representing 22 Gb) downloaded from NCBI (https://www.ncbi.nlm.nih.gov/assembly/?term=burkholderia) as available in September 2019. BLASTn v.2.2.26 (NCBI BLAST) was used for the comparison. Scaffolds with an e-value <10 e⁻²⁰ were kept. Assembly completeness was assessed using BUSCO v.5.4.3 (Seppey, Manni & Zdobnov, 2019) with the proteobacteria_odb10 database downloaded from https://busco.ezlab.org. The assembled draft genomes were annotated using Prokka v.1.14.5 (Seemann, 2014) and the annotations are available on Zenodo (Verstraete et al., 2023).

Assembly of the 16S rRNA genes

A custom pipeline was developed to assemble targeted genes from the raw Illumina reads (Mendez Silva et al., unpublished). In short, raw Illumina reads were mapped to a 16S rRNA reference gene (CP000010.1: c2677815–2676328 Burkholderia mallei) using Bowtie2 v.2.4.4 (Langmead & Salzberg, 2012). Mapped reads were subsequently filtered out and assembled with ABySS v.2.2.1 (Jackman et al., 2017).

Phylogenetic analysis of the endophytes

Three genes (16 rRNA, gyrB, and recA) were identified and used for phylogenetic analysis. A FASTA file with these sequences is available on Zenodo (Verstraete et al., 2022b). The 16S rRNA sequences were obtained from the raw Illumina reads, while the gyrB and recA sequences were obtained from the assembled scaffolds using the BLASTn tool. The gyrB nucleotide sequence (NC_006348.1: 3081–5549 Burkholderia mallei) and the recA nucleotide sequence (NC_006348.1: 290127–291197 Burkholderia mallei) were used as queries. Sequences were extracted using the extractseq function as implemented in EMBOSS (Rice, Longden & Bleasby, 2000).

The obtained sequences were combined with previously published datasets (Lemaire et al., 2011b; Lemaire et al., 2012b; Verstraete et al., 2013b; Danneels et al., 2023) to assess the phylogenetic position of the detected endophytes (Table S1). Automatic sequence alignment was performed with MAFFT v.7.490 (Katoh et al., 2002), followed by manual optimisation in Geneious R11. Possible incongruence among the different datasets was tested using a partition homogeneity test (implemented in PAUP*4.0b10a; Swofford, 2003). Due to sensitivity issues with the latter test (Barker & Lutzoni, 2002), resolution and support values of the different topologies were visually examined (hard versus soft incongruence; Johnson & Soltis, 1998). The best-fit nucleotide substitution model for each gene marker was selected using the Akaike information criterion in jModelTest v.2.1.10 (Posada, 2008). The model selection test showed that the GTR+I+G model is the most optimal model for 16S rRNA, and that the GTR+G model is the best for gyrB and recA. Bayesian inference analyses were performed with MrBayes v.3.2.7 (Huelsenbeck & Ronquist, 2001) on three individual data partitions and a combined data matrix under a mixed-model approach. Ten million generations were run, and parameters and trees were sampled every 1,000th generation. Chain convergence and ESS parameters were checked with Tracer v.1.7.2 (Rambaut et al., 2018). Bayesian posterior probability values above or equal to 0.95 are regarded as statistically supported (Alfaro, Zoller & Lutzoni, 2003).

Read filtering and assembly

Several species of Rubiaceae were recently sequenced using a whole genome sequencing approach with DNA extracted from leaves to establish their phylogenetic relationships (Ly et al., 2020). During the quality analysis steps, two types of raw reads with very different GC count per read (one peak at 39% and a second at 63%) were found, suggesting the probable presence of multiple organisms in the raw reads. To understand the origin of the different GC count peaks, a metagenomic approach was applied on the raw reads from Empogona and Tricalysia species.

For E. congesta, 37% of all reads could be assigned to a taxon name based on the NCBI non-redundant RefSeq protein database. About 30% of all reads or 80% of the named reads is assigned to Burkholderiaceae (Fig. S1A). For T. semidecidua, 30% of all reads is assigned to a taxon name and about 25% of all reads or 83% of the named reads is assigned to Burkholderiaceae (Fig. S1C). For T. lasiodelphys, 23% of all reads is assigned to a taxon name and about 18% of all reads or 75% of the named reads is assigned to Burkholderiaceae (Fig. S1E). For T. hensii, 20% of all reads is assigned to a taxon name and about 14% of all reads or 73% of the named reads is assigned to Burkholderiaceae (Fig. S1G). The Kaiju output can also be found in Data S1.

To confirm the presence of raw reads belonging to bacteria, the levels of k-mer coverage and GC count per distinct k-mer were analysed. The density plots of the k-mer coverage and the GC count indicate a low to medium k-mer coverage with a wide spread of the GC count (Figs. S1B, S1D, S1F and S1H, left part of the plots) suggesting the presence of reads with sequencing error and a medium coverage genome (i.e., the plant genome). However, the plots also show high coverage (approximately 700X) with GC counts of 15 to 20% (Fig. S1B, S1D, S1F and S1H, upper right of the plots). These unexpected reads with high coverage and high GC count were extracted from the set of raw reads using KAT.

The filtered bacterial reads from E. congesta were assembled into 632 scaffolds with an assembly size of 3.8 Mb, while the bacterial reads from T. hensii, T. lasiodelphys, and T. semidecidua were assembled into 769, 736, and 1,578 scaffolds with assembly sizes of 5.7, 7.8, and 4 Mb, respectively. The assembly N50 ranged from 9.8 to 48 Kb (Table 2).

Analysis of the assembled bacterial draft genomes

For the host species E. congesta, 612 of the 632 scaffolds (97%) resulted in a strong hit against the Burkholderia s.l. database (e-value cut-off 10 e⁻²⁰). For the draft genomes of the host species T. lasiodelphys, T. hensii, and T. semidecidua, there was a lower proportion of hits (89%, 50%, and 31%, respectively) (Table 3). Removing scaffolds without a strong hit against the Burkholderia s.l. genome database increased the N50 of the assemblies. However, this did not result in a large decrease in assembly sizes, indicating that only small-size scaffolds were discarded. BUSCO analysis indicated high completeness of the assemblies (>90%), except for the draft genome of the host species T. lasiodelphys (75.3%). In contrast to the other species, the large number of scaffolds and low N50 value for T. lasiodelphys suggest a fragmented and incomplete assembly. The genome assemblies should be considered as rough drafts, since bias might have been introduced when filtering the sequences by k-mer/GC count and BLAST. However, it is unlikely that host plant sequences remain present in the assemblies.

Phylogenetic analysis of leaf endophytes in Empogona and Tricalysia

The phylogenetic position of the endophytes found in Empogona congesta, Tricalysia hensii, T. lasiodelphys, and T. semidecidua were inferred from the 16S rRNA, gyrB, and recA sequences. The combined dataset demonstrated that all four non-nodulating leaf endophytes of Empogona and Tricalysia belong to Burkholderia s.l., more specifically, to the genus Caballeronia (Fig. 1). The non-nodulated endophyte of T. lasiodelphys is related to the nodulated endophytes of Sericanthe andongensis (Hiern) Robbr. and the non-nodulated endophytes of Psychotria psychotrioides (DC.) Roberty (BPP: 0.63). The non-nodulated endophyte of E. congesta falls within a highly supported clade of nodulated endophytes of several Pavetta species and non-nodulated endophytes of several Globulostylis species (BPP: 1.00). The non-nodulated endophyte of T. hensii is found as sister to the nodulated Candidatus Burkholderia kikwitensis (BPP: 1.00), nested within a clade of several other nodulated endophytes of Psychotria species. The non-nodulated endophyte of T. semidecidua falls in a clade with Caballeronia fortuita and C. novacaledonica (BPP: 0.97).

Detection of non-nodulated leaf endophytes in Empogona and Tricalysia

The majority of the studies on phylogenetic relationships within the Rubiaceae family to date has relied on phylogenetic approaches using individual nuclear or plastid DNA markers, or a combination of both (Wikström, Bremer & Rydin, 2020). However, phylogenomic approaches using more comprehensive amounts of genetic data are becoming more and more common, e.g., mitochondrial genomic data (Rydin, Wikström & Bremer, 2017), plastid genomes (Ly et al., 2020; Wikström, Bremer & Rydin, 2020), or a combination of hundreds of nuclear genes (Antonelli et al., 2021; Thureborn et al., 2022). The onset of the high-throughput sequencing methodology provides a novel tool to also detect leaf endophytes in Rubiaceae. High-throughput sequencing allows for the sequencing of total DNA, which is subsequently cleaned using bioinformatic filtering to only retain desired sequences, i.e., plant DNA sequences in most cases (e.g., Charr et al., 2020). For example, in the study of Ly et al. (2020), the objective was to obtain whole chloroplast genomes from 27 species in the Rubiaceae family. Before chloroplast genome assembly, the raw data was “checked in order to detect potential contamination”, with the unwanted reads being removed from further analysis. However, this contamination could be valuable on its own and possibly contain information on the presence of leaf endophytes. In fact, previous studies that used a DNA sequencing approach to detect leaf endophytes in Rubiaceae (e.g., Van Oevelen et al., 2001; Lemaire et al., 2011b; Verstraete et al., 2013a) also extracted total DNA but then eliminated the plant DNA by targeting bacterial DNA.

Our study is based on the unpublished raw data of Ly et al. (2020) (but made available here) and looks for evidence of leaf endophytes in the read contamination, specifically focusing on the genera Empogona and Tricalysia. These two genera belong to the Coffeeae tribe and are closely related to Sericanthe (Arriola et al., 2018), a genus known for its leaf nodulated symbiosis (Lemaire et al., 2011a). Unlike Sericanthe, Empogona and Tricalysia do not have visible nodules in their leaves; the detected leaf endophytes are therefore non-nodulated endophytes. In fact, by examining total DNA, we detected bacterial reads in E. congesta, T. hensii, T. lasiodelphys, and T. semidecidua (Table 3), indicating the presence of leaf endophytes.

Our reassessment of the original reads of the study of Ly et al. (2020) shows that we are dealing with leaf endophytes that are present in large proportions. Epiphytic contamination is unlikely because the leaves were cleaned with sterile water before the extraction of the total DNA. We have found that for each of the four investigated species, a large proportion of the reads is assigned to a limited taxonomic group. In T. semidecidua, this even reaches 83% of the named reads. This finding is in line with all previous studies about Rubiaceae endophytes. Although at this point, we cannot rule out that there might be more complex communities within the leaves, the fact remains that the particular Burkholderia-Rubiaceae interaction has been demonstrated for this new group of plants.

As such, the fact that leaf endophytes are detected in plants previously not implicated in leaf symbiosis is not unexpected. The wider occurrence of plant-bacteria interactions in Rubiaceae has been suggested before (Lemaire et al., 2012b; Verstraete et al., 2013a). It is therefore likely that additional hidden plant-bacteria interactions will be found when a systematic survey of leaf symbiosis in Rubiaceae is done.

The identity of leaf endophytes and their phylogenetic relationships

Our metagenomic analysis revealed that the bacterial leaf endophytes in Empogona congesta, Tricalysia hensii, T. lasiodelphys, and T. semidecidua belong to the family Burkholderiaceae. This is fully expected as so far, all leaf endophytes in Rubiaceae host plants have been identified as Burkholderia s.l. (e.g., Lemaire et al., 2012b; Verstraete, Janssens & Rønsted, 2017; Pinto-Carbó et al., 2018; Sinnesael, 2020; Georgiou et al., 2021).

After finding out the preliminary identity of the leaf endophytes, three genetic markers were extracted in order to achieve a more accurate identification, as well as to include the newly discovered leaf endophytes in a phylogenetic framework. Analysis of 16S rRNA, gyrB, and recA has already been extensively used to delineate species within Burkholderia s.l., as well as to unravel phylogenetic relationships at the generic level (Verstraete et al., 2011; Verstraete et al., 2013a; Lemaire et al., 2011a; Lemaire et al., 2011b). The use of these three markers particularly allows for comparison with other leaf endophytes and Caballeronia type strains. Even though the use of genomic data would be preferable and is common in recent literature about free-living Burkholderia s.l. (Bach et al., 2022), genomic information about leaf endophytes is often lacking. Such genomic data is usually extracted from pure cultures, but this is not possible for nodulated leaf endophytes, as they cannot be cultivated (Sinnesael et al., 2019). The study of the genomes of non-nodulated endophytes shows more promise but this has only just begun (Danneels et al., 2023).

The endophytes of Empogona and Tricalysia are identified as members of the genus Caballeronia (Fig. 1). The non-nodulated endophyte of T. lasiodelphys is most closely related to nodulated endophytes of Sericanthe and non-nodulated endophytes of Psychotria, while the endophyte of E. congesta is related to nodulated endophytes of Pavetta and non-nodulated endophytes of Globulostylis. The endophyte of T. hensii is most closely related to nodulated and non-nodulated endophytes of Psychotria. The endophyte of T. semidecidua falls in a clade with Caballeronia fortuita and C. novacaledonica. The type strain of C. fortuita was isolated from Fadogia homblei (Rubiaceae) rhizosphere soil in South Africa (Verstraete et al., 2014; Peeters et al., 2016), while the type strain of C. novacaledonica was isolated from Costularia (Cyperaceae) rhizosphere soil in New Caledonia (Guentas et al., 2016). Besides the fact that all these endophytes belong to the genus Caballeronia, there does not seem to be much of a phylogenetic pattern.

The study of Van Oevelen et al. (2001), which was the first to identify leaf endophytes in Rubiaceae host plants (i.e., in a few Psychotria species), found that the 16S rRNA sequences of the endophytes were similar to that of Burkholderia glathei. As a result, they assigned the Rubiaceae leaf endophytes to the genus Burkholderia (Van Oevelen et al., 2001). However, since then, several changes have been made to the taxonomy of this genus. First, all leaf endophytes were transferred to Paraburkholderia, when Burkholderia s.l. was split into a pathogenic group (Burkholderia s.s.) and a lineage of environmental bacteria (Paraburkholderia) (Sawana, Adeolu & Gupta, 2014). Later, Paraburkholderia was further subdivided and a new genus was created, Caballeronia (Dobritsa & Samadpour, 2016), which holds all nodulated endophytes as well as the non-nodulated endophytes of Globulostylis and Psychotria. The present study also designates the newly discovered non-nodulated endophytes of Empogona and Tricalysia as species of the genus Caballeronia (Fig. 1). The non-nodulated endophytes of the Vanguerieae genera (Fadogia, Fadogiella, Rytigynia, and Vangueria) remain in Paraburkholderia, except for those of Globulostylis (Data S2). This is in agreement with what was previously known (Verstraete et al., 2013b).

None of the investigated host plants has conspicuous bacterial leaf nodules in their leaves, and the endophytes are therefore non-nodulated endophytes. When looking at the phylogenetic tree of the leaf endophytes (Fig. 1), there is no apparent phylogenetic pattern for nodulation. The non-nodulated endophytes in Empogona and Tricalysia are not clustered, although they all belong to the genus Caballeronia. The newly found endophytes are related to other nodulated or non-nodulated leaf endophytes. Also, when analysing the results in a larger framework, no phylogenetic pattern is apparent for all non-nodulated endophytes in Rubiaceae since the majority of the Vanguerieae endophytes are situated within the genus Paraburkholderia (Data S2). However, we hypothesize that leaf nodulation should be considered as a character of the host plants, rather than of the leaf endophytes (see also Lemaire et al., 2012b).

Patterns in the host plants

In this study, we found Burkholderia s.l. endophytes in two genera of Rubiaceae that were previously not known to take part in leaf symbiosis. This brings the total number of genera in Rubiaceae for which leaf symbiosis is (molecularly) confirmed to ten: Psychotria (Psychotrieae tribe), Pavetta (Pavetteae tribe), Fadogia, Fadogiella, Globulostylis, Rytigynia, Vangueria (Vanguerieae tribe), and Empogona, Sericanthe, Tricalysia (Coffeeae tribe).

Finding phylogenetic patterns is however challenging. For the five Vanguerieae genera, the presence of Burkholderia s.l. endophytes is consistent at the genus level (Verstraete et al., 2013a) and the plants with leaf symbiosis only occur in Africa and Madagascar. For the genus Pavetta, it is generally assumed that most of the species have leaf nodules (Lersten & Horner, 1976; Miller, 1990; Lemaire et al., 2011b) and these are found in the Paleotropics (POWO, 2023). The presence or absence of nodules, as well as their form, has been used in the past to classify subgeneric taxa (e.g., P. series Enodulosae; Bremekamp, 1939). However, within Pavetta, the phylogenetic distribution of species without nodules is irregular (Bremekamp, 1934). Within the pantropical genus Psychotria, the situation is even more complex. The number of (known) species with nodules (ca 80 spp; Lemaire et al., 2012b) is rather limited compared to the total number of species in the genus (ca 1645; POWO, 2023, meaning that the nodulated form of leaf symbiosis is not a frequent character in Psychotria. Also, nodulating Psychotria plants are restricted to Africa and Madagascar. Unfortunately, a detailed list with the presence and absence of nodules is missing, so it remains uncertain to date to what extent leaf nodulation is present in Psychotria. The few nodulating Psychotria species that were included in molecular studies were not recovered as a monophyletic group (Razafimandimbison et al., 2014). Besides nodulating species, there are also some species without nodules but with non-nodulated endophytes (Lemaire et al., 2012b). Although these non-nodulating Psychotria plants were found in a clade (clade III in Lemaire et al., 2012b) separate from the nodulating species (clade II in Lemaire et al., 2012b), they also do not form a monophyletic group. Finally, the species Psychotria lucens Hiern was used in the past as a negative control (Van Oevelen et al., 2001) and later several other species without bacterial endophytes were found (Lemaire et al., 2012b). This means that all three conditions occur in Psychotria. However, the taxonomic delimitation of Psychotria has changed many times (Razafimandimbison et al., 2014) and finding evolutionary patterns within this megagenus remains challenging.

The genera Empogona, Sericanthe, and Tricalysia are closely related (Arriola et al., 2018), and the species of Empogona (Tosh et al., 2009) and Sericanthe (Robbrecht, 1978) were once included in Tricalysia. Perhaps it is therefore not so surprising to find leaf endophytes in these genera. The difference between Sericanthe on the one hand and Empogona and Tricalysia on the other, is that the former has leaf nodules, while the latter do not. A next step would be to investigate additional species of Empogona and Tricalysia to elucidate the true extent of leaf symbiosis in these two genera and to find out whether leaf symbiosis has any phylogenetic signal. Another observation worth investigating is that all tree genera are related to Diplospora (Arriola et al., 2018) and Discospermum Dalzell (Tosh et al., 2009). However, the major difference is that these two genera are found in (sub)tropical Asia, while the other three genera (Empogona, Sericanthe, and Tricalysia) are exclusively found in continental Africa and Madagascar (POWO, 2023). So far, the nodulated symbiosis is predominantly found in Africa (except for a few noduled Pavetta species in (sub)tropical Asia) and non-nodulated symbiosis is even restricted to that area. A broader screening of the Coffeeae tribe would therefore be useful to demonstrate the presence or absence of a geographic pattern in leaf symbiosis. However, for this, a new phylogenetic framework of the Coffeeae is needed, which shows the relationships between the different genera and onto which the character “leaf symbiosis” can then be plotted.