University of Groningen The ecology and evolution of bacteriophages of mycosphere-inhabiting Paraburkholderia spp. Pratama, Akbar Adjie

The ecology and evolution of bacteriophages of mycosphere-inhabiting Paraburkholderia spp.

Pratama, Akbar Adjie

IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.

Document Version

Publisher's PDF, also known as Version of record

Publication date: 2018

Link to publication in University of Groningen/UMCG research database

Citation for published version (APA):

Pratama, A. A. (2018). The ecology and evolution of bacteriophages of mycosphere-inhabiting Paraburkholderia spp. Rijksuniversiteit Groningen.

Copyright

Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).

Take-down policy

If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.

Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.

The dynamics and diversity of the

mycosphere virome

Akbar Adjie Pratama, Xiu Jia, Juliana Eschholz and Jan Dirk van Elsas In preparation for publication

Abstract

Viruses are abundant entities in most biomes on Earth and their role in ecosystem processes is essential. However, the viruses in soil ecosystems have been overlooked due to the challenging sampling and cultivation approaches. The current next-generation sequencing techniques have revolutionized the study of environmental viruses. Here, I hypothesize that soils contain high viral diversities, especially the mycosphere. The ecological role of viruses for mycosphere bacterial communities is, however, underexplored. In this study, I found that the mycospheres of two fungi, i.e. Russula ochroleuca and R. emetica, contain as-yet-undescribed viral communities. Most viruses found were predicted to infect Proteobacteria, including Paraburkholderia and Burkholderia spp. The analysis of the genomes of selected viruses showed that phages related to ϕ437/VC20 were abundant across two mycospheres. These phages, infecting mycosphere Paraburkholderia spp., were found to contain a gene, amrZ, that may encode a protein involved in modulation of biofilm formation. On another notice, a contig representing another phage, denoted VC14, from bulk soil populations, and predicted to infect Paraburkholderia spp., was found to encompass a gene encoding a phasin; this protein may advance host growth in nutrient-limited environments. Keywords: Soil, Mycosphere, virome, ecology

8

Introduction

Soil provides an environment in which organisms can live either apart or in close proximity provided colonisable surfaces are present. In the latter situation, they may establish interactions that can be manifold. The mycosphere provides a microenvironment in soil that offers both a colonisable surface and a ready abundance of nutrients – these factors accelerate microbial functioning, spur evolutionary adaptation and foster genetic exchange of its inhabitant organisms (Zhang et al., 2014). Paraburkholderia spp. - as mycosphere inhabitants - are frequently exposed to mobile genetic elements (MGEs) such as plasmids and phages. Thus, these organisms may serve their versatile lifestyle with genetic traits that favour their survival, for example, fungal- and plant- interactive traits (Haq et al., 2014; Warmink et al., 2011). Compared to marine and human gut ecosystems (Roux et al., 2016; Yutin et al., 2017), the role of soil viruses, especially of the mycosphere, has remained rather unexplored (Pratama and Elsas, 2017). The main reasons for this lack of attention and understanding lie in the physicochemical nature of the soil itself, the challenges posed with respect to required technology, the lack of spatiotemporal data and that of understanding bacterial-virus interactions at the microscale (Pratama and Elsas, 2017).

Metagenomics of the virome can accelerate the pace of virus discovery and thus provide the necessary foundation for the meaningful integration of viruses into ecosystem function theories (Brum and Sullivan, 2015). However, there is as-yet no study that aims to identify the virome in soil hotspots such as the mycosphere as related to key mycosphere inhabitants, e.g. Paraburkholderia. A recent study (Knowles et al., 2016) indicated that temperate phage dynamics is important in ecosystems with high microbial densities. Here, we explored the virome of the mycospheres of two selected fungi, Russula ochroleuca and R. emetica, versus the corresponding bulk soil, and hypothesized that (i) the mycosphere contains a high microbial diversity and an unexplored viral community (ii) mitomycin C (MMC) induction will excise phages from hosts and thus shift bacterial community compostions and (iii) phage communities have a significant role in mycosphere microbiomes. We set up a simple experiment, in which we selectively cultured natural Paraburkholderia communities in a semi-selective medium, PCAT. We then induced these communities to release phage populations, and subsequently examined these by transmission electron microscopy (TEM) and sequencing. The results of this study are presented here.

Material and Methods

Sampling and processing of mycosphere and bulk soil samples

