**Does Arbuscular Mycorrhiza Determine Soil Microbial Functionality in Nutrient-Limited Mediterranean Arid Ecosystems?**

**Neji Mahmoudi 1, Teresa Dias <sup>2</sup> , Mosbah Mahdhi <sup>3</sup> , Cristina Cruz 2, Mohamed Mars <sup>1</sup> and Maria F. Caeiro 4,\***


Received: 5 May 2020; Accepted: 4 June 2020; Published: 10 June 2020

**Abstract:** Arbuscular mycorrhizal fungi (AMF) are determinant for the performance of plant communities and for the functionality of terrestrial ecosystems. In natural ecosystems, grazing can have a major impact on mycorrhizal fungi and consequently on plant growth. The objective of this study was to evaluate the statements referred above in Mediterranean arid areas in Tunisia. Root samples and rhizosphere soils of five dominant herbaceous plants were studied at six distinct arid sites differing on soil proprieties and grazing intensity. At each site, chemical and dynamic properties of the soil were characterized as well as the AMF colonization intensity and the soil functionality. Results showed that the mycorrhizal frequency and intensity and spore density, varied between plants in the same site and, for each plant, between sites and evidenced a positive effect of mycorrhized plants on soil microbial activity. Grazing and soil properties strongly affected AMF composition and the soil microbial and biochemical dynamics, which presented the lowest values at the sites with the highest grazing intensities. In conclusion, these results demonstrate that AMF improve soil biological properties, supporting the hypothesis that mycorrhiza and grazing compete for plant photosynthates, and highlight the importance of mycorrhizal symbiosis towards soil functionality under arid conditions.

**Keywords:** arbuscular mycorrhizal fungi; arid areas; biological properties; conserved areas; grazing; mycorrhiza

#### **1. Introduction**

Arid and semi-arid regions of the world are considered as being particularly vulnerable to climate change [1]. They are already climatically stressed with high temperatures, low rainfall and long dry seasons. These ecosystems are highly dynamic, with bursts of productivity in the wet season of some years, and very low productivity in dry years. The arid and semi-arid regions of the world have been subjected to accelerated desertification due to increasing grazing intensity, decreased rainfall, higher temperatures and prolonged periods of drought [2,3]. These pressures, associated with increased anthropogenic impacts and climate changes, caused the decline of forests, regression and extinction

of many pastoral and forage species, and accelerated soil degradation and change of soil microbial communities [1].

In most cases, the degradation process starts with the disruption of functional networks that protect and alleviate stress, conferring to the organisms involved protection to global changes [4], including land use changes. Plants are well known for their symbioses with nitrogen fixing bacteria and arbuscular mycorrhiza fungi (AMF). These are just two examples of a vast range of symbiotic relationships that consist the plant microbiome and modulate plant phenotype and, thus, plant fitness. In this context, the microbial community of the rhizosphere [5] is of great importance to plant performance, playing a crucial role in ecosystem functioning. The microbial activities of the rhizosphere also determine the bioavailability of nutrients and, therefore, soil fertility. Plants interact with guilds of these beneficial microorganisms living in their rhizosphere, promoting their growth and development [6].

Arbuscular mycorrhiza involves reciprocal complementary resource exchanges between plant roots and soil fungi and is widespread in natural ecosystems [7]. When colonizing the roots, AMF develop intra- and extra-radical mycelium. The hyphae of the intra-radical mycelium colonize the cells of the root cortex and penetrate the periplasmic space where they develop vesicles and arbuscules that are the structures responsible for most of the exchanges between the fungi and the plant [8,9]. The extra-radical mycelium spreads its hyphae beyond the root surface and colonizes the surrounding substratum, increasing the volume of soil explored, and creating a privileged space for microbial development. As part of the interaction, the host plant provides the fungus with carbon in exchange for nutrients and water is taken up by the fungus [8,10]. Thus, mycorrhiza plays a crucial role in terrestrial ecosystem functioning, especially in arid or semi-arid areas, where root exudates are the major carbon source supporting soil microbial activities [11]. Apart from the nutritional benefits, mycorrhization tend to increase plant tolerance to other stress (biotic and abiotic) conditions [12].

Despite the major ecological importance of AMF, there is a lack of knowledge on the effect of AMF on soil microbial properties and biological activities, or the effect of ecosystem management (including grazing) on mycorrhization. The complex herbivore–plant–AMF is dependent on the species involved in these interactions, the intensity and frequency of grazing, vegetation, and topography [13,14]. Some studies indicate that grazing influences soil biochemical activity, usually through the degradation of soil structure by trampling [15] and the imbalance of soil chemical, microbial and biochemical properties [16].

There is limited information about the effects of grazing on AMF communities and on soil microbial composition and function, especially for low fertility soils, as found in most Mediterranean arid and semi-arid regions. This work addressed both questions in three arid and semi-arid regions of Tunisia, by evaluation of grazing intensity on AMF–plant interactions, and the consequent impact on soil microbial communities.

#### **2. Materials and Methods**

#### *2.1. Study Site*

This study was developed at three conserved natural areas (Figure 1) in the arid and semi-arid ecosystem of Tunisia, under Mediterranean climate: Bou-Hedma National Park (Sidi Bouzid coordinates, 34◦39 N, 94◦8 E), Zarat-Gabes protected area (Gabes coordinates, 33◦41 N, 10◦23 E) and Oued Dkouk Natural Reserve (Tataouine coordinates, 32◦37 N, 10◦18 E). The annual temperature range is very high in these areas, with minimal and maximal monthly temperature means of 17 ◦C (January) and 36 ◦C (August), respectively. The mean annual rainfall in these ecosystems varies between 100 and 260 mm.

**Figure 1.** Location of the three conserved natural areas in the arid Mediterranean ecosystem of Tunisia.

Experiments were carried out in three natural protected areas of Tunisia (Bou-Hedma, Zarat, and Oued Dkouk) (Table 1). At each natural park samples were taken in two sites, one inside and the other outside the protected area. The studied sites differed in grazing intensity and soil type (Table 1). The three sites inside the conserved areas were subjected to a light grazing, while the other three sites outside were subjected to more intensive grazing by domestic herds of sheep, goats and camels.


**Table 1.** Details of the six sampling sites at the three conserved natural areas.

#### *2.2. Roots and Soil Sampling*

Five plant species (common to all study sites) and three plants per species per site were analyzed for AMF colonization. The herbaceous plants collected were: *Lotus creticus* (Fabaceae), *Medicago truncatula* (Fabaceae), *Astragalus corrugatus* (Fabaceae), *Malva aegyptiaca* (Malvaceae), *Diplotaxis simplex* (Cruciferacea).

Plant roots were carefully collected in order to access the fine active roots where mycorrhiza colonization occurs. Simultaneously, rhizosphere soil of each plant was also collected. For each site, a pooled soil sample (composed of 5 soil cores with 10 cm diameter and 20 cm length) was collected in an area without vegetation and used as a control (bulk soil) for the influence of the plant in the dynamic soil characteristics. Soils were sieved (2 mm) to remove the remains of plants, gravel and earthworms, and stored at 4 ◦C for further analysis. Soil samples were analyzed in triplicates.

#### *2.3. AMF Colonization Status*

According to the methods of Phillips and Hayman [17], AMF colonization was evaluated after observation of 30 root fragments per plant. Of each plant species, root segments (1–2 cm length) were submerged in 10% potassium hydroxide (KOH) at 90 ◦C for 45 min. After bleaching and

