Metabarcoding reveals southern hemisphere fungal endophytes within wood of cultivated Proteaceae in Portugal

Endophytic fungal hitch-hikers have been difficult to detect in the past, and have potentially spread these latent pathogens via the global plant trade. The African genera Protea, Leucospermum and Leucadendron, commercially referred to as proteas, form the basis of a global flower production industry. The largest producers of proteas are Australia and South Africa, followed by Portugal and Spain. In the 1990s propagation material from South Africa was used to establish protea orchards in Portugal. We utilized metabarcoding to determine if this plant trade has carried host-specific fungal pathogens to a new environment. Wood samples collected from asymtomatic twigs from Portuguese farms, where propagation material had been imported from South Africa, was compared to material from South African farms that originally produced and supplied rooted and unrooted cuttings. DNA metabarcoding, using fungal-specific primers for the ITS2 gene region, produced 1237 OTUs. Focusing only on known pathogens of protea, we found that the Portuguese orchards contained fungal disease agents associated with Proteaceae or other plant families from the Southern Hemisphere. Our sampling technique could be used by agencies and applied to other plant material and pathogens to reduce the spread of pathogens.


Introduction
Trade of plant material is closely linked to the dispersal of plant pathogens Liebhold et al. 2012;Santini et al. 2013). There are numerous examples of fungal "hitchhikers" in the global trade of plant material (Santini et al. 2013). This can be a result of asymptomatic colonization of hosts and the abilty of the fungi to remain quiescent within the host tissue for long periods of time Crone et al. 2013;Hernandez-Escribano et al. 2018). This is a common phenomenon for species in the Botryosphaeriaceae; for example, the global movement of Diplodia sapinea was most likely caused by the import of infected pine seeds for the Pinus radiata plantation industry (Burgess and Wingfield 2002).
While the mechanisms of pest and pathogen introduction via plants has been addressed in the United States (Liebhold et al. 2012) and Europe Santini et al. 2013), little is known regarding the origin and main dispersal pathways of invasive fungi, especially endophytes. Such information is needed in order to estimate the relative risks related to the accidental introduction of these fungi. There have been some investigations considering endophytes (Sakalidis et al. 2013) and other cryptic plant pests conducted at high risk sites such as harbours, airports, and large nurseries in importing countries (Migliorini et al. 2015;Rassati et al. 2015). Yet there is a relative paucity of studies that address the issue of hitch-hikers in the plant trade.
South Africa has one of the greatest genetic resources of Proteaceae in the world (363 species), second only to Western Australia (697 species). The African genera Protea, Leucospermum and Leucadendron, commercially referred to as proteas, are the most commonly cultivated taxa in the family. Several cultivars of species in these genera have been developed and extensively propagated during the past 50 years. Germplasm of these plants has been used to develop the flower production industry in many countries, including Australia and Portugal (Crous et al. 2013).
DNA-based tools can be applied to detect fungal endophytes in specific hosts. They make it possible, to determine if the presence of these fungi could result from the importation of plant propagation material. This has been discussed for plant seeds (Cleary et al. 2019;Franić et al. 2019;Jimu et al. 2016) and ornamentals (Vettraino et al. 2015(Vettraino et al. , 2017.
Our investigation used DNA metabarcoding to detect fungi associated with asymptomatic twigs collected to determine the possible role of Proteaceae in the accidental introduction of non-native endophytes into Portugal. Samples were collected from Portuguese farms, where propagation plants had been imported from South Africa 20 years previously. These were compared with samples from two South African farms, one of which had dispatched rooted and unrooted cuttings to the Portguese nurseries. The samples collected were asymptomatic and we focused on the fungal endophytes within these samples.

Sites sampled
Four ornamental protea orchards were selected, two in southwestern Portugal, and two in South Africa, respectively in the Gauteng and the Western Cape provinces. The two South African sampling locations had very different climates; a Mediterranean-like climate with wet winter and dry summers in the Western Cape, and wet summers and dry winters in the Gauteng savannas. The two Portuguese orchards were in close proximity, and had the same climatic profile. The Portuguese orchards had been established approximately 20 years previously, using stem and rooted cuttings provided by the farm in the Western Cape of South Africa. Two sampling plots per orchard consisted of a 50-60 m transect with five sampling units equidistant along the transect, consisting of five plants each. The number of plant sampled was 50 per orchard for a total of 200 plants in the entire study. Common to all four orchards, the cultivars 'Pink Ice' (Protea compacta x Protea susannae P Matthews) and ' Safari Sunset' (Leucadendron salignum x Leucadendron laureolum; Stevens and Bell, 1962) were examined.

Processing of samples
Five twigs per sampling unit were collected. After the leaves were removed, samples (1 cm) were cut from the distal ends and surface disinfested for 1 min in 2% sodium hypochlorite, for 1 min in 70% ethanol, with washes in sterile distilled water for 1 min after each of the treatments. To collect material for DNA extraction, sterlie secateurs were used to excise a 1 cm portion from the distal end of each sterile stem sample after which the bark was removed from these pieces and discarded. Shavings were taken from the remaining section of wood using a sterile scapel. Shavings from the same sampling unit were bulked and allowed to dry naturally in paper bags for 7 d before being transferred, under sterile conditions, to 2 mL microcentrifuge tubes (Pirouet) and stored at −80°C.

DNA extraction, library preparation and sequencing
Approximately 200 mg of wood shavings were ground with TissueLyser (Qiagen-manufactured by Retch) for 4 min using two sterile 0.5 mm dia. Iron balls. Approximately 50 μg of this ground wood was transferred into the extraction kit tubes, and total genomic DNA was extracted with the isolation kit (Mo Bio Laboratories, Power Plant® Pro DNA Isolation Kit Cat# 13400-50) following the manufacturer's protocol. Extracted DNA was amplified using the fungal specific primer fITS7 (Ihrmark et al. 2012) and eukaryotic primer ITS4 (White et al. 1990) as the reverse primer which amplified the ITS2 gene region the universal genetic barcode for fungi (Ihrmark et al. 2012). PCR reaction mixtures (25 μl) contained 2 μl of genomic DNA, 1 μl of each 10 μM primer, 12.5 μl of KAPA Taq DNA Polymerase (Roche Applied Biosystems, Nutley, NJ, USA). Amplifications were carried out using the following PCR conditions: 5 min at 94°C; 32 cycles of (30 s at 94°C; 30 s at 57°C; 30 s at 72°C); 7 min at 72°C. Two extraction controls and two PCR controls were included in the analysis.
Libraries for the MiSeq sequencing platform were prepared and indexed with the Nextera®XT Index kit following the manufacturer's instructions (Illumina Demonstrated Protocol: 16S Metagenomic Sequencing Library Preparation). To estimate the amount of amplified DNA, PCR products were separated by electrophoresis on gels containing 1% (w/v) agarose (Molecular Biology Grade; Fisher Biotec). Band lengths were estimated using an Axygen 100-30,000 bp ladder DNA Marker (Fisher Biotec). All extraction and PCR controls were negative. DNA quality was re-evaluated by electrophoresis on 1% (w/v) agarose LE gels (Genespin) stained with GelStain (Invitrogen-Sybr®Safe DNA gel stain) and quantified with Qubit® 2.0 Flurometer (Invitrogen). Forty seven uniquely indexed PCRs were pooled to create the library for the sequencing run, which was performed on an Illumina MiSeq Reagent Nano Kit V2 (250 reads either direction) following the manufacturer's recommendations.

Data analysis
Raw Illumina MiSeq reads were analysed using the QIIME 1.9.1 pipeline (Caporaso et al. 2010). After demultiplexing, paired end reads were joined, and reads <200 bp with a Phred score < 20 were filtered out. Removal of chimeric sequences was done using Usearch 6.1 (Edgar 2010;Edgar et al. 2011) and the UNITE database (Kõljalg et al. 2013). Operational Taxonomic Units (OTUs) were selected using Uclust v1.2.22q (Edgar 2010) and clustered with 99% sequence similarity. The cluster seed (the first sequence of the cluster) was picked as the representative sequence per OTU cluster. Identification of operational taxonomic units (OTUs) was conducted through BLAST searches in UNITE. An OTU table was created and singletons, no blast-, Plant-and Protista-hits were removed. Representative sequences were obtained using Geneious version 10.0.9 (Kearse et al. 2012). Taxonomic identity of representative sequences was checked using the NCBI nucleotide BLAST website ("https://blast.ncbi.nlm.nih.gov/Blast.cgi?P") and the best hit was exported. In order to be sure that the OTU being considered was a known associate of the Proteaceae, or a known Southern Hemisphere plant pathogen, a species identity was inquired (Table 1). Thus, all sequences with an identity value lower than 97% of the best hit, were removed. Taxonomic identities were assigned to OTU's, firstly by blast searches in Genbank, and secondly by phylogenetic comparison to related species previously published in scientific articles. The data were then filtered further to include only canopy-associated fungi endemic to the Southern Hemisphere and fungi known to be associated with Proteaceae. Using the representative OTU (n = 144), identities were assigned to a species (16 species in total, ie more than one OTU corresponded to the same species) with phylogenetic analysis (Figs. S1-S1, 4 Table S1). The number of sequence reads per OTU was obtained from the filtered OTU table.

Statistical analyses
Statistical analyses were performed to determine if provenance (country, environment) and plant variety influenced the communities of fungi known to be associated with proteas. The presence or absence of the selected fungal species was used in a permutational multivariate analysis of variance (permanova) performed with the 'adonis' function in the vegan R package (Oksanen et al. 2018). The Jaccard index was used to calculate a dissimilarity matrix and the analysis was run with 9999 permutations. Post-hoc comparisons with a Bonferroni P value adjustment were performed using the pairwise Adonis R package (Martinez Arbizu 2017). The 'betadisper' and 'permutest' functions in the vegan R Ramularia spp. complex consists of a group of sixty-six reads with an equal similarity of 99% with GenBank deposited sequences of species R. proteae and R. stellenboschensis. R. proteae was isolated on Protea longifolia in Tasmania (Crous et al. 2000), R. stellenboschensis was isolated in Western Cape on Protea sp. (Crous et al. 2011b). In the text they are referred to as "Ramularia complex" and considered as a unique taxonomical group. This did not affect our results because both two suggested taxa have been exclusively isolated in Southern Hemisphere countries c Rachicladosporium eucalypti, within the fungal species selected in this work, is the only one not associated to Proteaceae that was described in the Northern hemisphere (Addis Abeba Botanical Garden). Due to the fact that this species is known from only a limited geographic area its presence in the sampling location was considered of a particular interest and it was included in the results package tested the assumption of homogeneity of multivariate dispersions between factor levels. Results were considered statistically significant at P values <0.05. The fungal communities were represented using unconstrained ordination, non-metric multidimensional scaling (NMDS) from the 'metaMDS' function in the vegan R package. The NMDS was generated with the Jaccard index and two dimensions. The results of the NMDS were extracted, and plotted using the GGplot2 R package (Wickham 2016). Ellipses were added to plots with the default 'stat_ellipse' function, and all points were randomly jittered (± 0.1 vertically and horizontally) to display the abundance of samples plotted at the same coordinates. All analyses were performed in R version 3.5.1 (R Core Team 2018).

Identification of fungal species
Data from the Portuguese farms were merged due to their proximity, while the data from Gauteng and Western Cape farms were kept separate (Table 1). Sixteen species (16,651 reads) of known association with Proteaceae, or representing known pathogens of Southern Hemisphere woody plants, were identified (Table 1). Phylogenetic assessment of the OTU's taxonomy is provided in Fig. S1-S14 and Table S1. Western Cape and Portugal orchards had similar numbers of species identified (11 and 10, respectively) while only five species were identified in the Gauteng orchard samples.
The highest similarity coefficient of the selected fungal taxa occurred between the two South African farms (31.3%), while Portugal and Western Cape farms had a higher similarity (23.8%) than Portugal and Gauteng farms (13.3%).
Five species were common to both the Portuguese and South African orchards. These included Aureobasidium leucospermi, Catenulostroma protearum, Neodevriesia fraserae, Leptosphaeria proteicola and Ramularia spp. complex (Table 1). Four of these species have been reported in Western Cape in the past 20 years. Five species (Batcheloromyces leucadendri, Exophiala capensis, Neocladosporium leucadendri, Neofusicoccum protearum and Neodevriesia knoxdaviesii) were exclusively detected in the Portugese samples (Table 1). Of these five species, four have been reported in Western Cape within the last 12 years.

Community analyses
NMDS2 analysis showed that the strongest driver of variation was between the fungal communities in the two counties (Fig. 1a). While there was some overlap in fungal communities associated with proteas between the two countries (Fig. 1a), there remained a significant difference (p = 0.0001; Table 2). A significant interaction (p = 0.0428) between protea variety and country was driven by the differences between the fungal communities associated with Portugal and South Africa (Table 2; Fig. 1b). Within each country, the fungal communities associated with the varieties were not significantly different (p = 0.1428) and the covariate accounted for little variation. Although all environments were significantly different to each other (p = 0.0001, Table 2; Fig. 1c), there was greater overlap between Gauteng and Portugal, than between the Western Cape and Portugal.

Discussion
Known fungal endophytes, including putative or known pathogens of South African Proteaceae, were detected in samples from Portugese protea orchards. This supports the primary hypothesis of this study that these pathogens could have been transported with propagation material between continents. The most notable of these is Neofusicoccum protearum, a well recognized canker pathogen of proteas in Southern Africa, Australia and Northern Hemisphere countries where proteas are cultivated (Crous et al. 2013). As previously reported by Denman et al. (2003), N. protearum is host-specific and its exclusive association with South African Proteaceae (subfamily Proteoideae) suggests that it is native and common in South Africa, even if it was not found in the South African samples in the present study.
Neofusicoccum protearum is sporadically reported in Australia, whereas other Neofusicoccum species are more commonly found associated with proteas (Burgess Tan et al. 2019). Hence, it can be assumed that it has most likely been introduced in other countries through the trade of South African protea plants. N. protearum causes lesions in cultivated Proteaceae mostly localized in the cambial tissues; necrosis often begins in the leaves passing to the twigs through the petioles and eventually killing branches (Crous et al. 2013). Since N. protearum is a canker pathogen that enters the woody stems through the leaves, we may have Fig. 1 Non-metric multidimensional scaling ordination of the fungal communities associated with Proteas grown in Portuguese and South African orchards. Points and ellipses are categorised by a) country, b) variety of Protea, country and c) environment. Small clusters of points represent the same fungal community as positions were slightly jittered to display the abundance of samples plotted with the same NMDS coordinates missed this pathogen when only examining wood of asymptomatic young twigs, tissues not commonly colonized by this fungus.
Fungi found in Portugese samples included eight species known to be pathogens of Proteaceae in South Africa. These were B. leucadendri, N. leucadendri, N. protearum, N. knoxdaviesii, A. leucospermi, C. protearum, L. proteicola and R. proteae. Of these, only N. protearum and A. leucospermi have previously been reported from the Northern Hemisphere (Crous et al. 2008(Crous et al. , 2013Crous et al. 2011a;Crous et al. 2011b). Previously, C. protearum has been collected from Australian Proteaceae in South Africa (Crous et al. 2009), while N. leucadendri and N. fraserae have been found in Australia (Bezerra et al. 2017;Crous et al. 2010). It is highly unlikely that these fungal species preexisted in Portugal prior to the introduction of proteas and thus they have probably been introduced with plants imported from South Africa.
Four species of fungi detected in Portugese samples have previously been reported from woody plants in familes other than the Proteaceae in South Africa (B. leucadendri and E. capensis) and Australia (N. leucadendri and N. fraserae). These findings suggest that they may have originally been transported to Portugal with infected plant material from the Southern Hemisphere. Because some of the known hosts are not Proteaceae, those fungal species were possibly transported on a host species belonging to a different plant family. They could then have become naturalized in Portugal, and subsequently jumped to Protea species, as new hosts. This would be consistent with the numerous examples of such host jumps (Slippers et al. 2005) between congeneric tree species (Gladieux et al. 2015) or those at a family level such as in the case of Austropuccinia psidii from native Myrtaceae to Eucalyptus species in South America (Coutinho et al. 1998) or between families such as found for the canker agent Chrysoporthe cubensis in Myrtaceae, Melostomataceae and Lythraceae (Gryzenhout et al. 2004(Gryzenhout et al. , 2006Rodas et al. 2005). Alternatively, although these four fungal species have not been reported as pathogens of South African Proteaceae, they could have been present as asymptomatic endophytes in the material sent to Portugal.
Metabarcoding has recently been used to illustrate how known pathogens are distributed as endophytes in woody plants like chestnut trees (Castanea sativa) (Fernandez-Conradi et al. 2019) and grape vines (Vitis vinifera) (Dissanayake et al. 2018). From an ecological prospective, it has been used to describe the taxonomic variation of plant-associated cryptic fungi linked to changes in environmental variables, for example the effect of altitude and leaf biochemistry on the diversity and composition of the leaf mycobiome in beech (Fagus sylvatica) (Unterseher et al. 2016). Other examples include the effect of the elevation gradient in the variation of phyllosphere fungi in the Chinese shrub Mussaenda shikokiana (Qian et al. 2018), the effect of temperate forest tree species richness and phylogenetic diversity in the variation of foliar fungal diversity (Griffin et al. 2019) and the influence of different forest tree stands on soil fungal composition (Jimu et al. 2017). Consistent with these findings, the present study contributes to the validation of the use of metabarcoding techniques, introducing its application to study of invasive plantassociated fungi. In particular, the possibility to confirm the presence of certain fungal species, demonstrates that protea orchards in Portugal harbor many potential pathogens, the majority of which are known only in the Southern Hemisphere. With a degree of resolution unattainable using any other technique, metabarcoding allows the detection of fungal endophytes present in extremely low quantities inside the collected specimens.
All fungal species in the current study were detected from asymptomatic twigs, and all (with the exception of N. protearum which is also a canker pathogen) have previously been described as leaf-associated fungi. As efforts to detect harmful microorganisms are likely to be based on ureliable visual inspections of plants (Brasier 2008), latent disease agents present as endophytes in plant propagation material can easily enter a new country undetected, where they may potentially establish and spread (Slippers and Wingfield 2007). With further development, metabarcoding might enhance the ability of biosecurity to detect particular fungal pathogens Paap et al. 2017;Wingfield et al. 2015). Existing technologies of DNA metabarcoding, and an improved access to, and support of, metadata could complement phytosanitary inspections to assist in reducing particular plant pathogens linked to the global trade of plants. This study provides an applicable technique to allow detection of a greater number of fungal hitckhikers than can be achieved using classic techniques. For the present, we suggest that rather than importing new propagation material from other continents, European protea growers should produce the required plant material locally. The use of tissue culture plants for export purposes would not completely eliminate fungal endophytes (Barrow et al. 2004;Pirttila et al. 2002Pirttila et al. , 2003, but could reduce the risk of spreading pathogens between continents.
Acknowledgments We are grateful for the assistance provided by the European Cooperation in Science and Technology (Action FP1401: A global network of nurseries as early warning system against alien tree pests -Global Warning). We also thank Kay Howard for assistance with the manuscript preparation and Chris Shaw for his work on the statistics.
Authors' contributions DM, TB and MW have made major contributions to the conception of the study and the writing of the manuscript, DM collected all the samples. DM and MM were responsible for the acquisition, analysis, and interpretation of the data. AS, AR and PT facilitated sampling and were involved in the critical editing for important intellectual content.
Data availability The data that support the findings of this study are openly available in NCBI at https://www.ncbi.nlm.nih. gov/sra/PRJNA548865, SRA reference number PRJNA548865. Temporary Submission ID: SUB5845073.

Declarations
Ethics approval (Not applicable)

Consent for publication (Not applicable)
Conflicts of interest/competing interests The authors declare no conflicts of financial or personal interest.