All samples were taken in October 2016 across a forest soil area in Noordlaren, Drenthe, Netherlands. Samples of Russula spp. were taken based on the procedures outlined in Zhang et al. (Zhang et al., 2015). In total, 37 samples were collected from two different Russula spp., including 12 individuals of yellow mushroom (MY; Russula ochroleuca) and 18 of red mushroom (MR; Russula emetica). Complete mushroom bodies including soil adhering to the hyphal feet were collected. Moreover, seven bulk forest soil samples (B) were obtained from sites approximately 0.5-1 m away from the individual mushrooms. The 12 individuals of MY, 18 of MR and seven of B were each randomly pooled into four replicates per sample type (MY: 3 samples per pool; MR: 6 samples per pool; B: 1 or 2 samples per pool). Plate counts (using R2A and PCAT media) yielding colony-forming unit (CFU) numbers were used to estimate the viable bacterial abundances in each sample.

From each replicate of each sample type, (mycosphere or bulk) soil was obtained. Thus 0.3-g (mycosphere soil) and 0.5-g (bulk soil) portions were obtained, which were transferred to sterile tubes containing 6 and 10 ml of sterile 0.85% NaCl. After vigorous shaking for three times (1 min each time), the samples were centrifuged at 100xg for one min and left to settle for 10 min. Supernatant were then obtained, serially

diluted (10-1_{, 10}-2_{, 10}-3_{, 10}-4_{, 10}-5_{) and subsequently plated onto the “Burkholderia}

semi-selective” medium PCAT (azelaic acid 2 g/L, tryptamine 0.2 g/L, MgSO₄.7H₂O

0.1 g/L, KH₂PO₄ 4g/L, K₂HPO₄ 4g/L, yeast extract 0.02 g/L, agar 15 g/L; pH 5.7) (Salles

et al, 2006), supplemented with cycloheximide to 50 ug/mL. The 100_{, 10}-1 _{and 10}-2

dilutions were used. The other dilutions (10-3_{, 10}-4_{, 10}-5_{) were plated onto R2A medium}

(Beckton Dickinson, NJ, USA). Plates were incubated at 27o_{C for up to 13 days.}

Colony-forming unit (CFU) counts were obtained at regular time intervals, after which the size of the culturable bacterial communities per g dry-weight soil was determined.

Bacterial 16S rRNA gene sequencing and analysis

Nucleic acid was extracted from all samples using a modified (Powersoil) DNA isolation kit (MOBio Laboratories Inc., Carlsbad, CA, USA). The modification consisted of adding glass beads (200-500 μL) to the cultures to spur mechanical cell lysis. This extraction method is a rapid way to produce highly pure DNA from bacterial cultures. Moreover, DNA was extracted from the cultured Paraburkholderia communities obtained from homogenized PCAT plate colonies using the UltraClean Microbial DNA Isolation Kit (MOBio Laboratories Inc., Carlsbad, CA, USA). The extracted DNAs were

8

The 16S rRNA amplicon library preparation was prepared using barcoded primer

