Characterization of Boxwood Shoot Bacterial Communities and Potential Impact from Fungicide Treatments

Agrochemicals are important tools for safeguarding plants from invasive pathogens, insects, mites, and weeds. How they may affect the plant microbiome, a critical component of crop health and production, was poorly understood. ABSTRACT Phyllosphere bacterial communities play important roles in plant fitness and growth. The objective of this study was to characterize the epiphytic and endophytic bacterial communities of boxwood shoots and determine how they may respond to commonly used fungicides. In early summer and early fall, shoot samples were collected immediately before and 1, 7, and 14 days after three fungicides containing chlorothalonil and/or propiconazole were applied to the canopy. Total genomic DNA from shoot surface washings and surface-sterilized shoot tissues was used as the template for 16S rRNA metabarcoding, and the amplicons were sequenced on a Nanopore MinION sequencer to characterize the epiphytic and endophytic communities. The bacterial communities were phylogenetically more diverse on the boxwood shoot surface than in the internal tissue, although the two communities shared 12.7% of the total 1,649 identified genera. The most abundant epiphytes were Methylobacterium and Pantoea, while Stenotrophomonas and Brevundimonas were the dominant endophytes. Fungicide treatments had strong impacts on epiphytic bacterial community structure and composition. Analysis of compositions of microbiomes with bias correction (ANCOM-BC) and analysis of variance (ANOVA)-like differential expression (ALDEx2) together identified 312 and 1,362 epiphytes changed in abundance due to fungicide treatments in early summer and early fall, respectively, and over 50% of these epiphytes were negatively impacted by fungicide. The two chlorothalonil-based contact fungicides demonstrated more marked effects than the propiconazole-based systemic fungicide. These results are foundational for exploring and utilizing the full potential of the microbiome and fungicide applications and developing a systems approach to boxwood health and production. IMPORTANCE Agrochemicals are important tools for safeguarding plants from invasive pathogens, insects, mites, and weeds. How they may affect the plant microbiome, a critical component of crop health and production, was poorly understood. Here, we used boxwood, an iconic low-maintenance landscape plant, to characterize shoot epiphytic and endophytic bacterial communities and their responses to contact and systemic fungicides. This study expanded our understanding of the above-ground microbiome in ornamental plants and is foundational for utilizing the full benefits of the microbiome in concert with different fungicide chemistries to improve boxwood health. This study also sets an example for a more thorough evaluation of these and other agrochemicals for their effects on boxwood microbiomes during production and offers an expanded systems approach that could be used with other crops for enhanced integrated pest management.

B acteria are the largest microbial group in the plant phyllosphere and are indispensable to plant development. An increasing body of studies has documented that many phyllosphere bacteria are beneficial to plant health and growth (1)(2)(3)(4)(5). These bacteria are generally divided into two groups based on where they live: bacteria living on plant surfaces are considered epiphytes, while those living in internal tissues are regarded as endophytes (6,7). Some endophytes can originate from the ingress of epiphytes into the plant tissue and vice versa, forming a phyllosphere continuum (8). Both epiphyte and endophyte colonization and succession are subject to constraints imposed by the plant and climate (9,10). In crop production systems, phyllosphere bacterial communities are also influenced by agricultural practices, including fertilization (11)(12)(13), irrigation (14), tillage (15,16), and agrochemicals (17,18). Understanding how various agricultural activities influence the plant microbiome is pivotal for establishing a systems approach to crop health and production.
Boxwood is an iconic and low-maintenance landscape plant and the number one evergreen nursery crop sold in the United States (19)(20)(21). Several studies have explored boxwood leaf endophytes with the goal of improving plant growth and health. For example, Burkholderia sp. SSG, an endophyte isolated from boxwood leaves, is an effective biological control agent against diverse pathogens (22) and a bio-fertilizer promoting boxwood development (3). Community-wide, Kong et al. (23) demonstrated that the tolerance of the ordinarily highly susceptible English boxwood (Buxus sempervirens Suffruticosa) to the destructive disease boxwood blight was linked to the culturable endophytic communities. These studies reveal the potential benefits that a microbial component could bring to boxwood health and production. However, little is known about the diversity and composition of the microbial communities in the boxwood phyllosphere on a larger scale and how they may respond to different agrochemicals, such as fertilizers and pesticides used for increasing crop yields and protecting plant health.
Like all pesticides, fungicides are generally placed into two categories: systemic and contact, according to their ability to penetrate plant surface barriers and to be distributed throughout the plant through xylem or phloem. For example, propiconazole is a systemic triazole fungicide with a specific mode of action: inhibiting fungal cell demethylase activity during ergosterol biosynthesis (24). This chemistry can translocate through the plant cuticle into the xylem and reach leaf tips and margins (25). In contrast, chlorothalonil is a broad-spectrum contact fungicide that works on the plant surfaces and interacts with multiple thiol-dependent enzymes in respiration and metabolism in fungal cells (26,27). In the field, these two types of fungicides are often used in mixtures or rotated one after another to provide better crop protection and lower the risk of developing fungicide resistance.
Propiconazole and chlorothalonil also have direct or indirect actions against bacteria. Propiconazole is reported to inhibit overall bacterial activity (28) and alter soil bacterial composition for 75 days when applied at a high concentration of 100 mg/kg (29). However, propiconazole may not act directly on bacteria, as they lack sterols in their cell structure, and the observed effects of propiconazole on bacterial communities are instead considered to be associated with microbial competition (30). Chlorothalonil, on the other hand, can bind to glutathione and deactivate thiol-dependent enzymes in bacteria (31) and thus has shown a direct impact on various bacterial communities in soils (32,33), pollinator guts (34), and aquatic systems (35). In particular, Chen et al. (36) showed that soil bacteria involved in nitrogen cycling were affected by chlorothalonil applied at the label rate in laboratory experiments. Similarly, Díaz Rodríguez et al. (33) demonstrated that chlorothalonil at a concentration of 4.3 g/L inhibited indoleproducing bacteria isolated from the wheat rhizosphere, represented by the genera Bacillus, Lysinibacillus, Paenibacillus, Pseudomonas, and Enterobacter. However, only limited studies have documented the effects of these two fungicides on foliar epiphytic and endophytic bacteria. Under greenhouse conditions, Yadav et al. (37) found that propiconazole had an inhibitory effect on the growth of the culturable leaf epiphytic and endophytic bacteria in tomato, and with high-throughput sequencing, they also showed that propiconazole reduced the endophytic bacterial community diversity and altered the community composition (37). In contrast, Doherty et al. (18) found that chlorothalonil did not have an effect on the culturable bacterial population on creeping bentgrass leaves. As foliar bacteria can play important roles in promoting plant health and productivity, a better understanding of whether and how foliar epiphytic and endophytic bacteria may be affected by these fungicides at a larger community scale under field conditions is fundamental to minimizing the negative impacts of fungicides and exploiting the beneficial microorganisms for better boxwood protection and production.
In this study, we characterized boxwood shoot epiphytic and endophytic bacterial communities and examined their responses to propiconazole (trade name: Banner Maxx) and chlorothalonil (trade name: Daconil Weather Stik, here referred to as Daconil) and their mixture (Concert II), three major fungicides registered for control of boxwood blight and other fungal diseases (20). We collected boxwood shoots at several time points before and after fungicides were applied to the canopy (cover spraying) in early summer and early fall. The Nanopore MinION platform was used to sequence the 16S rRNA amplicon for identification of epiphytic and endophytic bacteria.