acidification steps with hydrogen peroxide (H2O2) and hydrogen chloride (HCl), respectively, root pieces were colored with 1% Trypan Blue solution. The duration of each step varied among plant species, according to the respective root diameter and surface root characteristics. For each plant species from each site, 30 root fragments were placed on slides and preserved with lactoglycerol. In all the stained roots, the presence of hyphae, vesicles, and arbuscules inside the root were viewed through a microscope (Nikon, Tokyo, Japan) at 400× magnification. The frequency (F%): number of colonized roots/total number of observed roots) and intensity of mycorrhization (M%): proportion, in percentage, of the root colonized by AMF) were calculated using the Mycocalc program (http://www2.dijon.inra.fr/mychintec/Mycocalcprg/download.html)

#### *2.4. Quantification of AMF Spore Density*

Using the wet sieving method from Gerdemann and Nicolson [18], AMF spores occurring in soil samples were extracted and quantified. Three nested sieves of 1000, 100, and 32 μm were used. For each soil sample, quantities of 100 g were submerged in 1 L of tap water. After a stirring step, the supernatant was sieved through the nested sieves. The spores that hold on to the two sieves of 100 and 32 μm were recovered in 5 mL centrifuge tubes. After a step of centrifugation on a viscosity gradient (sucrose solution at 60%), the supernatant retained was rinsed with distilled water to remove the sucrose solution. Retrieved AMF spores of each soil sample were counted under a stereomicroscope (40× magnification) and average numbers were calculated per 100 g of dry soil.

#### *2.5. Soil Analysis: Physical and Chemical Properties*

Soil pH and electrical conductivity were determined in a 1:1 (*v*/*w*) water: bulk soil suspension using a pH meter (Matest, Treviolo, Italy) and a conductivity meter (Bibby Scientific, Bibby Scientific, Staffordshire, UK), respectively [19]. Soil texture, which represents the granulometric distribution of its constituents (the proportion between small particles: clay, silt and sand), was calculated using the Robinson's pipette method [20]. The other main physicochemical soil characteristics (organic matter, total nitrogen, total phosphorus and water content) were determined by conventional analyses performed by the Soil Analysis Laboratory in the Regional Commissariat for Agricultural Development, in Gabes. The determination of the organic matter was carried out indirectly, starting from the determination of the organic carbon content of soil. The determination of total phosphorus was evaluated by a degradation acid reaction step followed by a dosing step carried out in an automated spectrophotometer (Shimadzu, Kyoto, Japan). The total nitrogen was determined following the Kjeldahl method.

#### *2.6. Microbiological and Biochemical Properties*

For the impact of AMF on microbiological parameters, the carbon of the microbial biomass (Cmic) present in the plant rhizosphere and bulk soil were evaluated following the "fumigation–extraction" technique [21]. This method is based on three essential steps: a fumigation step with chloroform, incubation in a 10-days fumigation period, and an extraction step with ninhydrin-N reactive and potassium chloride (KCl). This technique has been used to provide rapid and accurate measurements of soil biomass -C and -N.

To evaluate the effect of AMF on biochemical properties, alkaline phosphatase and β-glucosidase activities were calculated and evaluated according to the method of Caravaca et al. [22] in a spectrophotometer (Shimadzu, Japan) at 398 nm. The dehydrogenase activity was determined as described by Garcia et al. [23] in a spectrophotometer (Shimadzu, Japan) at 490 nm.

#### *2.7. Statistical Analyses*

Analyses of variance (ANOVA) for repeated measures using the XLSTAT (v2010.5.04) software (Addinsoft, New York, NY, USA) were ascertained to test the effect and the significant difference between the studied parameters. Least significant difference values at the 5% levels of significance (*p* ≤ 0.05) were calculated to assess differences between different values.

To evaluate the effects of the grazing parameter (explanatory variable) on mycorrhizal properties (F%, M% and number of spores) (response variables), ANOVA and Canonical Correlation Analysis (CCorA) were applied.

To evaluate the relationships between soil physical and chemical parameters and mycorrhizal properties (F%, M% and number of spores), a principal-component analysis (PCA) was applied.

Pearson's correlation was used to determine relationships between variables: mycorrhizal properties and microbiological and biochemical parameters.

To model the relationships between mycorrhizal properties and parameters of soil microbial and biochemical activities by linear regression, XLSTAT (v2010.5.04) software was used.

#### **3. Results**

#### *3.1. Physical and Chemical Properties of Soils*

The more resilient soil characteristics such as soil texture (Table 2) and pH did not show big discrepancies among sites inside and outside the protected areas. However, differences were obtained for the more responsive soil dynamic characteristics: electrical conductivity (E.c), total nitrogen (T.N), total phosphorus (T.P), organic matter (Org. Mat), and water content (Wat. Con). In general, the studied sites had an alkaline pH ranging from 8.0 to 8.4. The highest percentages of soil organic matter were observed in Site 1 (2.6%) as well as inside the other protected areas. Identical patterns were observed for the total nitrogen, the highest level (194 ppm) also being observed in Site 1, followed by Sites 3 and 5 (inside the protected areas). For the total phosphorus content, it was the contrary: the highest values were registered in the sites outside the protected areas (15.3 ppm was the highest value, registered in Site 6). Water content (varying between 1.3 to 3.3%) and electrical conductivity (ranging from 1.3 to 2.5 s·m<sup>−</sup>1) always presented the highest values in the sites inside each protected area.


**Table 2.** Physical and chemical properties of the six sampling sites from the three conserved ecosystems.

E.c: electrical conductivity; T.N: total nitrogen; T.P: total phosphorus; Org. Mat: organic matter; Wat. Con: water content. Letters a–f: significant differences (*p* < 0.05); mean and standard error values (*n* = 3).

#### *3.2. AMF Colonization of Plant Roots*

Three plants and 30 root fragments of each one of the five plant species were analyzed per site. Direct observation of the root samples under the microscope showed that the roots of all studied plants, except those of *Diplotaxis simplex*, were colonized by AMF. All the structures characteristic of root colonization by AMF (intracellular aseptate hyphae, vesicles and arbuscules) were observed. The highest mycorrhizal frequency (F%)) was observed for *M. truncatula* in all the studied sites. At each site, the mycorrhizal frequency (F%)) varied among the plant species (Figure 2). The plants with the highest AMF root colonization belong to the Fabaceae family (*M. truncatula*, *A. corrugatus* and *L. creticus*) in all the studied sites, always with higher values in the site inside each conserved area. There were significant differences (*p* < 0.001) in the mycorrhizal intensity (M%) among sites, again with highest M% for *M. truncatula* as well as, for each plant species, for the plants in the sites inside the conserved areas (Figure 2).

**Figure 2.** Mycorrhizal frequency (**A**) and Mycorrhizal intensity (**B**) of the five herbaceous plants in the different studied sites. Data are reported as mean (±SE) of three replicates per sample.

#### *3.3. Densities of Spore Populations in the Studied Soils*

The density of AMF spores isolated from the rhizosphere of the sampled plants varied between 856 (*M. truncatula*) in Site 1 and 81 spores/100 g of soil (*D. simplex*) in Site 6 (Figure 3). Therefore, spore density varied significantly (*p* < 0.001) among the studied sites; the maximum values were recorded in the rhizosphere of plants from sites 1, 3 and 5 (inside the conserved areas and lightly grazed) and the minimum values in sites 2, 4 and 6 (outside the conserved areas and intensively grazed). Bare areas (bulk soil), followed by *D. simplex* rhizosphere soil, always presented the lowest values for the six sites (Figure 3).

**Figure 3.** Distribution of AMF spores in the different rhizosphere soils. Data are reported as mean (±SE) of three replicates per sample.

There was a clear relationship between the intensity of root cortex AMF colonization (M%), mycorrhizal frequency (F%), and the density of AMF spores in the rhizospheres (Figures 2 and 3). Furthermore, the highest spore density was registered in the rhizosphere of the plant species with higher mycorrhiza frequency and intensity (*M. truncatula*).

#### *3.4. E*ff*ect of the Grazing Intensity on the Di*ff*erent AMF Parameters*

The grazing intensity strongly affected the mycorrhizal colonization and density of spores (Table 3), in accordance with the lowest values always registered, either for the plant mycorrhizal status or for AMF spore populations, in samples from sites outside the protected areas (Figures 2 and 3). This explains the negative effect of the grazing intensity on the different mycorrhizal parameters (Figure 4).

**Table 3.** Two-factor ANOVA analysis of the impact of grazing intensity on mycorrhizal parameters.


\*\*\* Significant at *p* < 0.001, \*\* Significant at *p* < 0.01, ns: non-significant.

**Figure 4.** Results of a Canonical Correlation Analysis (CCorA) for the relationships between grazing intensity and AMF properties. Sp. num: number of spores; F%: mycorrhizal frequency; M%: mycorrhizal intensity.

#### *3.5. E*ff*ect of Soil Properties on the Di*ff*erent AMF Parameters*

According to the results of the Principal-component analysis (PCA) (Figure 5), the highest values of AMF colonization and number of spores were found in Site 1 at Bou-Hedma followed by Site 3 at Zarat. The major physical and chemical parameters of these sites (Org. Mat, E.c, Wat. Con, T.N) were in a strongly positive correlation with the different mycorrhizal properties. In contrast, the total phosphorus (T.P) available negatively affected all the mycorrhizal parameters and was considered a limited factor to these parameters (Figure 5).

**Figure 5.** Results of a principal-component analysis (PCA) showing the relationships between soil properties and AMF parameters. E.c: electrical conductivity; T.N: total nitrogen; T.P: total phosphorus; Org. Mat: organic matter; Wat. Con: water content; Sp. num: number of spores; F%: mycorrhizal frequency; M%: mycorrhizal intensity.

#### *3.6. Microbiological and Biochemical Properties of the Soils and Impact of AMF on Soil Microbial Communities*

We found significant effects (*p* < 0.001) of the mycorrhizal plants on Cmic (Figure 6), this meaning that plants are affecting soil microbial biomass and consequently potential soil dynamics. The Cmic values were lower in bare areas (control soil) than in rhizosphere soils, the maximum values being recorded in the rhizospheres of the three Fabaceae plants *M. truncatula, A. corrugatus* and *L. creticus*. Cmic also markedly decreased in the intensively grazed sites (outside the protected areas), following the same trend of mycorrhizal frequency and intensity in root samples and of number of spores in rhizosphere soils (Figures 2, 3, 6 and 7).

**Figure 6.** Microbial biomass carbon (Cmic) in the different sampled soils. Data are reported as mean (±SE) of three replicates per sample.

**Figure 7.** Phosphatase (**A**), dehydrogenase (**B**), β-glucosidase (**C**) activities in the studied soils. PNP: p nitrophenol, INTF: iodonitrotetrazolium formazan. Data are reported as mean (±SE) of three replicates per sample.

The highest activities of the three enzymes evaluated (phosphatase, dehydrogenase and β-glucosidase) (Figure 7) were observed in the rhizospheres of the three Fabaceae plants *M. truncatula*, *A. corrugatus* and *L. halophilus*, following the same trend of Cmic for each plant species and each site (Figure 6) as well as the values of mycorrhizal frequency and intensity and number of spores (Figures 2, 3 and 7), evidencing activities consistently greater in lightly grazed sites in comparison with those that were intensively grazed (Figures 2, 3, 7 and 8). The ANOVA analysis evidenced significant differences in these values.

**Figure 8.** Schematic representation of the AMF effects on soil fertility and ecosystem stability under lightly grazed sites (inside the protected areas) versus intensively grazed sites (outside the protected areas), for arid Mediterranean ecosystems.

A clear relationship was observed between the AMF parameters (frequency and intensity of colonization of roots and number of AMF spores in the rhizosphere) and microbial activity (microbial biomass and enzyme activities) (Table 4, Figures S1 and S2). Based on the Pearson correlation coefficient, positive correlations with a highly significant value (*p* < 0.001) were observed between all AMF parameters (F%, M% and number of spores) and microbial activity (microbial biomass and enzyme activities: dehydrogenase, phosphatase and the β-glucosidase). A lower significant value (*p* < 0.05) was only observed in Site 6 (outside Oued Dkouk Natural Reserve) between the frequency of AMF colonization and two microbial parameters: microbial biomass and activity of β-glucosidase (Table 4).


**Table 4.** Pearson correlation coefficient between frequency of mycorrhization (F%), intensity of mycorrhization (M%) and density of AMF spores (Sp. num) and microbial activity expressed by microbial biomass (Cmic), phosphatase activity, dehydrogenase activity, and β-glucosidase activity.

\*\*\* Significant at *p* < 0.001, \* Significant at *p* < 0.05.

The results shown above, concerning the effect of grazing on AMF colonization (Table 3, Figure 4) and the impact of AMF on microbial communities and activities (Table 4, Figures S1 and S2), are summarized in Figure 8, which evidences the importance of AMF on soil functionality under arid Mediterranean ecosystems. Significant mycorrhizal colonization of roots was registered in lightly grazed sites, leading to high levels of microbial communities expressed by high values of biochemical activities. All these necessarily enhances soil fertility, better-adaptation and growth of plants and finally, in the stability of the ecosystem (Figure 8).

#### **4. Discussion**

Mediterranean arid and semi-arid ecosystems are characterized by high temperatures and drought for most of the year. These conditions limit plant establishment and growth and accelerate soil degradation and change microbial communities [1]. Mycorrhiza form communication pathways between plants and soil, influencing plant nutrient cycling, and restoring and maintaining soil fertility, thus influencing the microbial communities of the rhizosphere and extending the influence of plants to the soil. Several works emphasized the role of AMF in sustaining plant cover in semi-arid and arid ecosystems [4,8,10] as is the case of the following Mediterranean conserved areas in Tunisia: Bou-Hedma National Park, Zarat protected area and Oued Dkouk Natural Reserve.

Under natural conditions, about 90% of the terrestrial plants are mycorrhized, and AMF are found in all climates and ecosystems [10]. Due to their role on plant nutrition and defense, and to the importance of the extra-radical mycelium in the establishment of biological networks, AMF are determinant for the establishment and sustainability of plant communities and environment functioning [7]. However, the diversity of the AMF community and the intensity of arbuscular mycorrhizal colonization of natural vegetation is dependent on the availability of AMF spores and the mycorrhizal dependency of the plant species [24,25], as well as on the soil structure and management [26]. All these statements agree with the results of the present study, which evidenced colonization by AMF for all the herbaceous species, except for the Cruciferaceae *Diplotaxis simplex* (Figure 2). This was not surprising since Cruciferaceae are usually recognized as non-mycorrhizal plants [27]. However, distinct AMF colonization rates and intensities were observed for each plant species, which may be related to different levels of mycorrhizal dependence [28], and/or the availability of AMF spores. Higher mycorrhization rates were observed in the legume species *M. truncatula, A. corrugatus* and *L. creticus*, which presented high levels of mycorrhizal intensity, independent of the sampled site (Figure 2). This may highlight their mycorrhizal dependency and high demand for phosphorus (P) in comparison with plants from other families such as Poaceae [29].

In general, and particularly in semi-arid and arid ecosystem, AMF vary greatly with soil characteristics. Several biotic and abiotic factors may contribute to the distinct mycorrhizal intensities observed for the same plant species [30,31]. One abiotic factor that severely interferes with mycorrhization is the concentration of available phosphorus [10,32,33]. In the present study, soil phosphorus concentrations varied between sites, ranging from 7 to 15.3 ppm, and could explain the distinct mycorrhizal intensities observed for the same species in distinct sites (Figures 2 and 5), supporting that its availability is a crucial driver of mycorrhizal communities and activities [34].

Another factor frequently cited as responsible for lower levels of mycorrhiza formation is the availability of AMF spores, which is known to greatly vary in the ecosystems [35]. Sporulation in AMF occurs when the development of the mycelium begins to be limited by nutrients and is a highly carbon demanding process. This may explain why the number of AMF spores in the rhizosphere varies among plant species and, for the same plant species, among sites [36–38]. Apart from the rate of spore formation, the number of AMF spores in a soil also depends on the rates of spore germination and degradation [39]. As AMF are obligate biotrophs, the number of spores and propagules tends to be higher in the rhizosphere than in the bulk soil [12,40], and higher in the rhizosphere of plants with a higher intensity of AMF colonization (Figure 3). These differences are particularly evident in arid and semi-arid soils with high organic matter turnover rates and low organic matter content [41].

In low fertility soils, mycorrhiza are highways for nutrient and water transport, expanding the plant root system and the volume of soil exploited by the plant. However, similarly to roots, hyphae are leaky and lose nutrients into the hyphosphere, which will promote the development of selected microorganisms. In this context, it is expected that AMF play a crucial role in the biological characteristics of the rhizosphere [22,42] (Figure 6). The high values of Cmic observed in the rhizosphere soil of mycorrhizal plants imply that mycorrhiza contribute to improve the availability of carbon substrates to the microbial community of the rhizosphere [43,44]. These results support the hypothesis that AMF establish unique interactions with plant roots, conferring special characteristics to the rhizosphere [45], where several by product-based symbiosis and microbial loops may be assembled, contributing to improved carbon use efficiency. In this particular respect, the importance of AMF in promoting the development of the microbial community is confirmed by the low and comparable levels of microbial carbon observed in the soils without plant cover (control) and in the rhizosphere soil of the non-mycorrhizal plant *D. simplex* (Figure 6).

AMF increase the diversity of the carbon sources available to the microorganisms in the rhizosphere [46]—which is partly due to their nutritional mode; the excretion of catabolic enzymes to the surrounding medium, and to the direct access by the AMF to the plant carbon. Phosphatase and β-glucosidase are two of those catabolic enzymes, and their activities were in fact higher in the rhizosphere of the mycorrhizal plants than in that of the non-mycorrhizal plant studied or in the bulk soil (Figure 5). The importance of the soil microbial activity in association with several enzymes' activities was highlighted by the similarity in the activity patterns of the two hydrolytic enzymes (phosphatase and β-glucosidase) as well as of dehydrogenase, indicators of microbial activity [47], which can be inferred as being decreased under high grazing intensity (Figures S1 and S2).

At this point, it is clear that plant and AMF species are important modulators of rhizosphere characteristics (Figures 6–8). What is not clear is why biological indicators of soil characteristics are so distinct in the three sites outside the conserved areas. One key factor common to these sites may be limiting and related to mycorrhizal colonization. Grazing may affect mycorrhization (Table 3, Figure 4) directly in nutrient-limited ecosystems through a direct competition for carbon (Figure 8). Grazing was consistently associated with lower soil organic matter and increased P concentrations showing a strong impact on soil chemical properties, able to negatively influence mycorrhization in a direct way through the increase of P bioavailability. But in nutrient-limited ecosystems grazing per se may also affect mycorrhization (Table 3, Figure 4).

Plants are the primary producers of ecosystems, obtaining their biomass and energy from the carbon fixed by photosynthesis (Figure 8). Depending on the plant species and on the growth conditions, 20–50% of the newly photosynthetically fixed carbon may be transported to the roots and lost as rhizodeposition. In arid soils, this carbon is the main source of energy and biomass building material for the microbial community. Therefore, increasing grazing intensity decreases the carbon availability from plants and, consequently, also decreases the carbon available for mycorrhization [13,14]. An increase in grazing intensity has impacts on soil biological characteristics, with a significant decrease in the mycorrhizal intensity, as evidenced in the present study, for the three mycorrhizal plant species sampled in the sites outside the conserved areas (Table 3, Figure 5; Figure 8).

#### **5. Conclusions**

AMF is of vital importance in the arid Tunisian ecosystem, where soils are generally poor in organic matter and nutrients. High mycorrhization levels and high spore densities were revealed across the different sites, with particular incidence in the sites inside the conserved areas. Therefore, grazing intensity affected the potential beneficial influence of AMF on soil microbial processes. At this point, the importance of AMF to ecosystem dynamics is clear. What is not clear is how AMF are adapted to grazing intensity to allow ecosystem stability, particularly for situations that may be represented by the three sites outside the conserved areas. Therefore, future research should be focused on determining the ways through which AMF adapt to grazing intensity. This may be a key feature of

ecosystem management, considering that AMF could be an important alternative for sustaining soil quality, and could be exploited as potential inoculants for rehabilitation and restoration programs in Mediterranean ecosystems.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/6/234/s1, Figure S1: Linear regression to model the correlations between mycorrhizal proprieties (mycorrhizal frequency, mycorrhizal intensity, and number of spores) and the microbiological parameter microbial biomass carbon. The analyses included all available data: from the rhizospheres of the four plants and from bulk soil, Figure S2: Linear regression to model the correlations between mycorrhizal proprieties (mycorrhizal frequency, mycorrhizal intensity, and number of spores) and biochemical activities: phosphatase (in <sup>μ</sup>g PNP. g−<sup>1</sup> soil·h<sup>−</sup>1), dehydrogenase (μg INTF. g−<sup>1</sup> soil·h<sup>−</sup>1) and <sup>β</sup>-glucosidase (μg PNP. g−<sup>1</sup> soil·h<sup>−</sup>1). The analyses included all available data: from the rhizospheres of the four plants and from bulk soil.

**Author Contributions:** Conceptualization, M.M. (Mohamed Mars); methodology, N.M.; validation, M.M. (Mosbah Mahdhi) and T.D.; formal analysis, N.M.; investigation, N.M.; data curation, N.M.; writing—original draft preparation, N.M.; writing—review and editing, C.C. and M.F.C.; visualization, C.C. and M.F.C.; project administration, M.M. (Mohamed Mars); funding acquisition, C.C. and M.F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Ministry of High Education and Research Development-Tunisia and by FCT/MCTES for the financial support to cE3c (project UIDB/00329/2020) and CESAM (projects UIDP/50017/2020+UIDB/50017/2020), through Portuguese national funds.

**Acknowledgments:** We thank the personal staff of Bou-Hedma National Park, of Zarat-Gabes Protected area and of the Natural reserve of Oued Dkouk, for facilities and support for root and soil sampling. Thanks, are also due to Francisco Caeiro for providing helpful comments and English revision.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Fungal Diversity in the Phyllosphere of** *Pinus heldreichii* **H. Christ—An Endemic and High-Altitude Pine of the Mediterranean Region**

#### **Jelena Lazarevi´c 1,\* and Audrius Menkis <sup>2</sup>**


Received: 30 March 2020; Accepted: 21 April 2020; Published: 28 April 2020

**Abstract:** *Pinus heldreichii* is a high-altitude coniferous tree species naturaly occurring in small and disjuncted populations in the Balkans and southern Italy. The aim of this study was to assess diversity and composition of fungal communities in living needles of *P. heldreichii* specifically focusing on fungal pathogens. Sampling was carried out at six different sites in Montenegro, where 2–4 year-old living needles of *P. heldreichii* were collected. Following DNA isolation, it was amplified using ITS2 rDNA as a marker and subjected to high-throughput sequencing. Sequencing resulted in 31,831 high quality reads, which after assembly were found to represent 375 fungal taxa. The detected fungi were 295 (78.7%) Ascomycota, 79 (21.0%) Basidiomycota and 1 (0.2%) Mortierellomycotina. The most common fungi were *Lophodermium pinastri* (12.5% of all high-quality sequences), *L. conigenum* (10.9%), *Sydowia polyspora* (8.8%), *Cyclaneusma niveum* (5.5%), Unidentified sp. 2814\_1 (5.4%) and *Phaeosphaeria punctiformis* (4.4%). The community composition varied among different sites, but in this respect two sites at higher altitudes (harsh growing conditions) were separated from three sites at lower altitudes (milder growing conditions), suggesting that environmental conditions were among major determinants of fungal communities associated with needles of *P. heldreichii*. Trees on one study site were attacked by bark beetles, leading to discolouration and frequent dieback of needles, thereby strongly affecting the fungal community structure. Among all functional groups of fungi, pathogens appeared to be an important component of fungal communities in the phyllosphere of *P. heldreichii*, especially in those trees under strong abiotic and biotic stress.

**Keywords:** needle pathogens; high altitude forests; DNA metabarcoding; Montenegro

#### **1. Introduction**

*Pinus heldreichii* is a high-altitude (grows at ca. 1200–2000 m) conifer tree species with a discontinuous and restricted distribution in the Mediterranean region. Its forests are naturaly regenerated and consists of small and disjuncted populations located in high mountain areas influenced by the Mediterranean climate in the Balkans and southern Italy. Although in the past *P. heldreichii* formed a continuous forest belt in the Balkans, currently its forests are scattered and largely isolated. Growing primary on shallow calcareous soils, it inhabits typical tree line habitats, often on steep ridges, mountain sides and screens. Such habitats are nutrient poor, exposed, dry and cold during the winter [1,2]. Larger *P. heldreichii* forests can still be found on mountain plateaus or in valleys situated at altitudes of ca. 1200–1300 m that are characterised by more developed soils (cambisols). Those forests are considered to be climazonal. i.e., characterised by permanent and stable vegetation forests of *P. heldreichii*, that are growing at ecological optimum for the species [3]. Being a tertiary relic, *P. heldreichii* evolved to survive

severe winter frosts, and short, dry and warm summers with intensive sun radiation. Due to the short growing season, *P. heldreichii* grows very slowly, reaching ca. 15 m height after about 150 years [1]. It develops thick and branched roots, which penetrate deep in cracks of calcareous stones [2]. In the last few centuries, *P. heldreichii* forests were affected by intensive livestock grazing, exploitation and forest fires [4]. Intensive grazing, as well as the short and dry growing season have also resulted in very limited natural regeneration [2–4], though recent observations show some regeneration in the abandoned mountain areas [5].

*Pinus heldreichii* is a protected species both in Balkan countries and in Italy due to the key importance of nature conservation, protection against gravitational natural hazards, landscape conservation and recreation. Hence, *P. heldreichii* requires special attention, i.e., the development and application of conservation measures [1–4,6].

Fungi represent the largest microbial component associated with forest trees. They play key roles in forest ecosystems, especially in pines forest that are obligatory mutualistic, and are important contributors to the primary production and carbon, nutrient and water cycling [7–9]. Pathogenic fungi may negatively affect health and growth of forest trees [10,11], while fungal endophytes and epiphytes support ecological adaptations of host plants and constitute an important component of microbial biodiversity [12,13]. However, information about fungal communities associated with *P. heldreichii* is limited. Previous studies on needle pathogens of *P. heldreichii* have focused on either those that were affecting natural regeneration [14–16], or potentially invasive ones such as *Dothistroma septosporum* [17]. Information about the occurrence of *D. septosporum* in stands of *P. heldreichii* was not known until recently [17], what indicated the need for a wider assessment of fungal diversity, including fungal pathogens associated with the phyllosphere of this tree species.

The aim of this study was to assess the diversity and composition of fungal communities in living needles of *P. heldreichii* specifically focusing on fungal pathogens. Needles were sampled across the natural distribution range of *P. heldreichii* in Montenegro, including sites under different environmental conditions. This was expected to demonstrate potential site-specific effects of environmental conditions on fungal diversity and community composition. By using high-throughput sequencing of fungal ITS2 rDNA, we examined fungal communities in 2–4 year-old living needles of *P. heldreichii* from six sites situated in four mountain areas in Montenegro.

#### **2. Materials and Methods**

#### *2.1. Study Sites*

Samples were taken at six sites that were situated in four mountain areas, representing the *P. heldreichii* distribution range in Montenegro (Figure 1). The sampling sites were Orjen (ORJ), Prekornica (PRE), Kuˇcka Korita (KKO), Kuˇcka korita North (KKN), Kuˇci MT (KMT) and Prokletije (PRO) (Figure 1, Table 1).

The sampling sites differed in altitude, morphology of terrain and soil properties (Table 1, Figure 2). Sites at ORJ and KMT were on altitudes of ca. 1800 m, with trees growing on stone ridges, high inclinations, and lithosols. Sites at PRE, PRO and KKO were on an altitude of ca. 1250 m, on flat terrain and more developed soil that was leptic cambisol. A site at KKN was in close proximity to the KKO, but it was on a slope with a thinner soil layer (molic leptosol). Trees at the KKN were in groups heavily infested by *Tomicus* sp. bark beetles (Scolitinae). Trees at ORJ and KMT were growing at the upper tree line, and thus, these sites were considered as sites with harsh growing conditions, while other sites had moderate growing conditions.

**Figure 1.** A map showing (**a**) distribution of *Pinus heldreichii* in the Mediterranean region; a grey area was added and represents an additional *P. heldreichii* site; the map is adapted from [18] and (**b**) sampling sites in Montenegro (ORJ-Orjen; PRE-Prekornica; KKO-Kuˇci Mt, Kuˇcka korita; KKN-Kuˇci Mt., Kuˇcka korita North, KMT- Kuˇci mountain, Momonjevo, PRO- Prokletije); a map was adapted from [19].


**Table 1.** Characteristics of *Pinus heldreichii* sampling sites.

<sup>a</sup> WRB soil clasification system [20,21], <sup>b</sup> Koopen climate classification [22,23], **<sup>c</sup>** Health status of the trees.

**Figure 2.** (**a**) *Pinus heldreichii* shoots with needles; trees growing at (**b**) Kuˇci Mountains, Momonjevo (KMT), and (**c**) Kuˇci Mountains, Kuˇcka korita (KKO) sites.

Climate at the study sites is characterised as a humid warm temperate climate type (Cf), which is represented by subtypes Cfsb, Cfs"b and Cfb [22,23]. Summers are short, dry and chilly, and winters are cold and windy. At different study sites, the mean annual air temperature is between 1 ◦C (at PRO) and 6 ◦C (at ORJ). The mean daily summer maximum is between 7 and 11 ◦C, and the mean daily winter minimum is between −5 and 0 ◦C. The absolute minimum is in winter at ca. −30 ◦C and the maximum is in summer, at ca. 30 ◦C. The mean annual precipitation is highest at ORJ with ca. 3800 mm, and lowest at PRO with ca. 2000 mm. The rainfall reaches maximum in late autumn and early winter, while the minimum is during the summer months that is often followed by 40–70 day-long periods of droughts. The mean summer precipitation is 300 mm at ORJ and ca. 220 mm on other sites, making ca. 10% of the total annual precipitation [22].

The bedrock at the study sites is solid chalk limestone, which contains a small proportion of insoluble residues. In general, soils are poorly developed, very water porous, skeletal leptosols (rendzina), having only A horizon, characterized by accumulation of humus (A-R profile). According to the classification of texture, the soil is a sandy loam and the structure is powdery. Soil is faintly differentiated down in the soil profile. Soil is rich in humus content (10–25%), poor in calcium-carbonate and has a weak acid reaction (pH ca. 6). On more flat terrains, leptosol develops into the leptic cambisol, which contains initial B horizon. Brown B horizon improves the water retention capacity of soil [20,24].

#### *2.2. Experimental Design and Sampling*

At the sampling sites, stands of *P. heldreichii* were healthy-looking (Figure 2a). An exception was the KKN site, where trees were damaged by *Tomicus* sp. bark beetles, resulting in discolouration and fungal infection of needles. At all sites, the shoot dieback was occasional, but with higher rates at the KKN site.

Sampling was carried out in May 2015, i.e., before the beginning of the growing season. At each site, experimental design included sampling of five twigs with needles from five different trees. Hence, at each site, five mature *P. heldreichii* trees situated at a distance of ca. 50 m from each other were selected and five twigs of up to 15 cm long and up to 2 m from the ground were sampled from different parts of the crown using secateurs. This sampling approach was used in order (i) to get a joint representative sample per each site (all samples per site were amplified with the same barcode, see below); (ii) to compare fungal communities among different sites, but not within the same site. Collected twigs were placed in plastic bags and transported to the laboratory.

In the laboratory, twigs and needles were assessed for the presence of disease symptoms. The current-year needles were typically green and without disease symptoms. Symptomatic needles were generally two years old or older and symptoms included changed colour that was yellow, orange, pale green or yellowish with red bands. Other symptoms were dying needle tips and a necrotic base of needles. In the case of insect attack (KKN site), defoliation and necrotic needles were common. In order to sample the potentially entire fungal community associated with needles of *P. heldreichii*, both healthy looking needles and needles with disease symptoms were selected. No surface sterilization was carried out. Following morphological examination, up to 5 representative needles that were 2–4 years old, were selected per twig. Selected needles were cut into 0.5–1 cm long segments, containing a random mixture of both asymptomatic and symptomatic parts of needles, placed in 2 mL screw cap tubes, and stored frozen at −20 ◦C before DNA extraction.

#### *2.3. DNA Isolation, Amplificationand Sequencing*

DNA extractions were done from 150 samples (6 sites × 5 trees × 5 needle samples), which previously were freeze-dried for 48 h. For isolation of total DNA, needles were homogenised in a Fastprep machine (Precellys, Montigny-le-Bretonneux, France). The extraction was completed using CTAB protocol [25]. After extraction, DNA samples were purified using a JetQuick DNA purification kit (Genomed GmbH, Leinfelden, Germany). The DNA concentration of each sample

was determined using a NanoDrop™ One spectrophotometer (Thermo Scientific, Rodchester, NY, USA) and adjusted to 1–10 ng/μL. Amplification by PCR of the ITS2 rDNA region was done using barcoded fungal-specific primer gITS7 [26] and barcoded universal primer ITS4 [27]. All 25 samples from the same site were amplified using primers with the same barcode, resulting in six different barcodes representing each site. Amplification of several samples with the same barcode was done to get a broader representativeness of fungal communities per site. Amplifications were performed using the Applied Biosystems 2720 thermal cycler (Foster City, CA, USA). An initial denaturation step started at 95 ◦C for 2 min, followed by 27 amplification cycles of denaturation at 95 ◦C for 30 s, annealing at 55 ◦C for 30 s, and extension at 72 ◦C for 60 s. The thermal cycling was ended by a final extension step at 72 ◦C for 7 min. The PCR products were analyzed using gel electrophoresis on 1% agarose gels stained with Nancy-520 (Sigma-Aldrich, Sweden). PCR products were purified using a sodium acetate protocol [28]. Purified PCR products were quantified using a Qubit fluorometer 4.0 (Thermo Fisher Scientific, Waltham, MA, USA), and an equimolar mix of all PCR products with six barcodes (one per each site) was used for high-throughput sequencing using a Pacific Biosciences platform (Menlo Park, CA, USA) at the SciLifeLab (Uppsala, Sweden).

#### *2.4. Bioinformatics*

Principles of bioinformatics followed that specified in [29]. The sequences obtained from the six samples that represented six sampling sites (Table 2), were subjected to quality control and clustering in the SCATA NGS sequencing pipeline [30]. The initial procedure started with quality filtering of the sequences that included the removal of sequences shorter than 200 bps, sequences with low read quality, primer dimers and homopolymers, which were collapsed to 3 bps before clustering. Only sequences containing a barcode and primer were retained. Then, the primer and sample barcodes were removed from the sequence, but information on the sample and sequence association was stored as meta-data. A single-linkage clustering based on 98% similarity was used to cluster sequences into different taxa. For each cluster, the sequence of the most common genotype was used for taxonomic identification. For clusters containing only two sequences, a consensus sequence was produced. The taxa were taxonomically identified using the GenBank database and the Blastn algorithm [31]. The following criteria were used for identification: sequence coverage >80%; 94–97% similarity to genus level and >98% similarity to species level. Sequences deviating from these criteria were identified only to a high taxonomic rank and were given unique names as shown in Table 3 and Table S1. Representative sequences of fungal nonsingletons are available from GenBank under accession numbers MT241905–MT242268.

**Table 2.** Generated high-quality ITS2 rDNA fungal sequences and detected diversity of fungal taxa at different sampling sites of *Pinus heldreichii*. Within the column *No. of Fungal Taxa*, values followed by the same letter in chi-square test do not differ significantly at *p* > 0.05.


**Table 3.** Relative abundance of the 25 most common fungal taxa detected in needles of *Pinus heldreichii* at six sites in Montenegro. Sites are as in Table


\*

A—Ascomycota,

B—Basidiomycota.

 1.

#### *2.5. Statistical Analyses*

Rarefaction analysis was carried out using Analytical Rarefaction v.1.3 [32]. Differences in richness of fungal taxa in different study sites of *P. heldreichii* were compared by nonparametric chi-square testing [33]. As each of the datasets was subjected to multiple comparisons, confidence limits for *p*-values of chi-square tests were reduced the corresponding number of times as required by the Bonferroni correction [34]. The Shannon diversity index, qualitative Sorensen similarity index and nonmetric multidimensional scaling (NMDS) in Canoco 5 [33,35,36] were used to characterize the diversity and composition of fungal communities.

#### **3. Results**

The phyllosphere fungi of *P. heldreichii* were examined from four different mountain regions, represented by six sampling sites across the distribution range of *P. heldreichii* in Montenegro (Figure 1). Amplification of fungal ITS2 rDNA from 150 needle samples, PacBio sequencing and quality filtering resulted in 31,829 high quality reads. Sequence assembly and BLASTn analyses showed that the fungal community in the phyllosphere of *P. heldreichii* was composed of 375 fungal taxa (Table 2, Table S1). Nonfungal taxa and singletons were excluded. The detected fungi were 295 (78.7%) Ascomycota, 79 (21.0%) Basidiomycota and 1 (0.2%) Mortierellomycotina. Identification at least to genus level was successful for 254 (67.7%) out of 375 fungal taxa (Table S1), and those represented 78.5% of all fungal sequences.

In different study sites, the number of fungal taxa varied between 104 and 192 (Table 2). The chi-square test showed that the largest difference in richness of fungal taxa was between KKO and the remaining sites (Table 2). Rarefaction showed that fungal taxa detected in all sites did not reach the species saturation (Figure 3).

**Figure 3.** Rarefaction curves showing the relationship between the cumulative number of fungal taxa and the number of ITS2 rDNA sequences from needles of *Pinus heldreichii* from six different sites.

Information on the 25 most common fungal taxa representing 79.0% of all fungal sequences is in Table 3. The most common fungi in the phyllosphere of *P. heldreichii* were *Lophodermium pinastri* (12.5% of all fungal sequences), *Lophodermium conigenum* (10.9%), *Sydowia polyspora* (8.8%), *Cyclaneusma niveum* (5.5%) and Unidentified sp. 2814\_1 (5.4%) (Table 3). The most common fungal pathogens of pine needles were *S. polyspora* (8.9%), *C. niveum* (5.5%), *Neocatenulostroma germanicum* (2.3%)*, Allantophomopsiella pseudotsugae* (2.1%) and *Cenangium acuum* (0.7%) (Table 3). A pine needle pathogen, *Dothistroma septosporum* (0.25%) was also detected, but at lower relative abundance (Table S1). The other detected pathogenic species that are mainly known as wound pathogens of deciduous trees and/or agricultural

crops, were *Phaeomoniella* 2814\_15 (3.0%), *Ramoconidiophora euphorbiae* (1.5%), *Collophorina sp.* (0.9%)*, Geastrumia sp.* (1.4%) and *Athelia acrospora* (0.7%). The detected fungal endophytes that are known to produce antimicrobial metabolites were *Phaeosphaeria pontiformis* (4.4%), *Microsphaeropsis olivacea* (2.2%)*, Lachnellula calyciformis* (2.0%) and *Mollisia ligni* (0.7%). Ubiquitous saprotrophs were rare (Table S1).

The community composition of the phyllosphere fungi varied among different sites (Figure 4). Dothideomycetes dominated fungal communities at PRO and KKN sites, while Leotiomycetes dominated at KMT and PRE sites (Figure 4). At ORJ and KKO sites, the relative abundance of these two classes was similar (Figure 4). Sordariomycetes showed higher relative abundance at ORJ (22.6%), while all other fungal classes in different sampling sites were less abundant (Figure 4). Agaricomycetes (Basidiomycotina) were rare and their relative abundance at different sites varied between 0.03% and 6.3% (Figure 4).

**Figure 4.** Relative abundance of different fungal classes in needles of *Pinus heldreichii* from six different sites in Montenegro. Classes comprising <1.0% of relative abundance for a given set of sequences were combined and shown as other. Sites are as in Table 1.

The NMDS of fungal communities associated with needles of *P. heldreichii* (Figure 5) showed that KKO, PRO and PRE sites clustered together and were separated from the remaining sites along axis 1. Fungal communities at the ORJ and KMT showed a closer proximity (Figure 5), while fungal communities at the KKN site differed from all remaining sites and were separated along both axis 1 and 2.

At ORJ and KMT sites that are characterized by harsh growing conditions, the fungal community was dominated by *Lophodermium pinastri*, *Sydowia polyspora*, *Phaeomoniella* sp., *Allanthopomopsiella pseudotsugae*, *Lachnellula calyciformis*, *Ramoconidiophora euphorbiae*, *Collophorina* sp., *Geastrumia* sp. and *Mollisia ligni* (Table 3). At PRE, PRO and KKO sites that are characterized by moderate growing conditions, fungal community was dominated by *Lophodermium conigenum*, *L. pinastri*, *Cyclaneusma niveum*, *Phaeosphaeria punctiformis*, Unidentified sp. 2814\_11, *Cenangium acuum* and *Dothistroma septosporum*. At the KKN, with insect attacks and partly necrotic needles, fungal community was dominated by *S. polyspora*, *Microsphaeropsis olivacea*, *Neocatelunostroma germanicum*, Uncultured sp. 2814\_10, *Chaetothyriales* 2814\_18 and *Geastrumia* sp. (Table S1).

**Figure 5.** Ordination diagram based on nonmetric multidimensional scaling (NMDS) of fungal communities associated with needles of *Pinus heldreichii* from six different sites in Montenegro, with 56.5% variation on Axis 1 and 28.2% on Axis 2. Each point in the diagram represents a single site and the size of the point reflects richness of the fungal taxa as in Table 2. Sites are as in Table 1.

The Shannon diversity index of fungal communities in different study sites was between 2.7 and 3.4 (Table 2). Among different sites, Sørensen similarity index of fungal communities was low to moderate as it ranged between 0.28 and 0.54 (Table S2).

#### **4. Discussion**

A recent development of high-throughput sequencing methods provides powerful tools to explore fungal diversity. Such tools enable identification of complex fungal communities and individual community components directly from environmental samples. Besides, while providing detailed and semi-quantitative information, these tools enable studying the effects of different factors on fungal diversity and community composition [29,37,38]. By using PacBio sequencing, we detected nearly 400 fungal taxa, which were associated with needles of *P. heldreichii* including fungi present in very small abundances (Table S1). Many fungal taxa remained unidentified, which remains a major challenge in fungal taxonomy [39]. However, while using high-throughput sequencing, we need be aware of potential limitations that can include methodological biases, limitations of markers and bioinformatics challenges [40,41].

The observed diversity of the phyllosphere fungi can be considered being high and comparable with similar studies on different *Pinus* species. For example, there were 446 and 260 fungal taxa recorded during studies of needle-associated fungi of *Pinus sylvestris* in Sweden and Poland, respectively [39,42]. Furthermore, by using fungal culturing, there were 118 fungal taxa associated with needles of *Pinus taeda* in USA [43], 35 taxa associated with needles of *Pinus halepensis* [44] and *P. sylvestris* in Spain [45]. According to [45], the species diversity of fungal endophytes varied from 21 in *Pinus monticola* to 49 in *Pinus nigra*. At different study sites, the growing conditions of *P. heldreichii* had in general a minor effect of the absolute richness of fungal taxa, which was similar among different sites (Table 2). However, with respect to the number of sequences obtained from samples of each site, lower richness of fungal taxa was at the KKO (Table 2), what was likely due to the dominance of several fungal taxa such as, e.g., *Lophodermium* species, *Cylaneusma niveum,* Unidentified sp. 2814\_11, *Phaeosphaeria pontiformis* and Unidentified sp. 2814\_1. Several fungal taxa of the present study were recorded for the first time in the Balkan region such as, e.g., pathogens *Neocatelunostroma germanicum* and *Allantophomopsiella pseudotsugae*. The study also revealed a number of fungal pathogens that were previously known

from agricultural crops as, e.g., *Phaeomoniella*, *Ramoconidiophora*, *Geastrumia* and *Athelia.* These could have been overlooked in other studies due to latent occurrence [46]. As *P. heldreichii* belongs to a group of forest trees that grows under harsh enviromental conditions, the possibility should not be excluded that both the host and the specific environment influence the composition of associated fungal communities. Studies on plant–endophyte associations in high stress habitats have revealed that at least some fungal endophytes can contribute to stress tolerance of host plants [12]. Indeed, the results of this study revealed site-specific differences in fungal communities associated with needles of *P. heldreichii*. For example, fungal communities were more similar among sites situated at lower altitudes, i.e., with moderate growing conditions (PRE, PRO, KKO) as compared to those sites at higher altitudes, i.e., with harsh growing conditions (ORJ, KMT) (Table 1, Figure 5). Further, the fungal community at the KKN site was largely different from the remaining sites, what was likely due to insect damage to trees, resulting in mainly symptomatic and partly necrotic needles. These observations demonstrate that the composition of fungal communities and their succession in needles of *P. heldreichii* can be determined by different abiotic and biotic factors. Different environmental factors have earlier been shown to play an important role in shaping fungal communities both in the phyllosphere and in the rhizosphere of forest trees [5,42,47,48].

Fungal pathogens represent an important biotic factor that may negatively affect health and growth of forest trees [11]. In the present study, a fungal pathogen *Sydowia polyspora* was among the most dominant fungi (Table 3). *Sydowia polyspora* has a wide geographical range [48] and is common in Europe. The pathogenicity of *S. polyspora* to young conifers (genera *Thuja*, *Abies*, *Tsuga*, *Larix*, *Picea* and *Pinus*) was previously reported [49]. In the Balkan region, it has been reported in forest plantations of *P. nigra* and *P. sylvestris* in Serbia, occurring on needles damaged by drought or frost [50]. It was also detected on needles of *P. halepensis* in Italy and Spain [44,51], *P. sylvestris* in Poland and Lithuania [29,39] and on *P. ponderosa* in North America [52]. *Sydowia polyspora* has recently been reported as a pathogen that dominates fungal communities vectored by bark beetles associated with *P. radiata*, *P. nigra* subsp. *salzmannii* and *P. sylvestris* in Spain [53] and *Pinus yunnanensis* in China [54]. *Sydowia* symptoms include needle necrosis and shoot dieback. The fungus is favoured by a warm climate, especially if the host is stressed by summer drought or insect or mite attack [53]. Indeed, in the present study the highest relative abundance of *S. polyspora* was at the KKN (Table 2), the site that was subjected to insect attack. Besides, it appears that *S. polyspora* had a negative effect on the relative abundance of *L. pinastri* (Table 2), which is in agreement with Behnke-Borowczyk et al. [39], who have observed the similar pattern on *P. sylvestris*. *Cyclaneusma niveum* was also among the dominant fungi (Table 3). It is one of the two fungal pathogens causing Cyclaneusma needle cast, which is an important needle disease reported from many pine species including *P. nigra* and *P. sylvestris* in the Balkan region and in Crimea [50,55,56], and reported from *P. halepensis* in Spain [44]. It was suggested that *C. niveum* is more frequently associated with *P. nigra* under warmer climate conitions, while *Cyclaneusma minus* is favoured by wet, humid, above-freezing conditions, and thus, more commonly infects *P. sylvestris* grown in central and northern Europe [39,50,57]. Our study is in agreement, as on *P. heldreichii*, *C. niveum* dominated the fungal community, while *C. minus* was recorded at a low frequency (Table S1). Interestingly, *C. niveum* was more common on sites with prevailing moderate than harsh enviromental conditions. *Neocatenulostroma germanicum* is a recently identified fungal pathogen causing needle blight on *P. mugo*, *P. sylvestris* and *P. nigra* subsp. *pallsiana* in Lithuania, Ukraine and Poland, where it was found to be commonly associated with *Dothistroma*, *Lecanosticta acicola* and Cyclaneusma needle cast infections [39,58]. In this study, it was present in the needles of *P. heldreichii* from five sites, with considerably higher relative abundance on symptomatic needles at the KKN site and the neighbouring KKO site as compared to the remaining sites. *Allantophomopsiella pseudotsugae*(syn. *Phomopsis pseudotsugae*) was a pathogen commonly detected on needles of *P. heldreichii* grown on sites under harsh growth conditions (ORJ and KMT) (Table 3). While being a pathogen of conifers, that is, mainly infecting pines, it develops on young shoots [59] and has been reported in several European countries [60]. The present finding of the fungus represents a new record for this part

of Europe, indicating that *A. pseudotsugae* could be among important pathogens of *P. heldreichii* growing under harsh conditions. *Cenangium accum* was recorded on *P. heldreichii* needles at the PRE site (Table 3). Previously it was demonstrated as a pathogen of weakened *P. nigra* and *P. sylvestris* in the Balkan region [50,55]. It has also been shown to be associated with Cyclaneusma needle cast in *P. sylvestris* in western Poland [39]. *Cenangium accum* develops predominantly on needles damaged by frosts or drought and its development is favoured by a high humidity. Similarly, *Dothistroma septosporum*, was recorded at a low relative abundance and only on the PRE (Table S1). Investigation of *Dothistroma septosporum* accomplished using PCR and species-specific primers, has demonstrated the presence of this potentially invasive pathogen across the *P. heldreichii* distribution range in Montenegro [17]. It was suggested that enviromental conditions present at *P. heldreichii* sites suppress the development of Dothistroma needle blight and that *P. heldreichii* is only slightly susceptible to *D. septosporum.* However, the infection level by *D. septosporum* may vary in different years, especially after rainy periods [10,17].

Among the other fungi detected, there were fungal endophytes that are known to be commonly associated with different tissues of forest trees [61,62]. The relationship between endophytic fungi and plants is not clearly understood and may change depending on the health status of the plant. However, for at least part of their life, they colonise plant tissues asymptomatically [12,63,64]. In the present study, *L. pinastri* and *L. conigenum* were the most commonly detected fungi on *P. heldreichii* needles. They are globaly distributed and commonly associated with pines [42,65]. It was shown recently that *L. pinastri* colonises healthy needles latently as an endophyte, initiates active growth at the begining of needle senescence and sporulates after needle fall. It is a dominant coloniser of dying needles and a saprotroph contributing to a needle decomposition [39]. In the Balkan region, *L. pinastri* was frequently reported from *P. nigra* and *P. sylvestris* grown in forests, forest nurseries and plantations [50,55]. *Lophodermium conigenum* is known as a coloniser of damaged needles and can also form fruitbodies on broken branches [65]. In comparison to *L. pinastri*, it is less frequently reported on both *P. nigra* and *P. sylvestris* in Serbia [50] and on *P. sylvestris* in Spain [45] and the UK [66]. In the present study, *L. pinastri* was present on all investigated sites, while *L. conigenum* showed a higher relative abundance on sites characterised by moderate growth conditions (PRO, KKO and PRE). The relative abundance of *L. pinastri* was at least 10 times lower on the KKN site than on any other site, what was likely due to its exclusion by pathogenic fungi which dominated on this site.

Some fungal endophytes of conifers can produce bioactive secondary metabolites, supporting ecological adaptations of host plants owing strong antimicrobial activities [63,67]. *Microsphaeropsis olivacea* (syn. *Coniothyrium olivaceum*) is a fungal endophyte, which was shown to produce bioactive compounds that are considered as promising antibacterial and antifungal agents [68]. *Microsphaeropsis olivacea* was among the most abundant species on the KKN, where *Sydowia polyspora* and *Neocatenulostroma germanicum* were also abundant. *Mollisia ligni* has been shown to produce mollisin, a compound that is known for a strong fungicidal activity against *Sydowia polyspora* [69]. A high frequency of *Mollisia ligni* may result in decreased frequency of *S. polyspora* on *P. sylvestris* [39]. On *P. heldreichii*, *N. germanicum* was absent in the cases where *M. ligni* was present (KMT).

The study has also detected a number of fungi previously known as pathogens of deciduous trees and/or of agricultural crops including *Phaeomoniella* 2814\_15, *Ramoconidiophora euphorbiae* and *Geastrumia sp.,* which were for the first time recorded on conifers, possibly as latent endophytes. Fungi from genus *Phaeomoniella* are important pathogens of grapevines, causing grapevine trunk disease [70], *Ramoconidiophora* and *Collophorina* species have been reported from necrotic and symptomless wood and leaves of *Prunus*, *Castanea, Vitis,* and from roots of *Calluna* [71], while *Geastrumia sp.* is known as one of the species involved in disease complex of sooty blotch and flyspeck of apple (*Malus domestica*) [72]. *Phaeomoniella sp., Ramoconidiophora, Collophorina* and *Geastrumia* were more abundant on sites with harsh growth conditions (ORJ and KMT). Among other pathogens, *Lachnellula calyciformis* was abundant on sites under harsh growing conditions, while *Phaeosphaeria pontiformis* was on sites characterised by moderate growing conditions. Corticoid fungus from genus *Athelia*, which typically occurs as saprotroph [73,74] was present on sites with moderate growth conditions.

Fungi may modulate stress tolerance, enhance growth and increase reproduction [7,47,75]. Habitats subjected to a high abiotic stress can be inhabited by a species-specific fungal community [76,77]. Our results demonstrated that needles of *P. heldreichii* constitute a habitat for a species-rich community of fungi, the composition of which was found to be largely driven by environmental conditions and/or health status of host trees. Indeed, similar patterns were also observed for fungal communities in roots of *P. heldreichii* [5], indicating that without the habitat-adapted fungal symbionts, plants are hardly capable of surviving in high stress habitats.

Fungal communities evolve together with host plants [12], suggesting that relic and endemic pine species such as *P. heldreichii* [5,54], can be associated with a specific fungal community. For example among the 50 mushroom species, which were recorded in *P. heldreichii* forests to date [78–81], many were rarely observed and exclusively recorded in *P. heldreichii* forests. Hence, *Chalciporus ammarelus*, *Geastrum minimum*, *Hygroporus gliocyclus*, *Hygroporus hypothejus*, *Rhizopogon roseolus* and *Morchella esculenta* s.l. are listed on the preliminary red list [79], and thus, are protected [82,83]. Furthermore, recently described species *Erioscyphella curvispora* [84] was discovered in *P. heldreichii* needle litter and a newly described genus *Perzia*, accommodated by a type species *Perzia triseptata*, gen. nov. was discovered on the xeric bark of *P. heldreichii* [85]. Other findings include newly described species *Velutarina bertinscensis* [86], *Peziza montrivicola* [87], *Cenangiopsis ragvanii*, *C. junipericola* [88] and rare species *Trichophaea flavobrunea* [89] that were determined on other tree species growing on *P. heldreichii* sites. Those findings represents a special value to science, and together with overall recorded fungal diversity demonstrate a high value of *P. heldreichii* habitats as being unique and biodiversity hotspots in high-altitude mountain areas.

In summary, the results demonstrated that needles of *P. heldreichii* were associated with a species-rich community of fungi, the composition of which was found to be largely dependent on environmental conditions and/or health status of host trees. Needle pathogens appeared to be an important component of fungal communities associated with needles of *P. heldreichii*, but caused limited damage likely due to interaction with other fungal species and/or due to plant defence mechanisms. In order not to become a threat to the health and growth of forest trees, the present and potentially new invasive pathogens should be monitored regularly. Further research is also needed to understand patterns of different fungal species coexistence.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/5/172/s1, Table S1: Relative abundance of fungal taxa detected in needles of *Pinus heldreichii* at six different sites in Montenegro, Table S2: Sørensen similarity index of the phyllosphere fungal communities among the six *Pinus heldreichii* sampling sites.

**Author Contributions:** Conceptualisation, J.L. and A.M.; methodology, J.L. and A.M.; formal analysis, J.L. and A.M.; resources, J.L.; data curation, J.L. and A.M.; writing—original draft preparation, J.L.; writing—review and editing, J.L. and A.M.; funding acquisition, J.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministry of Science of Montenegro through the grant INVO HERIC No: 01-646. A.M. was supported by the Swedish Research Council Formas (grant no. 2019-00597).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Phylogenetic Characterization of** *Botryosphaeria* **Strains Associated with** *Asphondylia* **Galls on Species of Lamiaceae**

#### **Beata Zimowska 1, Sylwia Oko ´n <sup>2</sup> , Andrea Becchimanzi <sup>3</sup> , Ewa Dorota Krol <sup>1</sup> and Rosario Nicoletti 3,4,\***


Received: 20 December 2019; Accepted: 20 January 2020; Published: 21 January 2020

**Abstract:** In the last decade, *Botryosphaeria dothidea* has been steadily reported as an associate of gall midges (Diptera, Cecidomyiidae) in a variety of host plants and ecological settings. This cosmopolitan fungus is well-known for its ability to colonize many plant species, as both a pathogen and an endophyte. Thus, the shift from this general habit to a lifestyle involving a strict symbiotic relationship with an insect introduces expectancy for possible strain specialization which could reflect separated phylogenetic lineages. Considering the recent taxonomic revision concerning species of *Botryosphaeria*, we evaluated the phylogenetic relationships among strains recovered from *Asphondylia* galls collected on several species of Lamiaceae in Poland and in Italy, and all the currently accepted species in this genus. A number of strains previously characterized from gall samples from Australia and South Africa, whose genetic marker sequences are deposited in GenBank, were also included in the analysis. As a result, full identity as *B. dothidea* is confirmed for our isolates, while strains from the southern hemisphere grouped separately, indicating the existence of genetic variation related to the geographic origin in the association with gall midges.

**Keywords:** *Asphondylia*; *Botryosphaeria*; *B. dothidea*; DNA sequencing; gall-associated fungi; Lamiaceae; phylogenetic relationships; symbiosis

#### **1. Introduction**

Although the nature of their symbiotic relationship has not been clearly ascertained, the occurrence of *Botryosphaeria dothidea* as an associate of many gall midge species (Diptera, Cecidomyiidae) is steadily reported, regardless of host plants and ecological contexts. For a long time, the identity of the fungal symbiont has been controversial, by reason of inherent difficulties in the isolation procedure, and of several taxonomic reassessments. In fact, several other fungi, such as *Cladosporium* spp. and *Alternaria* spp., have been frequently reported as gall associates, basically in connection with their saprophytic aptitude, which occasionally makes them conceal the real symbiont during the isolation procedure [1–4]. On the other hand, nomenclatural inconsistency, which only recently has been resolved after the epitypification of *B. dothidea* [5], may account for some previous incorrect reports referring to *Macrophoma*, *Diplodia*, and *Dothiorella* [2,3,6,7]. Indeed, taxonomists recommend a careful interpretation of past literature concerning this fungus [8,9].

Like several species in the Botryosphaeriaceae, *B. dothidea* is well-known for its cosmopolitan distribution and ability to colonize a high number of plants, either as a pathogen or as an endophyte [8–10]. The involvement in cecidomyid-gall formation on a variety of plants confirms it as a very adaptive species. Although the shift from association with plants to a lifestyle characterized by a strict symbiotic relationship with an insect suggests possible strain specialization, which could reflect separated phylogenetic lineages, data resulting from previous studies did not provide evidence for this hypothesis [1,11]. However, in the last decade, the evolution in fungal taxonomy boosted by the application of DNA sequencing has provided a remarkable contribution in view of a better resolution of the *Botryosphaeria* species aggregate. Six species presenting a *Fusicoccum* anamorph, namely *B. corticis*, *B. dothidea*, *B. fabicerciana*, *B. fusispora*, *B. ramosa*, and *B. scharifii*, were recognized based on a detailed phylogenetic analysis; another species, *B. agaves*, is included in this group, although its anamorphic stage has never been described so far [8]. In addition to this basic set, several new species have been more recently identified, mostly based on isolations from tree plants in China (Table 1).

Species of Lamiaceae are widespread in the Mediterranean region, where they seem to represent a diversity hotspot for gall midges. In fact, two new species of the genus *Asphondylia* have been recently described from galls collected on host plants such as *Coridothymus capitatus* [12] and *Clinopodium nepeta* [4], and two more are in course of characterization from *Micromeria graeca* and *Clinopodium vulgare* (Viggiani, personal communication). However, their distribution appears to reach Central Europe, following the geographical spread of some hosts, as documented in the case of *A. serpylli* and *A. hornigi*, respectively associated with *Thymus* spp. and *Origanum vulgare* [3,13]. In the course of our investigations on a number of species of Lamiaceae, we had the opportunity to collect isolates from galls of *Asphondylia* spp. in two European countries with different climatic conditions. Despite a certain variation in biometric characteristics, sequence homology of internal transcribed spacers of ribosomal DNA (rDNA-ITS) confirmed *B. dothidea* as the fungal symbiont of gall midges in both contexts. However, the relative unreliability of data available in GenBank for this species, which is to be taken into account after the recent rearrangements within *Botryosphaeria*, prompted us to undertake a more accurate study of the phylogenetic relationships with the new taxa and other isolates from *Asphondylia* galls from other plant species/countries, in order to assess whether or not strains adapted to this particular symbiotic relationship are homogeneous in taxonomic terms.


**Table 1.** *Botryosphaeria* species with *Fusicoccum* anamorph described after 2013.

#### **2. Materials and Methods**

#### *2.1. Isolation and Morphological Observations*

Isolations of fungal associates of gall midges were carried out on potato dextrose agar (PDA) amended with 85% lactic acid (1 mL·L−1) in 90 mm diameter Petri dishes. After removing the outer residues of the flower calyx, fragments from gall walls were cut and transferred onto the agar medium. Plates were incubated in darkness at 25 ◦C. Hyphal tips from the emerging fungal colonies of the botryosphaeriaceous morphotype were transferred to fresh PDA plates, for morphological observations and storage of pure cultures at 4 ◦C. Production of pycnidia was induced in cultures prepared in plates containing 2% water agar (WA), topped with sterilized pine needles, which were kept at room

temperature under near-UV illumination [21]. Observations were carried out under a Motic BA 210 microscope (Xiamen, China), and images were taken through a 1 MP Motic camera and Scopelmage 9.0 software. Minimum, maximum, mean values, and the length/width (L/W) ratios of 50 conidia from each isolate were measured for a comparison with previously annotated reference sizes for *Botryosphaeria* species [8,17].

#### *2.2. DNA Sequencing and Phylogenetic Analysis*

Selected strains recovered from galls on several species of Lamiaceae collected in Italy and Poland (Table 2) were sampled from the surface of PDA cultures with a scalpel. The mycelial matter was transferred to 1.5 mL Eppendorf tubes, for DNA extraction. DNA isolation was performed by using a DNA easy plant and fungi isolation kit (EurX, Gda ´nsk, Poland), according to manufacturer's protocol. DNA concentration was estimated on 1.5% agarose gel, compared with GeneRulerTM DNA Ladder Plus (Thermo Scientific, Waltham, MA, USA), and measured by using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). DNA samples were diluted to a concentration of 20 ng·μL−<sup>1</sup> and stored at 20 ◦C. Amplification of loci currently considered in taxonomy of *Botryosphaeria* [8] was carried out, using primers ITS1 and ITS4 for the rDNA-ITS region [22], and primers EF1-728F and EF1-986R for the translation elongation factor 1-alpha (TEF1) region [23]. PCR reaction mixtures contained 20 ng of genomic DNA, 0.2 mM of dNTP, 0.2 mM of each primer, 1 × Taq buffer mM buffer (10 mM of Tris-HCl, 1.5 mM of MgCl2, and 50 mM of KCl), and 1 U of Taq polymerase, and were adjusted to a final volume of 25 μL with sterile distilled water. PCR was conducted in a Biometra T1 thermocycler (Analytik Jena, Jena, Germany). The following reaction profile was applied: 95 ◦C—5 min, 35 cycles (95 ◦C—45 s, 52 ◦C—45 s, and 72 ◦C—45 s), with final elongation at 72 ◦C—5 min. PCR products were separated in 1.5% agarose gels containing EtBr in TBE buffer, at 140 V, for 1 h.


**Table 2.** Isolates of *B. dothidea* from *Asphondylia* galls collected on Lamiaceae used in this study.

ITS: internal transcribed spacer; TEF1: translation elongation factor 1-alpha.

After checking and determining the size of the resulting PCR products, we submitted samples to Genomed (Warsaw, Poland), for sequencing. The obtained nucleotide sequences were compared with reference strains of *Botryosphaeria* spp. from GenBank. All sequences were checked and manually edited by using CLC Main Workbench 8.1.2 software (QIAGEN, Aarhus, Denmark) where necessary. Besides our original sequences, additional sequences of *Botryosphaeria* isolates from *Asphondylia* galls were searched in GenBank for inclusion in the phylogenetic analysis, where a strain of the species *Botryobambusa fusicoccum* was used as outgroup (Table 3). The combined and single ITS and TEF1 sequences were aligned by using Muscle [24] and manually adjusted with AliView software [25], where necessary. The phylogenetic analysis was conformed to a recent protocol [26]. Congruence between the different datasets was tested by using the partition homogeneity test in PAUP software version 4.0b10 [27]. Gaps were treated as missing characters. Phylogenetic analyses of the concatenated and single-sequence data for maximum likelihood (ML) were performed by using RAxML software version 8.2.12 [28] with GTR+G model of nucleotide substitution and 1000 bootstrap replications. Concatenated sequences were also analyzed for maximum parsimony (MP) by using PAUP, under the heuristic search parameters with tree bisection reconnection branch swapping, 100 random sequence additions, maxtrees set up to 1000, and 1000 bootstrap. Posterior probabilities of the concatenated dataset were determined by Markov Chain Monte Carlo (MCMC) sampling in MrBayes version 3.0b4 [29]. MCMC chains were run for 4000,000 generations, sampling every 100, with a 25% burn-in discarded. Phylogenetic trees were drawn by using FigTree software [30]. Both the alignments and the trees of concatenated dataset were deposited in TreeBase (http://purl.org/phylo/treebase/phylows/study/TB2:S25558).


**Table 3.** Reference strains used in the phylogenetic analysis.

#### **3. Results**

Cultures on PDA of *Botryosphaeria* isolates recovered from *Asphondylia* galls on Lamiaceae displayed a sparse to moderately dense aerial mycelium, with diverse colors, from white-cream to gray to olivaceous-black, darkening with age, occasionally with narrow or wider columns of mycelium (Figure 1A). Pycnidial conidiomata developed after 10–14 days on PDA, or 8–12 days on pine needles in WA (Figure 1B). These fruiting bodies released buff, and, respectively, black (Figure 1C) or cream masses of spores containing typical *Fusicoccum* conidia, one-celled or with one septum (Figure 1D). They were smooth, hyaline, mostly with granular content, fusiform or irregularly fusiform, wider in the middle to upper third, base-truncate or subtruncate with rounded apex, and quite variable in size (Table 4). Muriform conidia referable to the synanamorphic stage *Dichomera* were never observed, unlike what previously resulted in subcultures of strains from galls collected on *T. vulgaris* directly prepared from the isolation plates [3]. This failure was assumed to possibly derive from prolonged storage at 4 ◦C of the stock cultures of our strains. Likewise, no isolate produced ascomata throughout the observation period.

**Figure 1.** (**A**) Variable morphology of isolates from *Asphondylia* galls collected in our study. (**B**,**C**) Production of pycnidia, respectively, on pine needles on WA (water agar) and on PDA (potato dextrose agar). (**D**) *Fusicoccum* conidia.


**Table 4.** Morphological features of *B. dothidea* isolates from Lamiaceae examined in this study.

Because of the above unreliability of morphological characters, DNA sequence homology was fundamental for an accurate taxonomic identification. PCR products of approximately 560 bp for ITS region and 300 bp for TEF1 region were amplified and successfully sequenced in 14 isolates from *Asphondylia* galls considered in this study. All the nucleotide sequences obtained were deposited in GenBank (Table 2), and blasted against the ex-epitype strain of *B. dothidea* (CMW8000/CBS115476) [5]. Identity was 99.79% for all strains but one (MgBt1, 98.97%) for ITS sequences, while for TEF1 sequence identity was 100% except four isolates (Th/g2017, OvFs/g, OvdF3e, and again MgBt1) matching at 99.16%.

This identity was highlighted in the subsequent phylogenetic analysis considering reference strains of all the recognized *Botryosphaeria* spp. producing *Fusicoccum* conidia, along with previously identified strains of *B. dothidea* collected from *Asphondylia* galls, in other contexts, worldwide. Although several contributions have been published in recent years on the subject of gall midges and associated fungi, a search in GenBank showed that sequences of both ITS and TEF1 are only available for some isolates from *Acacia* spp., collected in Australia and South Africa, which were the subject of a previously mentioned study [1].

The trimmed and manually adjusted alignment of concatenated locus contained 58 strains (including the outgroup) and consisted of 490 and 255 bp for ITS and TEF1, respectively. The best scoring RAxML tree (Figure 2) had a final likelihood value of −1961.235515. The matrix had 172 distinct alignment patterns, with 5.74% of undetermined characters or gaps. The ML tree of ITS alone showed poor resolution compared to ML trees based on TEF1 and concatenated sequences (Figure S1). ML trees of TEF1 and ITS + TEF1 showed almost the same topology, except for *B. qingyuanensis*, which is included in the same clade of *B. wangensis-B. sinensis-B. qinlingensis* in the phylogram based on TEF1 alone (Figure S2). Parsimony analysis yielded 1000 equally parsimonious trees (tree length = 174 steps; consistency index = 0.902; retention index = 0.915; relative consistency index = 0.826; homoplasy index = 0.098). Of the 745 characters used, 76 were parsimony-informative, 65 were variable and parsimony-uninformative, and 604 were constant. The same clades were supported in MP and ML analyses, except for *B. minutispermatia*, which resulted in being less closely related to the Australian isolates in the MP tree (Figure S3).

**Figure 2.** Phylogenetic tree based on maximum likelihood (ML) analyses of concatenated ITS and TEF1 sequences from strains considered in this study. Original isolates from Lamiaceae are in bold. Bootstrap support values ≥ 60% for ML and maximum parsimony (MP) are presented above branches as follows: ML/MP, bootstrap support values < 60% are marked with '-'. Branches in bold are supported by Bayesian analysis (posterior probability ≥ 90%). *Botryobambusa fusicoccum* MFLUCC 11-0143 was used as the outgroup reference.

#### **4. Discussion**

The wide morphological variation observed in our study is in line with recent findings that biometric data of *Botryosphaeria* species may overlap and are no longer relevant for the identification of *B. dothidea* [8,9,17]. Even *Dichomera* conidia are infrequently observed in this species [8,14,21], indicating that this character cannot be reliably taken into account for taxonomic purposes. Hence, reports from *Asphondylia* galls primarily referring to this anamorphic stage [1] should be considered with caution, because of the common co-occurrence of saprophytic *Alternaria* producing similar phaeodictyospores.

Therefore, nucleotide sequencing has become the primary identification method for *B. dothidea* through the assessment of homologies with the ex-epitype strain (CMW8000/CBS115476) [5], particularly considering the loci of ITS and TEF1. In fact, although ITS was shown to clearly distinguish *B. dothidea* from its closest relatives, the recent discovery of several cryptic species within the Botryosphaeriaceae makes it necessary to combine with TEF1 for a more reliable identification [8,31,32].

Phylogenetic analysis disclosed a clear identity with *B. dothidea* of our heterogeneous strain sample from Lamiaceae, regardless of their origin from two different climatic contexts. In fact, both Polish and Italian isolates grouped together with the type strains of this species. Conversely, Australian and South African strains formed two distinct groupings, which indicated a divergence from the most recent common ancestor. Particularly, the former was associated to type strains of the species *B. minutispermatia* in a more comprehensive clade, while the latter exhibited a slightly higher phylogenetic distance. Considering that both Australian and South African isolates were collected in association with gall midges on *Acacia* spp., this finding could be interpreted as being in agreement with a previous analysis

of the bulk *B. dothidea* sequences deposited in GenBank, showing a population structure which is shaped by geographical distance rather than host-plant preference [9]. In fact, lower identity with *B. dothidea* strains in the GenBank database resulted for these isolates, particularly those from South Africa whose TEF1 sequence homology was not higher than 96.83%. The observed variation particularly concerning this genetic marker requires further assessments to establish if these clusters should be interpreted as separated lineages within *B. dothidea*, or if they may represent distinct species.

The case of *B. auasmontanum*, a species characterized from a single strain (CBS121769) collected in Namibia [14], deserves further consideration. In fact, besides the holotype, our analysis included two more isolates from Italy ascribed to this species [26], whose ITS and TEF1 sequences are available in GenBank. Unexpectedly, they did not cluster with the holotype in the phylogenetic tree (Figure 2). A BLAST search in GenBank with sequences by these two strains clearly shows their 100% identity with *B. dothidea*, while ITS and TEF1 sequence homology with CBS121769 is lower (95.21% and 90.28%, respectively). The phylogenetic separation of the holotype of *B. auasmontanum* is supported by two large gaps resulting from the sequence alignment (Figure 3), which means that the claimed evidence of possible synonymy between the two species [26] applies only to the two Italian isolates. Thus, their inclusion in phylogenetic studies as representatives of *B. auasmontanum* should be avoided, and their incorrect taxonomic identification should be taken into account. Figure 3 also describes the substantial similarity among all strains of *B. dothidea*, including Australian and South African isolates, and the closely related *B. minutispermatia*, which essentially differ for nucleotide substitutions at single definite positions.

**Figure 3.** Alignment of ITS and TEF1 sequences of strains of *B. dothidea*, *B. minutispermatia*, and *B. auasmontanum.* Gaps in the sequences of the holotype of the latter species (CBS121769) are evident, along with single nucleotide substitutions in TEF1 sequences of Australian and South African isolates, and the type strains of *B. minutispermatia*.

#### **5. Conclusions**

Characterization of *Botryosphaeria* strains from different plant species and environmental contexts goes through the finding of novel species, with the expectation that additional undescribed taxa may be discovered in the future [17]. This concept involves strains associated with cecidomyids, inducing gall formation on a wide variety of plant species, which so far have been overlooked in assessments concerning taxonomy of such a highly adaptive and widespread fungus. Results of our investigation showed full identity with *B. dothidea* of isolates from galls collected from Lamiaceae, while a possible separation from this species should be verified for isolates previously recovered from *Acacia* in Australia and, particularly, South Africa. Indeed, a more adequate definition could be obtained by integrating these findings with data concerning cecidomyid-associated *Botryosphaeria* strains from other countries. Therefore, a more active cooperation among researchers currently working on this topic worldwide is to be encouraged, in order to shed further light on this unique biological association.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/1424-2818/12/2/41/s1. Figure S1: Maximum-likelihood tree of ITS sequence. Figure S2: Maximum-likelihood tree of EF sequence. Figure S3: One of the 1000 most parsimonious trees resulting from the analysis of concatenated ITS and TEF1 sequences.

**Author Contributions:** Conceptualization, B.Z. and R.N.; methodology, A.B., B.Z., E.D.K., and S.O.; formal analysis, A.B. and S.O.; writing—original draft preparation, B.Z., R.N., and S.O.; writing—review and editing, B.Z., E.D.K., and R.N.; funding acquisition, E.D.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Authors thank Sarah Lucchesi (University of Southern Maine, Portland, USA) for revising the English style of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