set adapted for the Illumina HiSeq2000 and MiSeq (Caporaso et al., 2012). Each 25 µL PCR reaction contained 9.5 µL of MO BIO PCR Water (Certified DNA-Free), 12.5 µL of QuantaBio’s AccuStart II PCR ToughMix (2x concentration, 1x final), 1 µL Golay barcode tagged Forward Primer (5 µM concentration, 200 pM final), 1 µL Reverse Primer (5 µM concentration, 200 pM final), and 1 µL of template DNA. The conditions for PCR are as follows: 94 °C for 3 minutes to denature the DNA, with 35 cycles at 94 °C for 45 s, 50 °C for 60 s, and 72 °C for 90 s; with a final extension of 10 min at 72 °C to ensure complete amplification. Amplicons are then quantified using PicoGreen (Invitrogen) and a plate reader (InfiniteÒ 200 PRO, Tecan). Once quantified, volumes of each of the products were pooled into a single tube so that each amplicon is represented in equimolar amounts. This pool is then cleaned up using AMPure XP Beads (Beckman Coulter), and then quantified using a fluorometer (Qubit, Invitrogen). After quantification, the molarity of the pool is determined and diluted down to 2 nM, denatured, and then diluted to a final concentration of 6.75 pM with a 10% PhiX spike for sequencing on the Illumina MiSeq. Amplicons are sequenced on a 151bp x 12bp x 151bp MiSeq run using customized sequencing primers and procedures (Caporaso et al., 2011, 2012). DNA sequence data were then generated using Illumina paired-end sequencing at the Environmental Sample Preparation and Sequencing Facility (ESPSF) at Argonne National Laboratory (Lemont, IL, USA). Specifically, the V4 region of the 16S rRNA gene (515F-806R) is PCR amplified with region-specific primers that include sequencer adapter sequences used in the Illumina flow cell (Caporaso et al., 2011, 2012). The analysis of 16S rRNA was performed with QIIME2 (in the Peregrine platform). Sequences were quality filtered at fasqc program and demultiplexed. Sequences were assigned to closed reference operational taxonomic units (OTUs) at a 97% identity threshold using the SILVA database (https://www.arb-silva.de). Alpha diversity OTU richness and UniFrac distance were calculated using the vegan R package using 500 permutations.

Mitomycin C (MMC) induction

Induction with mitomycin C (MMC) (Sigma-Aldrich, St. Louis, Mo, USA) was conducted according to Pratama and van Elsas, (2017). Briefly, the homogenized bacterial communities from the PCAT plates were introduced into 5 ml of LB medium and incubated overnight at 28 °C (shaking at 180 rpm). The resulting cultures were then transferred (1:100) into replicate Erlenmeyer flasks containing 40 ml of fresh LB medium and growth was monitored until the exponential growth phase was reached (after about eight h). Thereafter, cultures were split into two 20-ml cultures. MMC was added to the cultures, at final concentrations of 10 μg/mL (MMC-10), with the

corresponding culture serving as the control. The cultures were incubated and the

OD₆₀₀ was monitored for 24 h. Decreases of OD values were taken as indications of

progressive cell lysis and prophage release. The experiments were done with four biological replicates. The resulting crude lysates were finally filtered over Whatman 0.22 μm cellulose acetate filter (GE Healthcare Life Sciences, Pittsburgh, PA, USA) and stored at −20 °C until further analysis.

Phage precipitation, DNA extraction and sequencing

Viral particles were physically separated from bacterial culture using previously described methods (Pratama and van Elsas, 2017). Briefly, the suspensions produced as above were purified according to the poly ethylene glycol (PEG) method of

(Sambrook and Russell, 2001)with the following modifications as in (Pratama and

van Elsas, 2017). Suspensions were centrifuged at 11,000 xg for 15 min at 4 °C, and then supernatants were filtered over Whatman 0.22 μm puradisc syringe filter-cellulose acetate (GE Healthcare Life Sciences, Pittsburgh, PA, USA). NaCl (29.2 g) was dissolved in the 500 mL lysates to final concentration 1 M, which were then stored on ice for 1h. Solid polyethylene glycol (PEG) 8000 was added to the supernatant to a final concentration of 10% (w/v) and the mixture stored overnight at 4 °C to allow phage particles to precipitate. The PEG-precipitated lysate was then centrifuged at 11,000 xg for 10 min at 4 °C (Sorvall SLA-1500 rotor). The supernatants were

discarded to 20 mL and 10x SM buffer (10 mM NaCl, 50 mM Tris, 10 mM MgSO₄, and

0.1% gelatin) was added for storage and later analysis.

The virome nucleic acid extraction was performed with a virome DNA Isolation Kit (Norgen, Biotek Corp, ON, Canada) following manufacturer’s protocols, with slight modification, i.e. DNase I inactivation temperature was 80 °C for 10 min. In addition, 16S rRNA PCR amplification using 16SFP/16SRP universal 16S rRNA gene primer set

(Pereira e Silva et al., 2012)was performed to confirm the absence of genomic DNA

in the phage DNA extracts.

The Illumina HiSeq was used to sequenced virome DNA. 2500 paired-end by BaseClear (Leiden, Netherlands). The libraries were prepared using Illumina genomic Nextera XT libraries. The quality analyses of FASTQ sequence reads were done using the Illumina Casava pipeline version 1.8.3. The Initial quality assessment was based on data passing the Illumina Chastity filtering. Subsequently, reads containing PhiX control signal were removed using an in-house filtering protocol. In addition, reads containing (partial) adapters were clipped (up to minimum read length of 50 bp). The second quality assessment was based on the remaining reads using the FASTQC quality control tool version 0.10.0. The final quality scores per sample yielded 450,342 -1,6743,484 -reads, at 38.21 average qualities.

8

Real-time PCR quantification

The 16S rRNA gene (region V5-V6) was amplified by qPCR on an ABI Prism 7300 Cycler (ThermoFisher Scientific, Waltham, Mass, USA). The qPCR components consisted of 12.5 µL SYBR Green (Applied Biosystems, CA, USA), 0.75 µL each primer (Eurogentec, Liège, Belgium), 6 µL PCR-water and 2 µL extracted DNA sample (5 ng/mL) as a template and primers 16SFP/16SRP (), in total 25 µL reaction volume. The amplification and detection of ɸ437 product were preformed using ABI 7300 (ThermoFisher Scientific, Waltham, Mass, USA) with qPCR reaction conditions,

denaturation cycling at 95o_{C for 27 second, annealing cycling at 62}o_{C for 1 min and}

elongation cycling at 72o_{C for 30 second. Standard curves were generated over five}

orders of magnitudes, e.g. from 103 _{to 10}7 _{copies of template, using plasmid containing}

specific regions. The qPCR efficiency was 102 %. Data were calculated as log copy numbers per gram of dry-weight soil.

De novo assembly, virome identification and

taxonomic assignment

Virome raw reads were quality-filtered using FastQC. The high-quality reads were assembled using meta-SPAdes and the assembled reads were quality-checked using meta-QUAST in iVirus environment (Bolduc et al., 2017) in Cyverse (http://www. cyverse.org/). Viral sequences identification were predicted in assembled virome contigs (>1 kb) using VirSorter (Roux et al., 2015). VirSorter detects viral signals or regions based on some criteria: (i) the presence of viral hallmark genes (e.g. viral major capsid, tail), (ii) enrichment of gene with the best hits against viral reference database, (iii) the depletion of PFAM affiliated genes, (iv) enrichment of genes with no hits in PFAM and viral reference database, (v) the enrichment in short genes and (vi) the depletion in strand switching. Virome reference was used and to increase the confident on the identified viral regions, decontamination step was applied. Thus, only sequences with the very accurate viral region (category 1 and 2) and prophage (category 4 and 5) were kept for further analyses. InterproScan (Mitchell et al., 2015) was used to annotate the identified viral contigs.

In silico analysis to identify potential hosts with tRNA

Presumptive identification of hosts using tRNA matches was based on (Paez-Espino et al., 2016). The sequences of the tRNA genes were identified in the contigs of the identified viral and prophage sequences using ARAGORN (Laslett and Canback, 2004). To determine the possible host, the tRNA sequences were then compared against the genomic tRNA database (http://gtrnadb.ucsc.edu) and non-redundant database in BLASTN, using cut-off values of 90% coverage and 80% identity.

Results

Total and culturable bacterial communities in the mycosphere

and bulk soil

The average values of total bacteria on R2A medium were on the order of 106_{, 10}7

and 105 _{CFU per g dry soil for MY, MR and B samples, respectively. On the other hand,}

the average values of bacteria growing on PCAT medium were about 106_{, 10}6 _and

102 _{per g dry soil for MY, MR and B, respectively. Thus, the CFU counts showed a}

significant increase of bacterial abundances in the mycosphere samples compared to bulk soil ones on both media (ANOVA, n=4, p<0.05) (Figure 8. 1A). Moreover, Paraburkholderia species types growing on PCAT may have been enriched by the tweo mycospheres, making up considerable aprts of the respective communities, versus very low parts (< 0.1%) in bulk soil.

In contrast, the 16S rRNA gene based qPCR of the bacterial communities from the

mycosphere and bulk soil samples were higher, with copy numbers ranging from 108

to 109 _{per g dry soil. These demonstrated a trend similar to the CFU counts, however}

with no significant differences between the mycosphere and bulk soil values (ANOVA, n=4, p>0.05) (Figure 8. 1B). Thus, whereas the total bacterial community sizes were not strongly different between mycosphere and bulk soil samples, they significantly exceeded the numbers found on plates..

We further used PCR-DGGE to evaluate the bacterial community structures from the mycosphere and bulk soil samples (Figure 8.1C and Figure 8. 1D). We observed higher numbers of bands in the bacterial communities in the mycosphere samples (MY and MR) than in the B samples. We observed 13-20 dominant bands in MY and MR samples versus about 10-15 in B ones. Interestingly, two to three unique bands were observed in the mycosphere samples (Figure 8. 1C). Principal components analysis (PCA) revealed close relatedness between all mycosphere populations, which were distantly related to the B samples (Figure 8. 1D). This suggested niche commonality and similar ecological drivers in both mycospheres as compared to the ones in the B samples.

Bacterial community dynamics and inducibility of the cultured

mycosphere bacteriomes

To assess the shift of microbial communities upon induction and inducibility of the virome, the cultures produced from mycospheres versus bulk soils were induced in mid-log phase with MMC. The two MMC concentrations, 5 and 10 μg/mL and three different induction times, 8h, 10h and 12h were tested. A significant decrease of the

8

40 20 0 PCO2 (19,9% of total variation) -20 -40

-60 -40 -20 0 20 40

PCO1 (35,5% of total variation)

d a b c dc d c a b a b Russula ochr oleuca

Russula emetica Bulk soil

70 80 90 100 89.2 84.2 83.9 87.9 90.1 94.7 79.8 70.3 75.5 Bulk_A

Bulk_B Bulk_C Bulk_D MR_B MR_A MY_D MR_D MY_B MY_A

A

C

B

D

Growth medium

logCFU per gram mycosphere or bulk soil 1

2 3 4 5 6 7 8 PC AT R2 A Ru ssu la oc hr ole uc a Ru ssu la em eti ca Bu lk so il a a ab ab abc Sample

logBacteria 16S rRNA copy number per gram mycosphere or bulk soil

6 8 10 Ru ss ul a oc hr ol eu ca Ru ss ul a em et ic a Bu lk so il a Figure 8.1. Bacte riome structures from mycospheres of Russula ochroleuca and R. emetica and bulk soil samples. ( A) log CFU per g dry weight of mycos phere and bulk soils, (B ) log 16S rRNA gene copy number per g dry weight of mycosphere and bulk soil. (C and D) Denaturating

8

Figure 8. 2 . ( A) Mitomycin C (MMC) induction of PC AT cultures. ( B and

C) average relative abundance

of bacterial

phyla

and genera,

respectively.

Estimated (

D) Shannon index and (

E) principal coordinate analysis (PCoA) of PC AT selected community. (F )

Visualization of TEM of phage particles from

MY

and MR

samples

No such drop in the OD values was observed in the zero-MMC treatment (un-induced) across all cultures. Thus MMC had a bacterial population-reducing effect.

Then, 16S rRNA amplicon sequencing was used to assess the dynamics of the PCAT-cultured bacterial communities upon MMC induction. For this, we extracted total bacterial DNAs from the initial PCAT cultures, before MMC induction (8h) and after induction (24h).

The raw sequencing generated totals of 3,134,466 high-quality reads of the V4 region of the 16S rRNA gene (515F-806R). Collectively, the OTUs across all our cultured samples fell in eight prevailing bacterial phyla. These included Acidobacteria, Actinobacteria, Bacteroidetes, Chloroflexi, Firmicutes, Planctomycetes, Proteobacteria and Verrucomicrobia, next to others.

In all samples, Proteobacteria was the dominant phylum, making up >97% of the total bacterial communities (Figure 8.2B and Figure 8.2C). In the initial PCAT cultures, Burkholderia/Paraburkholderia (±40-58%) was the most dominant bacterial genus. However, the community structures clearly shifted in the 8- and 24-h cultures. In the B samples, Pseudomonas started to dominate (≥ 30-70%) the community. On the other hand, in both mycosphere samples (MY and MR), Yersinia (22-62%), next to Pseudomonas (3-26%) turned out to be a dominant type (Figure 8.2C and Figure 8.2D). However, Burkholderia/Paraburkholderia, although still present, decreased to 1% or below. We conclude that, next to Burkholderia and Paraburkholderia communities, other fast-growing bacteria are able to compete utilizing nutrients in LB medium. Moreover, principal coordinate analysis based on Bray-Curtis distance showed that bacterial communities – before induction - clustered by sample type, i.e. B, MY and MR. However, this clustering dwindled away upon inductions (Figure 8. 2D). Additionally, we also found that sample type significantly influenced the α-diversity (Two-way ANOVA, F = 0.0362, p< 0.05) values (Figure 8.2E), which confirmed our hypothesis that diversity in mycosphere and bulk soil communities were different. However, again the significance of this difference disappeared after MMC induction (Figure 8. 2E).

To examine the efficiency of MMC induction and to detect the morphology of virus-like particles (VLPs) in all cultures (incuded and un-induced) (Figure 8.2F), we used transmission electron microscopy (TEM). An extensive screening found no observable VLPs in the bulk soil samples (induced and un-induced). Moreover, two dominant morphological types, i.e. phages belonging to the Siphoviridae and Podoviridae, were observed in the MMC-induced samples from both mycospheres. This result suggests that the virome, indicated by the observable VLPs, was raised upon MMC induction, while the virome from bulk soil communities were below the detection limit of TEM.

8

Mycosphere virome structure

To characterize the mycosphere viromes, we performed shotgun sequencing of all (control and MMC-induced) PCAT cultures. Across the samples, we obtained 450,342-1,6743,484 reads (Supplementary Table 8. 1) that were assembled into an average of 32,219.14 contigs across the dataset. Of all contigs, (10,937) 2,2%, (5,088) 6,5% and (16,494) 10% were ≥ 1Kb (contigs from MY, MR and B, respectively). These were screened (using VirSorter) for the presence of viral signals/regions (see Materials and Methods) (Supplementary Figure 8. 1B). The results showed that 989 contigs (from total of 32,519 >1Kb contigs) had viral signals (VirSorter category 1, 2, 4 and 5). Next, to determine whether the current virome contigs revealed similarity to any Paraburkholderia prophages, we clustered the previously-found 26 viral contigs from Paraburkholderia genomes (Pratama et al., 2018) to the current 989 virome contig data. This resulted in 521 “viral clusters-VCs” (Supplementary Table 8. 2). After manual inspection, a final set of 201 VCs, which encompass contigs ≥10Kb and/or circular (assumed to be complete) ones, was produced. This 201 VC set (spread among the three sample types) (see Figure 8. 3B) will be the basis of the further analyses in this study. Interestingly, we found that phage ϕ437, that was previously isolated from P. terrae BS437, clustered with two almost identical mycosphere virome contigs (from both MY and MR). This clustering has similarity of 99-100%, suggesting that ϕ437 is exclusively of mycosphere origin. From the 201 VCs, about 55, 34 and 22 VCs were unique for MY, MR and B, respectively. We also found about 50 VCs to be shared between the MY and MR viromes. Further, 10 VCs were shared in all samples (MY, MR and B).

To determine the protein sharing among the 201 VCs, orthologous protein cluster

analyses were done using OrthoVenn (cut-off 10-5_{), with a total of 12,441 amino acid}

sequences (from 567 viral contigs in 201 VCs). The results showed that MY and MR shared the highest number of orthologous proteins, 55-338 shared proteins, between the control and MMC induction treatments, respectively. On the other hand, the two mycosphere samples shared low numbers of orthologous genes with the B samples, i.e. 22-45 proteins. This suggests that the two mycosphere viromes tended to have a higher similarity as evidenced by shared genes, compare to the bulk soil virome (Figure 8. 3A and Figure 8. 3B).

A B C B-control MR-control MY-control B-MMC MR-MMC MY-MMC 114 77 26 61 178 11 39 18 338 28 55 23 66 4 12 16 3 4 33 133 165 115 51 14 259 14 2 10 27 45 39 112 4 1 3 24 14 4 2 ₂₀ 6 9 2 1 2 0 16 6 14 3 7 3 1 1 4 1 0 0 0 1 11 22 48 50000 150000 250000 34 50 55 10 0 3 0 2 0 0 0 3 22 0 22 B MR MY Para-PP

Figure 8.3. Composition of the virome dataset. (A) X axis: Sequence size; Y axis: cumulative

coverage, (B) Shared viral clusters among samples and (C) Shared proteins among samples.

Shared protein analyses were done using OrthoVenn (cut-off 10-5_).

Use of tRNA sequences as proxies for host prediction

The tRNA is often the integration site of viruses in bacterial genomes (Williams, 2002). In a previous study in our lab, Paraburkholderia terrae strain BS437 prophage ϕ437 was shown to contain integrated host tRNA (Pratama and van Elsas, 2017). We investigated the possible hosts of the viruses detected as contigs, using tRNA signatures. Using a stringent cut-off (see Materials and Methods), we identified 131 tRNAs in 71 VCs (35.32% of all VCs). However, only 101 tRNAs could be assigned to genus level (Supplementary Table 8. 3). Phage ϕ437 like sequences, which were among the final 201- VC data set, indeed showed the presence of tRNA sequences similar to those of the host P. terrae BS437, as reported previously (Pratama and van Elsas, 2017). Taking this analysis as the basis, the dominant hosts of the different mycosphere viruses in our current data set clearly belonged to the Proteobacteria. Among these, we found just two hits to Paraburkholderia (VC14 and VC20) and two to Burkholderia (VC179 and VC183), in the midst of numerous hits to hosts belonging to the γ-Proteobacteria (Supplementary Table 8.2). Moreover, hits to Bacteroides, Firmicutes, Chloroflexi and Actinobacteria were also found (Figure

8

8.4A). Interestingly, the composition of potential hosts matched what was found by

the 16S rRNA analysis (Figure 8. 2B).

Achromobacter (CP017433) P . oxyphila (AB488693) P. phytofirm_ansPsJN (AY49747) P. ca riben sis D SM 13 236 T P. terrae BS437 P. gramin is (U9693 9.1) P. _sp .M F2_-27 P. b_ann ens_is (AB 561 874) P. terrae NBRC100964 P. c arib iens is M WAP 64 P. phytofirmans BS455 Burkh olderia (CP0 2349 9) P. terrae DSM 1 7804T P. n_odo sa (AY₇₇ 31₈₉ ) P . ph ym at um (A J3 02 31 2) P. glathei (U96935.1) P. _he leia (A B4₉₅ 12₃₎ P. x_eno vora_ns (U 863 73.1) 91 50 77 88 86 99 86 62 92 99 95 86 95 0.01 Paraburkholderia Burkholderia Achromobacter

B

Proteobacteria

A

β γ 5 α Actinobacteria Firmicutes Bacteroidetes 6 3

The number of metagenomic viral contigs assigned to each phylum

3 1 0.1 80 Chloroflexi 1

Figure 8. 4. (A) All predicted hosts based on tRNA analysis. (B) β-Proteobacteria dominated by

Paraburkholderia spp. and Burkholderia spp. hosts in virome data (blue circles).

How similar are the mycosphere VCs potentially infecting Paraburkholderia and Burkholderia spp. (VC14, VC20, VC179 and VC183) to the previously isolated prophage ϕ437 from P. terrae BS437? As mentioned above, viral cluster VC20 showed >99-100% identity with phage ϕ437. This suggests that Paraburkholderia terrae BS437-like organisms were growing out from the mycosphere samples, and, upon induction, released ϕ437-like phages. Interestingly, the other Paraburkholderia and Burkholderia associated viral clusters (VC14, VC179 and VC183) were singletons and originated from the B soil virome. BLAST-P analysis showed that VC14 is relatively close to Pseudomonas aeuriginosa phage H66 - AGC34680 (92% coverage, 32% identity). Phylogenetic analysis based on integrase proteins, however, indicated that VC14 was close to B. cenocepacia phage BcepIL02 (FJ937737), next to VC20/ϕ437. Moreover, VC179 and VC183 were closely related to one another (Figure 8. 4A); BLAST-P analysis showed best hits to B. cenocepacia phage vB BceM AP3 - KP966108 (95% coverage and 85 identity). This confirmed the comparative analysis, which showed that VC179 and VC189 had high synteny in phage lysis and structural genes (Figure 8. 4B).

The ecological significance of the two district phages that were induced from - and potentially infect - Paraburkholderia spp. in mycosphere and bulk soil communities could lie in the genomic potential offered to their hosts. In the case of VC20/ ϕ437, we previously surmised that moron gene amrZ may be involved in enhancing biofilm formation (Pratama and van Elsas, 2017); biofilm formation is a strategy of P. terrae to interact with soil fungi (Warmink and Van Elsas, 2009). Moreover, the Paraburkholderia bulk soil phage VC14 also carried some potential moron genes that could potentially enhance host fitness. One moron gene was identified as a gene encoding a “granule-associated protein”. BLAST-P analyses of this predicted protein showed 94% identity (94% coverage) to a poly hydroxyalcanoate granule associated protein (phasin) of Burkholderia sp. GAS332 (WP_074290515.1) and to a similar protein of Paraburkholderia xenovorans (92% coverage, 59% identity). In Pseudomonas, polyhydroxyalkanoate-granule-associated protein (PhaF) promotes bacterial growth in nutrient (nitrogen) -limited conditions (Prieto et al., 1999). This finding suggests that VC14 phage might help the growth of Paraburkoholderia spp. under limiting (nitrogen) in bulk soil; forest bulk soil indeed is often poor in nitrogen versu rich in carbon sources.

Furthermore, the high numbers of genes for hypothetical proteins found in the current virome data set (30%) (Supplementary Table 8.4) showed a potentially large genetic potential offered by the virome to the host. However, studies of these proteins still need to be done to fully understand the extent of the role of the here studied virome in the mycosphere system.

Discussion

There is a huge knowledge gap with respect to the ecological role of viruses, especially in soil ecosystems. However, recent viral metagenomics (virome) studies have found that marine viruses have a significant role in ocean’s nutrient cycle (Roux et al., 2016). One would expect that soil viruses are equally important to soil functioning and the

soil microbiome, considering the raised number of soil viruses, of about 109 _{per gram}

soil (Williamson et al., 2013) and the high diversity of soil bacteria (Breitbart, 2012). In this study, we focused on the bacterial dynamic and the diversity of its virome upon induction.

Our analyses showed that the Russula mycosphere soil samples have higher bacterial abundance and diversity compared to bulk soil samples. This was shown by the CFU counts, the 16S rRNA gene counts and the PCR-DGGE analysis (Figure 8. 1). To focus on the Paraburkholderia communities, we cultured the mycosphere samples on PCAT medium, which has previously been used to isolate (Para)burkholderia

8

from soil (Salles et al., 2002). The PCAT medium was found to effectively isolate the

target organisms, , but was not fully specific, as it also allowed growth of the gamma-Proteobacterium Pseudomonas (Figure 8. 2D).

To investigate the prophage structures of these culturable communities, we used MMC inductions. This was based on the notion that, in more bacterial-dense environments such as the mycosphere, soil phages tend to favor lysogeny over lysis. This is consistent with the fact that resistance to superinfection in lysogens becomes increasingly important in crowded situations and the energy costs of carrying resistance to infection should be less than through mutation (Knowles et al., 2016). Upon induction, we observed by TEM that Siphoviridae and Podoviridae dominated in the MMC-induced lysates from mycosphere samples, versus an absence of phage particle detection in bulk soil samples (Figure 8.3B). This was consistent with the fact that we identified more viral signals in mycosphere samples as compared to bulk soil ones (Figure 8.4A). Interestingly, among these identified viral contigs, we found Paraburkholderia terrae strain BS437 prophage ϕ437 to be prevalent in the communities grown and induced from the mycospheres of both soil fungi. This is in line with a previous study, where P. terrae was isolated from mycosphere samples (Warmink and van Elsas, 2008) and ϕ437 excised spontaneously and upon MMC induction (Pratama and van Elsas, 2017). Moreover, the positive detection of this phage/host combination across mycospheres of Lyophyllum and Russula spp appears to indicate that the combination is ‘broadly’ generalistic with respect to mycosphere niche space.

Interestingly, based on the virus-host connection analysis using tRNA similarity, we also found one viral contig infecting Paraburkholderia spp. (VC14) in bulk soil samples (Supplementary Table 8. 2) and two other viral contigs infecting Burkholderia (VC179 and 183). The BLASTP and phylogeny analysis based on the phage integrase showed that VC14 was related to Pseudomonas aeruginosa virus H66 and other Burkholderia viruses, i.e. B. cenocepacia phage BcepIL02 (FJ937737) and P. terrae phage VC20/ϕ437, respectively. On the other hand, VC179 and VC183 were closely related to Burkholderia phage vB BceM AP3 - AKA61172. The ecological role of the two district Paraburkholderia viruses might be linked to the ecological niche of their host. The mycosphere viruses ϕ437/VC20 possibly enhance biofilm formation as a strategy of P. terrae to interact with host fungi (Warmink and Van Elsas, 2009). In contrast, the VC14 viral type infecting bulk soil Paraburkholderia, as it was carrying poly hydroxyalkanoate-granule-associated protein (PhaF), possibly advances the fitness of its host, allowing growth under nutrient (nitrogen)-limited conditions in the forest floor.

Acknowledgements

We thank Maryam Chaib De Mares for helpful discussion.

Authors’ contributions

Collected the samples, isolated the bacteria, extracted metagenomics DNA: AAP. Performed MMC induction: JE. Viral precipitation, DNA extraction and sequencing: AAP. Viral metagenomics analysis: AAP. Bacterial amplicons 16S rRNA analysis: XJ. Wrote the manuscript: AAP and JDvE. All authors read and approved of the manuscript.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by a scholarship of the Indonesia Endowment Fund for Education (Lembaga Pengelola Dana Pendidikan, - LPDP, Departemen Keuangan, Republik Indonesia) to AAP. The China Scholarship Council (CSC) to XJ.

8

SUPPLEMENTARY INFORMATION CHAPTER 8

A novel inducible prophage from the mycosphere

inhabitant Paraburkholderia terrae BS437

The dynamics and diversity of the mycosphere virome

Akbar Adjie Pratama, Xiu Jia, Juliana Eschholz and Jan Dirk van Elsas In preparation for publication

Supplementary_files_chapter8 (shared link below)

Supplementary tables Chapter 8: Table S8.1. Virome metagenomics assemble statistics; Table S8.2. Viral clusters; Table S8.3. Possible hosts based on tRNA; Table S8.4. Annotation of 201 VCs.

Please download here: https://drive.google.com/drive/folders/1bHnLhS-GerhwWC EiFIK6g2unwFS7e8Gz?usp=sharing

> 1 Kb _{Dereplicate (95% Identity} , 80% coverage), > 10 Kb Taxonomy Soil ecology V

irome quality check

Trimmed Assembly Virsorter 989 contigs

+

Paraburkholderia

spp. prophages

26 contigs _Viral clusters

(201

VCs- 567 contigs)

B

Russula ochr

oleuca

Russula emetica Bulk soil

Mycosphere and bulk samples

Plate culture (R2A, PCA

T) DGGE CFU count 16S rRNA qPCR PCA T culture MMC induction TEM VLP precipitaion DNase treatment V iral DNA extraction Total DNA extraction initial 8h 24h (c on tr ol and MMC ) Library preparation

Next generation sequencing

A

Supplementary figure 8.1

. Overview of the research scheme. (

A) the sampling plans and (