RESULTS
Sequencing summary. A total of 7,844,581 and 13,033,542 16S rRNA amplicon raw reads were generated for epiphyte and endophyte samples, respectively (see Tables S1 and S2 in the supplemental material). About 13.2% and 99.4% reads from the epiphyte and endophyte samples, respectively, were excluded from analyses due to matching the sequences of the boxwood chloroplast. After further data cleaning and removal of the samples that did not meet the threshold, only 6,556,458 reads were retained for epiphytes and 71,109 reads for endophytes. There were differences in the epiphyte and endophyte sequence reads among the samples. For epiphytes, the Concert II-treated and nontreated control samples had the highest reads/sample, while those treated with Banner Maxx and Daconil had the lowest in early summer and early fall, respectively (Table S1). For endophytes, Banner Maxx-impacted samples had the least high-quality sequence reads in early summer but the most in early fall. Comparatively, Daconil-impacted samples had the most, while nontreated controls had the least in early summer and early fall (Table S2).
There also were differences in sequencing depth as measured by Good's coverage scores between epiphyte and endophyte samples. The coverage scores for epiphytes ranged from 0.975 to 0.995. However, that score was lower for the endophytes, ranging from 0.818 to 0.924. Rarefaction curves showed that sequencing reached near plateau for the epiphyte samples, indicating a sufficient depth for characterizing these communities and assessing how they were affected by fungicide treatments. However, sequencing depth was shallow for the endophyte samples (Fig. S1).
Overall epiphytic and endophytic community compositions. Boxwood shoot epiphytic and endophytic bacterial communities differed in the membership and relative abundance of individual members. A total of 1,646 bacterial genera from 39 phyla were identified from the epiphyte communities, but only 213 bacterial genera from 22 phyla were identified from the endophyte communities.
There were major differences at the phylum level between the epiphytic and endophytic bacterial communities. Although Proteobacteria was by far the dominant epiphyte and endophyte, accounting for 81.7% and 78.8% of the respective communities (Fig. 1A), there were significant differences in the relative abundance of other major phyla. Specifically, Firmicutes, Bacteroidetes, and Cyanobacteria were more abundant in the shoot internal tissue than on the surface. In contrast, Actinobacteria were more abundant on the shoot surface than in the internal tissue (Fig. 1A).
Similar differences were observed at the genus level between the epiphytic and endophytic communities. Although epiphytic and endophytic communities shared 210 bacterial genera (Fig. 1B), the former were much more diverse than the latter (1,646 versus 213 genera). The vast majority (1,436 of the 1,646) of the epiphytic bacteria were not detected in the endophytic communities. In addition, the predominant genera differed between the shoot surface and internal tissue (Fig. 1C). In particular, the most abundant epiphytic genera were Methylobacterium, Pantoea, Brevundimonas, Massilia, Sphingomonas, Lichenihabitans, Curtobacterium, Bacillus, Duganella, and Caballeronia (Fig. 1C). In contrast, the most abundant endophytic bacteria were Stenotrophomonas, Brevundimonas, Myroides, Pseudomonas, Bacillus, Staphylococcus, Gloeothece, Nostoc, Burkholderia, and Nitrosomonas (Fig. 1D). Notably, the relative abundances of these predominant genera varied between the two seasons.
Fungicide impact on bacterial community structures. Greater and more consistent differences were observed in the epiphyte compared to the endophyte community structure among fungicide treatments and among sampling times. Specifically, distinct epiphyte community structures were present among fungicide treatments and among sampling times in both early summer and early fall (Table 1 and Fig. 2). Although there were also differences in the endophyte community structure in early fall (Table 1), these differences were limited in extent-not resulting in any clear cluster among fungicide treatments or among sampling times (Fig. S2). Stronger fungicide impacts on the epiphyte community structure were observed in early fall than early summer. For example, fungicidal effect accounted for greater community structure variation in early fall than early summer (0.254 versus 0.086; Table 1). This difference is also visually evidenced by the boxwood shoot samples from different fungicide treatments forming distinct clusters in early fall but not in early summer ( Fig. 2A and B). Specifically, boxwood shoot samples treated with Daconil and Concert II appeared clustered together, and both were distant from those treated with Banner Maxx and nontreated controls (Fig. 2B). Additionally, the Banner Maxx-impacted samples were also distant from the nontreated controls.
Greater sampling time impacts on the epiphytic community structure, however, were observed in early summer than early fall. For instance, sampling time explained more community structure in early summer than early fall (0.445 versus 0.126). This is further evidenced by four clear clusters corresponding to the four sampling times in early summer but not in early fall ( Fig. 2C and D). Specifically, the early summer samples collected after fungicide treatment consistently differed from the pretreatment samples (boxplot in Fig. 2C); the samples collected 1 day after treatment also differed from those collected 7 and 14 days after (boxplot in Fig. 2C). Comparatively, none of the early fall samples collected from three posttreatment sampling times formed a clear cluster; nor did the pretreatment samples, indicative of a similar epiphytic community structure among the four sampling times.
Scope and nature of fungicide impact on bacterial communities. Both analysis of compositions of microbiomes with bias correction (ANCOM-BC) and analysis of variance (ANOVA)-like differential expression (ALDEx2) were used to identify bacterial genera with significant (relative) abundance changes at the community level, with ANCOM-BC being more sensitive than ALDEx2. Specifically, ANCOM-BC identified 178 epiphytes with differential abundance in early summer and 704 in early fall, while ALDEx2 identified only 6 and 40 in each season, respectively ( Fig. 3 and Fig. S3). However, only two of the endophytes identified by ANCOM-BC in each season differed in abundance due to fungicide treatments, while ALDEx2 showed none (data not shown).
Fungicide treatments affected a large number of epiphyte genera but only a few endophyte genera, as identified by ANCOM-BC. Among the 178 epiphytic genera identified in early summer, 28 were affected by all 3 fungicides, and 70 were impacted by 2 of the 3 fungicides, with the rest affected by only 1 fungicide (Fig. 3A). Likewise, among the 704 epiphytic genera identified in early fall, 183 were affected by all 3 fungicides, 238 were impacted by 2 of the 3 fungicides, with the rest affected by only 1 fungicide (Fig. 3B). Overall, Daconil affected the most bacterial genera, with 142 in early summer and 573 in early fall, followed by Concert II (125 and 525 genera), and Banner Maxx (37 and 210). On the other hand, Banner Maxx and Concert II affected two endophyte genera in early summer and early fall, respectively (data not shown).
Fungicide treatments suppressed more epiphytes than they enriched. Specifically, ANCOM-BC showed that fungicide treatments suppressed 233 and 1,087 genera in early summer and early fall, respectively, while they enriched only 71 and 221 genera during the same seasons (Table 2). Likewise, ALDEx2 showed that fungicide treatments   (Table 2). Again, these analyses showed that Daconil had an impact on the most bacterial genera, followed by Concert II and then Banner Maxx. Fungicidal effects on predominant epiphytes and endophytes. In addition to ANCOM-BC and ALDEx2, count regression for correlated observations with the beta-binomial (corncob) was used for differential abundance analysis of the predominant epiphytic and endophytic bacterial genera. For each genus, three fungicide treatments were individually compared to the nontreated controls. Likewise, three posttreatments were individually compared to the pretreatment.
Fungicide treatments significantly affected the abundance of several predominant epiphyte genera, and the affected genera differed greatly between early summer and early fall. In early summer, Pantoea and Curtobacterium were enriched on the surface of boxwood shoots treated with Daconil and/or Concert II as identified by corncob (P # 0.0010, Fig. 4A and C), ANCOM-BC (P # 0.0028, Table S3), and ALDEx2 (P # 0.0325, Table S4). Caballeronia was also enriched by Daconil and Concert II (P # 0.0010, Fig. 4A and C) or by Daconil alone (P = 0.0383, Table S3). In contrast, Methylobacterium was suppressed by Daconil and Concert II (P # 0.0010, Fig. 4A and C), while Bacillus was suppressed by Daconil alone (P = 0.0457, Table S3). In early fall, Brevundimonas and Lichenihabitans were enriched on the surface of boxwood shoots treated with Concert II as identified by corncob (P # 0.0010, Fig. 4A and C), ANCOM-BC (P # 0.0437; Table S3) and ALDEx2 (P # 0.0050, Table S4). These two bacterial genera were also enriched by Daconil (P # 0.0100, Fig. 4A and C). Additionally, Caballeronia was promoted by Concert II (P # 0.0100, Fig. 4C; P = 0.0479, Table S4) and by Daconil (P # 0.0100, Fig. 4A and C).   In contrast, both Sphingomonas and Bacillus were suppressed by Daconil (P # 0.0464, Table S3). Although fungicide treatments affected the (relative) abundance of predominant endophytes, the impact was limited in the number of genera and in extent, and the affected genera also differed between two seasons as identified by corncob. In early summer, Stenotrophomonas was suppressed by both Daconil and Concert II (P # 0.0100, Fig. 4B and D), while Brevundimonas was enriched by both Banner Maxx (P , 0.0500) and Daconil (P # 0.0100, Fig. 4B and D). In early fall, Stenotrophomonas, Pseudomonas, and Staphylococcus all were enriched by Banner Maxx (P , 0.0500. Figure 4B and D).
According to corncob, posttreatment abundance changes were observed in more epiphyte genera at greater rates in early summer than early fall. First, among the 10 predominant epiphytes, 9 changed their abundance in early summer, while only 6 changed in early fall ( Fig. 4A and C). Second, the number of significant changes (P , 0.0500) observed was 19 in early summer but only 9 in early fall. Third, 16 of the 19 observed changes were at a P value of #0.0010 in early summer, while only two of the nine changes were at the same level. Similar trends were identified by ANCOM-BC (Table S5) and ALDEx2 (Table S6).
Four major posttreatment abundance change patterns were observed in early summer. First, abundance continued to increase with increasing posttreatment time. Massilia, Sphingomonas, and Duganella seemed to have followed this pattern (Fig. 4C). Second, abundance continued to decrease with increasing posttreatment time as represented by Bacillus. Third, abundance fell at 1 day after fungicide treatment but bounced back at 7 days after and continued through 14 days after. Methylobacterium followed this pattern. Fourth, abundance jumped at 1 day after fungicide treatment but substantially dropped at 7 days after and continued through 14 days. Pantoea clearly followed this pattern in early summer and to a lesser extent in early fall. Similar trends were detected by ANCOM-BC (Table S5) and ALDEx2 (Table S6).

DISCUSSION
This study demonstrated that boxwood bacterial communities were phylogenetically very diverse and abundant, especially on the shoot surface, and they were profoundly impacted by fungicide treatments. These findings are fundamental to understanding the biology of boxwood, an iconic landscape plant in American gardens and the largest evergreen shrub crop in the United States, and to develop a systems approach to boxwood crop health and production.
The great bacterial diversity unveiled in the present study helps to understand boxwood as a low-maintenance crop. Biodiversity is a major determinant of ecosystem productivity and stability (38)(39)(40). Over 1,600 bacterial genera from 39 phyla were identified from boxwood shoot surface and internal tissue in this study. Some of these bacteria are commonly found in the aerial parts of other agricultural crops and trees (41)(42)(43)(44). This agreement suggests the ubiquitous distribution of these bacteria in air or soil (41) and their specialized mechanisms for adaptation and colonization on plant surfaces or inside the tissue (43,45). Comparatively, the bacterial communities were much more diverse on boxwood shoot surfaces than in internal tissue. Likewise, the predominant bacterial genera identified from the shoot surface and internal tissue also differed greatly. Specifically, Methylobacterium, Pantoea, Brevundimonas, Massilia, and Sphingomonas dominated the boxwood surface, while Stenotrophomonas, Brevundimonas, Myroides, Pseudomonas, Bacillus, Burkholderia, and Staphylococcus prevailed in the internal tissue (Fig. 1). These differences were not unexpected, as endophytic colonization is highly regulated by the plant host (46). DAA, each of the three fungicide treatments was compared to the nontreated controls, while each of the three posttreatment sampling times was compared to the pretreatment. The darkness of cell colors is scaled to model the coefficient value, with darker red indicating a higher relative abundance, while darker blue indicates a lower relative abundance compared to the control. The triangles indicate a significant taxon after FDR correction in a global test, and the significant covariate estimate (i.e., treatment or sampling time) is marked by asterisks: ***, P # 0.0010; **, P # 0.0100; *, P , 0.0500.
Similarly, several culturable endophytes, such as Pseudomonas, Bacillus, Staphylococcus, and Burkholderia, were also identified in English boxwood leaves (B. sempervirens Suffruticosa) (23). However, other dominant endophytes in the B. sempervirens cultivar Vardar Valley studied here were not found in the English boxwood. This could be due to physiological or morphological differences in the cultivars used in the two studies, as has been documented in other plant species (47), or to differences in methodology. The previous boxwood study used bacterial cultures of plant tissue and a genomic sequencing approach (23), while this study was molecularly based, with amplicon sequencing. As seen with the recently isolated boxwood endophyte Burkholderia sp. SSG that shows effective biocontrol (22) and biofertilizer activities (3), the newly identified endophytes expand the bacterial candidates that may be mined and evaluated for beneficial activities. For example, strains of Stenotrophomonas and Brevundimonas as endophytes have been shown to synthesize indole-3-acetic acid and solubilize phosphate (48)(49)(50). Our findings warrant further investigations of these endophytes to explore the potential for their use in boxwood improvement and disease management.
Several lines of evidence show that the three fungicides had broad and strong impacts on boxwood shoot epiphytes, while the effect was minimal on endophytes. First, fungicide treatments affected the abundance of a much larger variety of epiphytic than endophytic bacterial genera, according to both ANCOM-BC and ALDEx2 (882 versus 4). Second, they had a larger effect on the boxwood shoot surface bacterial community structure than on the community structure in the internal tissue. These differences were expected for the contact fungicides such as Daconil, as they work on the plant surface and do not penetrate into internal tissue and translocate internally. Even with the systemic fungicide Banner Maxx, after being sprayed onto plant foliage, it first acted on the surface bacteria until it penetrated into internal tissue; this effect was dependent upon the length of time required to penetrate and the percentage of material that entered the tissue. However, the limited effect of the systemic fungicide on boxwood endophytes is different from that reported by Yadav et al. (37). Under greenhouse conditions, they showed that 0.1% propiconazole had a suppressive effect on the growth of culturable bacterial endophytes isolated from tomato leaves. With 16S rRNA sequencing, they further showed an altering effect of propiconazole on the endophytic community structure and composition, and the fungicide enriched the relative abundance of Bacillus while it suppressed that of the genus Pseudomonas. The concentration of fungicide applied may have contributed to the different responses of endophytes to propiconazole in the two studies. We applied propiconazole at a rather low concentration (0.014%), following the label rate for effective boxwood blight management (51). Additionally, the amount of propiconazole deposited on the boxwood surface under field conditions may be less than the amount applied as described by Bai et al. (52) on wheat straw and leaves, resulting in an even lower concentration inside the boxwood tissue. This speculation is consistent with another previous study finding that the propiconazole impact on endophytic bacterial communities is chemical concentration dependent (29). The lack of fungicidal effect on boxwood shoot endophytes could also be attributed to the shallow sequencing depth due to large plant host contamination, particularly for the early summer samples. Concert II acted more like Daconil than Banner Maxx in this study; this was expected because it contains 70% chlorothalonil, the active ingredient (a.i.) in Daconil, while only 20% propiconazole, the a.i. in Banner Maxx.
The observation that Daconil and Concert II showed stronger impacts on the boxwood shoot epiphytic bacterial community than Banner Maxx was not unexpected. Chlorothalonil reacts with multiple thiol-dependent proteins, thus eventually building up toxicity in bacterial cells (26,31). In fact, chlorothalonil impacts on bacteria have been widely documented in soil, rhizosphere, pollinators, and amphibian skin (33)(34)(35)(53)(54)(55). In this study, the abundance of Pantoea, Brevundimonas, and Lichenihabitans was consistently elevated by Daconil and Concert II, indicating that Gram-negative bacteria may have some tolerance to chlorothalonil (56). In contrast, the mechanisms by which propiconazole might affect bacteria are unclear because bacteria lack the fungal cell wall structure for the fungicide to act on as an ergosterol biosynthesis inhibitor (57). However, propiconazole used at a higher concentration is known to adversely impact bacterial communities in plants (37) or soil (29). Our results indicate, however, that propiconazole used at the label rate has relatively little impact on the surface-colonizing bacteria, although this rate is effective for controlling the boxwood blight pathogen (51). This is in line with Mmbaga et al. (58), who noted that propiconazole fungicide effectively controlled powdery mildew in flowering dogwood without significantly impacting the culturable bacterial population on leaves.
We found that over 99% of the initial sequence reads for endophytes with our primer set (59,60) were from boxwood chloroplast instead of bacteria; this was rather disappointing, but important. This observation highlights the need to design better primers that will exclude chloroplast contamination from the plant host while amplifying the full length of 16S rRNA to use the full potential of the Nanopore long-read sequencing technology. It will be equally important to carefully reexamine the published sequencing data that were collected using the current primer set (59,60) and the Nanopore sequencing platform to make sure that they do not have the same problem, which may have led to erroneous conclusions in the literature.
In summary, this study advanced our understanding of the phyllosphere bacterial communities of this iconic landscape plant and important evergreen crop in the United States in a number of ways. First, boxwood phyllosphere bacteria were shown to be very diverse and abundant, especially on the shoot surface. Second, a commonly used contact fungicide chemistry, chlorothalonil, had a broader and stronger impact on boxwood shoot epiphytic bacterial communities than a commonly used systemic fungicide, propiconazole. Third, both chemistries had a limited effect on the endophytic bacterial communities. These results are the first critical step to exploring and utilizing the full potential of the phyllosphere microbiome and fungicide applications and develop a systems approach to boxwood health and production.

MATERIALS AND METHODS
Research site, boxwood plants, and site management. Boxwood shoot samples were collected from a field planting of 5-year-old Buxus sempervirens Vardar Valley at a nursery located in Ronda, North Carolina. The planting was initially established in 2017, and at the time of study, these plants stood at 60 cm tall and 40 to 50 cm wide. Fertilizer was applied early April, and herbicides Roundup (active compound: glyphosate, Bayer AG, Leverkusen, Germany) and Goal (active compound: oxyfluorfen, Nutrichem Co., Ltd., Beijing, China) were applied at label rates in late May for weed management. Pruning was performed once a year during winter.
Fungicide, application rate, and schedule. The three fungicides included in this study were Daconil Weather Stik (DL), Banner Maxx (BM), and Concert II (C2), and they were supplied by Syngenta (Greensboro, NC, USA). Each was applied at its label rate (Table 3), using a new 1-gallon handheldsprayer (HDX, Home Depot, USA), onto the boxwood canopy until runoff. Also included was a nontreated control. The four treatments, each having four replicate plants, were arranged in a randomized complete block design. As boxwood new growth starts in April, the initial treatment began on 12 April 2021. Thereafter, Daconil was reapplied every 2 weeks, while Banner Maxx and Concert II were applied every 3weeks until the first week of November. Each plant received the same treatment throughout the spraying period.
Boxwood sample collection schedule and protocol. In this study, boxwood shoot samplings were timed to take place when all three fungicides were applied on the same day, once in early summer and once in early fall. The early summer sampling started on May 26, and boxwood shoots were collected immediately before (0), and 1 day (1), 7 days (7), and 14 days (14) after fungicide application to examine the dynamics of the shoot microbiome as affected by fungicide. For each sampling, 10 to 15 7-cm-long boxwood shoots (about 20 leaves per shoot) were collected from the top and middle sections of each plant using a hand pruner (FELCO, Seattle, WA, USA), and placed in a new Ziploc bag. The pruner was sterilized with 70% alcohol and air dried between replicate plants. Sample bags were then placed in an ice chest with icepacks on the bottom. The 0-and 1-day samples were then transported to the lab at Virginia Tech's Hampton Roads Agricultural Research and Extension Center in Virginia Beach, VA, for processing the next day. The 7-and 14-day samples were collected in the same manner and shipped overnight to the lab in a cardboard box with icepacks. The early fall sampling started on 25 August 2021. All boxwood plants sampled were generally healthy at the times of sampling in both summer and fall; collected shoots samples did not have any disease symptoms. Sample processing and DNA extraction. (i) Sample preprocessing for epiphyte and endophyte collection. Epiphytes were collected from shoot samples by washing. Five symptomless boxwood shoots were arbitrarily selected from each replicate sample and then placed in a 50-mL conical tube (VWR, Radnor, PA, USA) prefilled with 45 mL sterile isotonic solution (0.15 M sodium chloride and 0.01% Tween 20). The tubes were then agitated for 30 min at 250 rpm on a rotary shaker (New Brunswick Scientific G24 environmental incubator shaker, Edison, NJ, USA) at room temperature. Following agitation, an ultrasonic cleaner (Branson Ultrasonics CPXH, Danbury, CT, USA) was used for 10 min to further dislodge surface microorganisms (61,62). The sonicated isotonic solution was retained for epiphyte analysis.
For endophytes, the five previously washed shoots were further processed to remove the remaining epiphytes following the protocol of Milazzo et al. (63) with minor modifications. Briefly, in a clean laminar flow hood, five shoots for each replicate sample were submerged and agitated for 30 s in each of the following solutions in sequence: 150 mL of 70% (vol/vol) ethanol, 150 mL of 5% (wt/vol) sodium hypochlorite, 150 mL of 70% (vol/vol) ethanol, and 300 mL of sterile distilled water. Each of the five shoots was individually dip-agitated in a series of three 50-mL conical tubes filled with sterile distilled water to wash off any remaining chemical residue or epiphytes for each replicate sample; each shoot was washed individually instead of five together. Then, the five cleaned shoots from each replicate sample were placed in a new sterile 50-mL conical tube and stored at 220°C for endophyte analysis.
(ii) Shoot washing concentration and epiphyte DNA extraction. To collect epiphytes, shoot washings were first centrifuged at 7,000 Â g for 10 min (Hermle Labnet Z-383K, Edison, NJ, USA). After the supernatant was carefully discarded, the pellet was collected and transferred to a clean 2-mL microcentrifuge tube. The pellet was further centrifuged at 15,000 Â g for 2 min (Eppendorf 5424, Hamburg, Germany) to remove any supernatant and then processed using the Qiagen DNeasy PowerSoil Pro kit (Qiagen, Hilden, Germany) to extract DNA according to the manufacturer's protocol with a minor modification. Briefly two 1-min homogenizations with a 30 s break in between were performed using a FastPrep-24 bead beating and lysis system (MP Biomedicals, Santa Ana, CA, USA) at the default speed of 4 m/s. DNA was eluted from the MB spin column membrane using nuclease-free water and then stored at 220°C until use for epiphyte analysis.
(iii) DNA extraction from washed boxwood shoots. Shoot DNA extraction was performed using sterilized utensils in a laminar flow hood. The five surface-sterilized shoots were first homogenized using mortars and pestles with liquid nitrogen. About 200 mg of the ground tissue was transferred to a zirconium bead tube (500 mm garnet and 6 mm zirconium, PFMM 500-100-25U) (OPS Diagnostics, Lebanon, NJ, USA) that was prefilled with 400 mL of AP1 buffer (provided with the Qiagen plant minikit). RNase A (4 mL; also provided with the Qiagen Plant minikit) was added to the tube, and then it was vortexed for 3 s. Tissues were homogenized further with an MP FastPrep-24 instrument at the speed of 4 m/s for 1 min. The rest of the steps followed the Qiagen DNeasy plant minikit protocols. DNA was eluted from the MB spin column membrane using nuclease-free water and stored at 220°C until it was used for endophyte analysis.
PCR amplification. The primer pair for amplifying the full length of the 16S rRNA gene was 27F and 1492R (59,60). Nanopore ligation adapters were added to the primers for downstream barcoding (Nanopore Technology, Oxford, UK). DNA extracts were diluted to 5 ng/mL, except that the epiphyte DNA samples from early summer were diluted to 0.5 ng/mL due to a low concentration of eluted DNA from extraction. The PCR mix included 0.2 mL of polymerase (TaKaRa Bio, Japan), 1 mL of primers, 5 mL 10Â buffer, 4 mL deoxynucleoside triphosphate (dNTP), and 10 mL and 1 mL DNA for early summer and early fall samples, respectively. PCR amplification was programmed as follows: 95°C for 2 min, followed by 29 cycles of 95°C for 30 s, 66°C for 45 s, and 72°C for 1 min and then one extension at 72°C for 10 min. PCR products were then purified using Promega Wizard SV gel and PCR cleanup system (Madison, WI, USA).
Nanopore library preparation and sequencing. The Nanopore SQK-LSK110 kit and EXP-PBC096 kit (Oxford Nanopore Technologies [ONT], Oxford, UK) were used to prepare the libraries. Briefly, 50 ng of each PCR product was diluted to a total volume of 24 mL for barcoding. Each reaction mix included 1 mL of ONT barcode mix, 24 mL of PCR product, and 25 mL of LongAmp Tag 2Â master mix (New England Biolabs, Ipswich, MA, USA). The thermal conditions for barcoding were programmed according to the SQK-LSK110 kit instructions. The barcoded samples were further cleaned up using AMPure XP beads (Beckman Coulter, Brea, CA, USA). In each sequencing run, equal molar amounts of 16 barcoded PCR products were pooled together to a total of 1 mg. DNA repair and end-prep were done using the NEBNext FFPE DNA repair mix and the Ultra II end-prep module (New England Biolabs) according to the SQK-LSK110 kit protocol. The NEBNext Quick T4 DNA ligase module (New England Biolabs) was used for library ligation, and the short fragment buffer (SFB) was selected for washing the library during ligation cleanup with AMPure XP beads. Then 40 fmol (34.6 ng of length 1,400 bp) of DNA library was loaded to a MinION R9.4 flow cell for sequencing following the priming and loading protocol. At the end of each run, the flow cell was washed using an EXP-WSH004 kit according to the manufacturer's protocol. The flow cell was then stored at 4°C until the next sequencing.
Bioinformatic and data processing. (i) Basecalling. The Nanopore proprietary software MinKNOW (core versions 4.3.12 and 4.4.3) and Guppy (GPU version 5.0.11) were used for base calling. The parameters for quality filtering were set at Q9 and for a read length between 1,000 bp and 2,000 bp. These base-called reads were demultiplexed by their barcodes and output in FASTQ format. The FASTQ reads were archived at the European Nucleotide Archive under study ID PRJEB56603. An in-house python package, NanoPrep (version 0.19.1, https://github.com/xpli2020/NanoPrep), was developed to facilitate further read processing, grouping, and format conversion. Briefly, NanoFilt (version 2.7.1) (64) was used to retain reads with quality of no less than 10 and remove barcode sequences by setting parameter headcrop and tailcrop to 50. These reads were then grouped by sample names and converted to FASTA format using seqtk (version 1.3) (https://github.com/lh3/seqtk) for bioinformatic analyses (see below).
(ii) Chimera sequence removal. The SILVA version 138.1 SSU NR99 (65, 66) database was used for chimera correction, as it contains robust reference sequences for 16S small subunit rRNA (67). This is less demanding on computational resources compared to using the whole-genome database compiled using Centrifuge (68) for taxonomic assignment. Because a small portion of the reference sequences still contain ambiguous bases and homopolymers (67), further processing was performed using the Quantitative Insights into Microbial Ecology (QIIME 2) (69) environment coupled with the REference Sequence annotation and CuRaIon Pipeline (RESCRIPt) (70). Briefly, the ambiguous bases (default minimal 5) and homopolymers (default minimal 8) in the SILVA database were first removed. The length of the database sequences was then filtered to retain a minimum of 1,200 bp for the SILVA database. Dereplication was followed to keep the unique sequences with different taxonomies (set p-mode "uniq"). Lastly, we trained a taxonomy classifier with the RESCRIPt evaluate-fit-classifier plugin (a wrap program to the QIIME 2 naive Bayes method). The trained classifier was then exported into FASTA format for chimera removal using Minimap2 (71) and Yacrd (72), following Cuscó et al. (73) with an in-house pipeline.
(iii) Database construction and taxonomy assignment. The software Centrifuge (68) was used for database construction and taxonomy assignment as suggested by Santos et al. (74). A reference database was compiled based on the completed bacterial whole-genome reference sequences from the National Center for Biotechnology Information (NCBI) (retrieved on 21 December 2022). There are a few advantages to using this database: (i) its taxonomy is annotated down to the species and strain levels; (ii) it is updated daily; (iii) Centrifuge supports the retrieval of the reference sequences and compilation of this database without extra code-and data-wrangling. Boxwood chloroplast sequences were also downloaded from the NCBI using the python script ncbi2kraken.py (75) searching with "(Buxus sempervirens[Organisms] OR Buxus microphylla[Organisms]) AND chloroplast[filter]". Thereafter, a local reference database consisting of a total of 30,810 bacterial and 65 boxwood chloroplast sequences was compiled and indexed with Centrifuge using the Burrows-Wheeler transform (76) and the Ferragina-Manzini indexing schemes (77). Centrifuge then classified the sequence reads and estimated their abundance using an expectation-maximization algorithm as described by Trapnell et al. (78) and Patro et al. (79). The software kraken-biom (80) was used to form the sequence abundance table for downstream analyses in R (version 4.2.2) (81). Notably, we used the term "operational taxonomic unit" ("OTU") to describe these classified full-length 16S rRNA amplicon sequences in this study, as opposed to the conventional definition based on sequence similarity clustering (82,83).
(iv) Data cleaning. First, sequence reads matching those of boxwood chloroplasts and OTUs with fewer than five sequences across all samples were removed. Second, samples with fewer than 1,000 reads for epiphytes and 100 for endophytes were excluded. Sequencing depth was evaluated with Good's coverage using the phyloseq_coverage function from the metagMisc package (84) and verified with rarefaction curves using the amp_rarecurve function from the ampvis2 package (85).
Statistical analyses. (i) Epiphytic and endophytic bacterial community structure variations. The Bray-Curtis (BC) dissimilarity (86) was used to quantify community structure variations between boxwood shoot samples. Epiphyte and endophyte abundance data were first transformed using the Hellinger method (87). Ordination was then performed on the BC dissimilarity using principal-coordinate analysis (PCoA). For statistical analysis of fungicide, sampling time, and their interactions' effect on the BC dissimilarity, the OTU abundance table was split into early summer and early fall, and these data sets were evaluated using the permutational multivariate analysis of variance method implemented by the adonis function from the vegan package with 1,000 permutations (88). Group dispersion homogeneity was also tested using the betadisper function of the vegan package with 1,000 permutations. Analysis of variance was used to evaluate the respective effects of fungicide and sampling time on the first two dimensions of the PCoA and was followed by Tukey's honestly significant difference (HSD) pairwise comparison test to elucidate whether the Bray-Curtis distance of one group was significantly different from that of the other (i.e., Banner Maxx versus nontreated, day 1 versus pretreatment). The significance level of all analyses was 0.05 unless otherwise stated.
(ii) Differential abundance analysis of fungicide impacts on bacterial communities. As suggested by Nearing et al. (89), analysis of compositions of microbiomes with bias corrections (ANCOM-BC) (90) and ANOVA-like differential expression (ALDEx2) (91)(92)(93) were used to assess whether and how fungicide treatments may have impacted the (relative) abundance of epiphytic and endophytic bacterial communities in each season. The ggvenn package (version 0.1.9) (94) was used to summarize the differentially abundant genera identified by ANCOM-BC and ALDEx2 in a Venn diagram. The changes in relative abundance of the top 10 predominant genera were also analyzed with count regression for correlated observations with the beta-binomial (corncob) (95). Note that ANCOM-BC evaluates absolute abundance, while ALDEx2 and corncob focus on relative abundance.
(a) ANCOM-BC The nonrarefied OTU table was used with the ancombc2 function to analyze the effect of fungicide and sampling time on the microbial communities at the genus level. Fungicide and sampling time were implemented as two covariates in the model formula (i.e., ; fungicide 1 sampling time), and the resulting P values were adjusted using the Benjamini-Hochberg (BH) false-discovery rate (FDR) correction (96). Default settings were used for other parameters.
(b) ALDEx2 The nonrarefied OTU table was first agglomerated to the genus level for each season. The same model formula in ANCOM-BC was used for the function aldex.clr with 1,000 Monte Carlo samplings. The denominator for the geometric mean abundance was calculated based on all the genera. A glm model was fit with the aldex.glm function, and the effect size was calculated with the aldex.glm.effect function. ALDEx2 currently uses the more conservative Holm-Bonferroni method (97) for the family-wide error rate correction.
(c) corncob The differentialTest function of the corncob package was implemented to test the fungicide and sampling time impacts on the top 10 predominant genera in each season. To test only differential abundance, overdispersion or abundance variability was controlled by applying the same formula in ANCOM-BC2 and ALDEx2 to the nonnull abundance model and both null and nonnull overdispersion models simultaneously. The Wald test was used for significance testing with 1,000 bootstraps, and the P values were corrected using the BH-FDR approach. A heatmap was made to show how significant taxa and their abundance change with respect to sampling time and fungicide application.
Data availability. The FASTQ reads are available at the European Nucleotide Archive under study ID PRJEB56603. R codes can be provided upon request.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.7 MB.