Influence of seasonal changes and salinity on spinach phyllosphere bacterial functional assemblage

The phyllosphere is the aerial part of plants that is exposed to different environmental conditions and is also known to harbor a wide variety of bacteria including both plant and human pathogens. However, studies on phyllosphere bacterial communities have focused on bacterial composition at different stages of plant growth without correlating their functional capabilities to bacterial communities. In this study, we examined the seasonal effects and temporal variabilities driving bacterial community composition and function in spinach phyllosphere due to increasing salinity and season and estimated the functional capacity of bacterial community16S V4 rRNA gene profiles by indirectly inferring the abundance of functional genes based on metagenomics inference tool Piphillin. The experimental design involved three sets of spinach (Spinacia oleracea L., cv. Racoon) grown with saline water during different seasons. Total bacteria DNA from leaf surfaces were sequenced using MiSeq® Illumina platform. About 66.35% of bacteria detected in the phyllosphere were dominated by four phyla- Proteobacteria, Firmicutes, Bacteroidetes, and Actinobacteria. Permutational analysis of variance (PERMANOVA) showed that phyllosphere microbiomes were significantly (P < 0.003) affected by season, but not salinity (P = 0.501). The most abundant inferred functional pathways in leaf samples were the amino acids biosynthesis, ABC transporters, ribosome, aminoacyl-tRNA biosynthesis, two-component system, carbon metabolism, purine metabolism, and pyrimidine metabolism. The photosynthesis antenna proteins pathway was significantly enriched in June leaf samples, when compared to March and May. Several genes related to toxin co-regulated pilus biosynthesis proteins were also significantly enriched in June leaf samples, when compared to March and May leaf samples. Therefore, planting and harvesting times must be considered during leafy green production due to the influence of seasons in growth and proliferation of phyllosphere microbial communities.


Introduction
The phyllosphere is the aerial part of plants and the most prevalent bacterial habitats on plants, as well as the primary center of photosynthetic activities. For leafy greens like spinach, it is the edible portion of the herbaceous plant with direct impact on human health due to the presence of potential human pathogens [1]. Many culture-independent studies have shown that the phyllosphere harbors many hundreds of unique bacterial taxa with bacterial community composition varying across plant species and more diverse than previously thought using culturedependent methods [1][2][3][4]. It has been shown that phyllosphere microbiomes are related to specific processes at the interface between plants, microorganisms and the atmosphere [5]. These authors noted that the phyllosphere is a major bacterial habitat on earth and these microbes are exposed to extreme, stressful, and changing environments. Also, the phyllosphere microbiota and the atmosphere present a dynamic continuum of interactions with volatile organic compounds and atmospheric trace gasses. Due to these interactions, bacterial communities on the phyllosphere are subjected to large changes in temperature, soil water content, UV radiation, relative humidity, and leaf wetness [6][7][8]. These changes combined with changes in soil physiochemical factors may result in stressful conditions for the growth of leafy green vegetables such as spinach. These environmental changes may influence the ability of human pathogens in association with other bacterial community members to persist in the phyllosphere [9], and the overall changes temperature, soil water content, UV radiation, and relative humidity can elevate the problems of salinity in soil [10].
Impacts of salinity on plants and microbes are generally through selective ion toxicity [11,12]. Changes in soil salinity may be one of the parameters impacting the use of recycled wastewater for agricultural irrigation [13]. In the next few years, the semiarid region of southwestern USA would likely experience increases in temperature, lesser precipitation, and more severe droughts [14,15]. Our understanding of leaf microbiota as providers of specific functions, for example, pathogen exclusion [16] and nitrogen fixation [17], relies on continued efforts to catalog the bacterial communities on plant foliage. Data analysis involving changes in environmental parameters influencing microbiome community changes often use the 16S rRNA gene sequencing methods while functional information on these changes requires analysis based on shot gun sequencing [18]. However, 16S rRNA data may be used to predict the functional attributes of bacterial assemblages based on new bioinformatic tools such as PICRUSt [19], Tax4Fun [20], and Piphillin [21]. The main objectives of this study were to determine bacterial composition at different seasons, and to examine functional capabilities of bacterial communities on the phyllosphere induced by seasonal changes and salinity.

Experimental treatments
These studies were conducted outside in closed recirculating tanks [22] filled with a mixture of loamy sand and peat moss (here fore designated as sand tanks). The system consists of 24 experimental plant growth units at the U. S. Salinity Laboratory in Riverside, CA. The experiments were conducted to determine the interaction of salinity and drought treatments on phyllosphere community composition. Experiments 1 started on 7 December 2012 and end in March 2013, and experiment 2 started on 14 March, 2013 and ended in May 2013 in large sand tanks as previously described [22]. Seeds were seeded and seedlings thinned to 25 plants per row in sand culture tanks as previously discussed [23] (S1 Fig). Experiment 3 was conducted in smaller sand tanks as previously described [24] and started in late April 2013 and ended in June. All tanks were irrigated with modified half-strength Hoagland's nutrient solution combined with various salinity levels. Starting nutrient solution was made up in Riverside tap water to flush the system. The three experiments utilized modified half Hoagland's solution with (in mM): 2.5 Ca (NO 3 4,7,9,12,15 dS m −1 by adding CaCl 2 , MgCl 2 , NaCl, Na 2 SO 4 to the base tap water-nutrient solution, and the salt concentrations used were based on EXTRACT Chem model [25] to predict the ion composition needed to achieve the target EC values. The control solution was without the added salinity nutrient and was maintained at an EC of 0.85 dS m −1 . Salinity treatments started after the first pair of true leaves was fully expanded on all the plants as previously reported [23] as well as tank irrigation. The three experiments were randomized design with three replications [22,23] including control with an EC at 0.85 dS m −1 and two different water types dominated by sulfate and chloride ion during the first experiment, while the second and third experiments only chloride ion water type. In these subsequent experiments we used only chloride water type since the first experiment did not show any statistical differences in spinach yields at any salinity levels between sulfate and chloride water types [22]. The concentrations of Na, K, Mg, Ca, and total-S [22] on the phyllosphere samples were as previously determined [22]. The average temperatures (˚C) and reference evapotranspiration (ET 0 ) that occurred during the experiment was acquired from the California Irrigation Management Systems (CIMIS) weather station no. 44 at University of California Riverside (S2 Fig), California [26]. About five leaf samples from one plant in the middle of the plot were cut above the soil surface with sterile blade and placed in stomacher bags and weighed. A total of five plants were sampled per plot, and the DNA from these plants were pooled into one tube per replicate for sequencing. A total of three replicates were used for each treatment. Total bacterial community DNA was recovered from the plant material by homogenization with 100 ml of PBS for 2 min at 260 rpm in a Seward Stomacher 400 Circulator (Seward Ltd., London, UK). Concentrated samples were used for isolation of the genomic DNA used for MiSeq sequencing.

DNA extraction and V4 16S Illumina MiSeq sequencing
Phyllosphere DNA samples were extracted with the Power soil DNA Kit (MoBio Laboratories, Solana Beach, CA) from triplicate plots and stored at -20˚C after further cleanup steps with DNA Clean and Concentrator (Zymo Research Corp-Irvine CA). DNA was quantified using a Nanodrop ND-2000 C spectrophotometer (Nanodrop Technologies, Wilmington DE) and Qubit HS kit (Fisher Scientific) and run on 1.0% agarose gel. DNA samples were sent to Second Genome (San Bruno, CA, USA) for V4 16S Illumina MiSeq sequencing. The V4 region of the 16S rRNA genes was amplified be Second Genome using fusion primers 515F (5 0 -GTG CCAGCMGCCGCGGTAA-3 0 ) and 806R (5 0 -GGACTACVSGGGTATCTAAT-3 0 ) [27], for sequencing using the MiSeq (Illumina, San Diego, CA, USA) instrument. Due to the presence of plant chloroplast DNA in surficial vegetable microbial metagenomes that may be confused with bacteria DNA, plant-derived chloroplast sequences were downloaded from Greengenes (May 2013 release) and removed before sequence analysis, and the sequence were rarefied to a standard number based on the lowest sequenced size [28] (Guron et al., 2019). Sequence data are available on this link to the general public: https://ars.usda.gov/ARSUSERFILES/ 20361500/PUBLIC%20DATA/DNA-SEQUENCE-OTU.DOCX. database [32]. Sequenced paired-end reads were processed using USEARCH [33], and clustered at 97% similarity by UPARSE (de novo OTU clustering). Clustering, chimera filtering, quality filtering was done as previously discussed [23]. Alpha-diversity was calculated to estimate sample richness and Shannon diversity and Beta-diversity metrics was calculated for the inter-comparison in a pair-wise fashion to determine dissimilarity score and store it in a distance dissimilarity matrix. Abundance-weighted sample pair-wise differences were calculated using the Bray-Curtis dissimilarity. All analyses were generated using Second Genome R package (vegan: R package version 2.2-1) as previously reported [23]. Univariate differential abundance of OTUs was tested using negative binomial noise model for the over dispersion and Poisson process intrinsic to this data, as implemented in the DESeq2 package [34], and as described for microbiome applications [35].

Inference of metagenomes
Piphillin was used to leverage the most up-to-date genome database for metagenome prediction from 16S rRNA sequence data [21], using the Kyoto Encyclopedia of Genes and Genomes (KEGG) (KEGG 70.1 database). A genome was inferred for each 16S rRNA OTU based on the sequence identity between an OTU's representative sequence and the nearest neighbor 16S rRNA sequence from the genome databases restricted to a minimum identity of 97%. OTU abundance was normalized by 16S rRNA copy numbers and then multiplied by the gene contents of each inferred genome to predict each sample's metagenome.

KEGG Pathways and genes significance testing
To identify differentially abundant pathways and genes, a Wilcoxon Rank Sum test was employed. Where samples could be paired across categories, a paired Wilcoxon Signed Rank test was employed. P-values were adjusted by Benjamini-Hochberg procedure to control for false discovery rates from multiple testing.

Result and discussion
Alpha diversity index (Shannon index of diversity) was significantly higher in June (Table 1). Our study agrees with a study of leaf age and seasonality of the tree species Quercus ilex in the Mediterranean forest which revealed an increase in the richness of epiphytic bacteria with increasing time of colonization in summer [36]. Other studies have suggested that the bacterial community composition in leafy greens change with season and not solely with leaf maturation [9]. Our study is one of the few studies that show seasonal effects of phyllosphere bacterial community changes on leafy greens and confirmed the findings of Williams et al. [9], where, for the first time, a seasonal effect was reported. The study by Williams et al. [9] was conducted in the Salinas Valley, CA, USA with Romaine lettuce and two irrigation regimes, whereas this study was conducted in southern California with drip irrigation and three seasons. Most of the leafy green productions in Salinas Valley are grown from June to September, whereas leafy green production in southern California, USA is generally grown during the winter months

Comparison of leaf microbiomes in response to salinity
Salinity did not significantly influence bacterial beta-diversity on the phyllosphere using weighted abundances (Fig 1).
Also, salinity did not affect community composition (PERMANOVA P = 0.501; Fig 1) based on the dimensional reduction of the Bray-Curtis distance between microbiome samples, using the PCoA ordination method. Samples were separated by primary axis (Axis 1) and secondary axis (Axis 2) by collection months with March and May samples separating from June samples. March and May samples had more heterogeneity than June samples, as June samples clustered closer to each other. Overall, 67.2% of the sample variation was explained by the major axes. At the family-level abundances Halomonadaceae (P = 0.0099) and Weeksellaceae (P = 0.0253) were significantly affected by salinity (Fig 2). There were no significant differences from other families based on Spearman's test (e.g. Enterobacteriaceae: P = 0.671; Pseudomonadaceae: P = 0.272; Sphingomonadaceae: P = 0.754; Hyphomicrobiaceae: P = 0.624; Geodermatophilaceae: P = 0.737).
Detail analysis was conducted to determine the effect of salinity on bacteria to the species level based on bacterial OTUs. This analysis identified 156 OTUs that were significantly (P � 0.05%) impacted by salinity (S1 Table). About 91 OTUs had significant decrease in their F statistics with their log-2fold-change also indicating negative statistics between the first experiment that ended in March and the last experiment that ended in June while 65 OTUs were significantly enriched with positive statistics. Some of these OTUs were unidentified or unclassified, with only thirty-three identified to the genus level and eight to the species level that were significantly enriched (P � 0.05%) due to the impacted of salinity during the study period (Table 3). This data provided additional insight into data presented in Fig 2 with Halomonadaceae and Weeksellaceae as two of the families that were significantly enriched by salinity. It also provided additional insight into the data that other bacterial species were also significantly affected by salinity during the study. For instance, specific species that were significantly enriched by salinity included Sphingobacterium multivorum (0.57-fold, P<0.00025) Acinetobacter Iwoffii (0.59-fold, P < 0.0004), Methylobacterium adhaesivum (-0.36-fold,

Comparison of leaf microbiomes in response to season
Phyllosphere microbiomes were significantly (P = 0.003) affected by season. Air temperatures in experiment 1, 2 and 3 were 15.9˚C, 17.9˚C, 20.2˚C, respectively [22,37]. Leaf surface temperature moderately affect phyllosphere microbiomes (P = 0.057), as well as daily mean evapotranspiration (ET0) significantly (P < 0.002) affected leaf surface microbiota during the three growing seasons. Daily mean evapotranspiration (ET0) values of each experiment and growing season as number of days was recently reported in a related study [22]. Seasonal water demand during the three experiments vary; ranging in very low in experiment 1, averaging 1.96 mm during the first ten days. During experiments 2 and 3 the minimum ET0 values were 4.43 and

PLOS ONE
6.36, respectively. It should be noted that the growth period became shorter in the successive seasonal experiments as temperature increased, thus cumulative potential evapotranspiration (PET) decreased as previously reported [22,37] from 276.3, 240.6 and 202.1 mm with increasing temperature for experiment 1, 2 and 3, respectively. High salinity caused yield loss and decreased all gas exchange and vegetative parameters as previously reported [22], and the high salinity may also select for bacteria population with high tolerant to salinity. The effects of temperature based on season on the most abundant families were analyzed and grouped according to their relative abundances (Fig 3). Wilcoxon signed rank sums and mean relative abundances for family-level taxa between temperature groups (months) are provided in Table 4. The influence of temperature on relative abundance at the family level was not significantly different for all the eight major families between the collection months of March and May based on Wilcoxon signed rank sums and mean relative abundances for family-level taxa (Table 4). No significant effects were also observed between March and June for the eight major taxa. However, seven of the top eight families tested were significantly different between May and June ( Table 4).
The effect of season which was associated with increases in temperate showed significant differences (P � 0.05%) in the number of OTUs impacted by changes in season due to

PLOS ONE
temperature changes (S2 Table). Very high numbers of OTUs (709) were significantly affected by season in comparison to only 156 that were impacted by salinity. However, only 255 OTUs were impacted negatively by season whereas 454 were significantly enriched. At the genus level, 101 bacteria were identified including 15 identified to the species level suggesting greater effect of season on phyllosphere microbiome than salinity (Table 5). This was also reflected by the high number of log2-folg change as reported in Table 5 in comparison to the smaller changes in salinity. As seen in Table 4 the significant effect of season was observed in six families between the month of May and June harvesting, but none during the first two planting. Our data suggest that more individual bacterial species increased in density than decreased due to season than salinity during the study. In this study salinity significantly affected the family Halomonadaceae (Gammaproteobacteria) (P < 0.009) and Weeksellaceae (P < 0.025) in the phyllosphere samples throughout the three experiments (Fig 3). Halomonadaceae family such as Halomonas elongata, is a wellknown example of a bacterium that can adapt to life over the whole salt concentration range    (halophilic) from near fresh water to halite saturation [38]. The unclassified Halomonas, Sphingobacterium multivorum and Acinetobacter Iwoffii identified in this study had strong positive log2-change, suggesting that these species were not negatively impacted by salinity. In this study, the effects of soil salinity, water content, and how changes in seasons due to warming trends impact bacterial communities on leaf surfaces may provide some guidance into future research in management of saline soils [22,37,39] and changes in leafy green surface microbiome. Therefore, understanding the effect of changes in salinity and water content on soil microorganisms is important for crop production, sustainable land use, and rehabilitation of saline soils [40]. The different ways in which various climate drivers such as temperature, precipitation, ET0 and their interactions with seasons, irrigation methods, and leaf age might affect phyllosphere microorganisms and their activities [9], presents a serious challenge in understanding the impact of drought and salinity on leaf surface bacterial communities. Our study agrees with the work of William et al. [9] emphasizing the difficulties in predicting with full certainty the effects of soil salinity and seasons on leafy green bacterial communities until we collect more long-term data from both control and field studies. As seen from our previous work, the impact of drought/temperature and salinity on bacterial communities or changes in climate drivers may be mostly associated with the rhizosphere [37]. Furthermore, it is expected that the southwestern USA will experience increases in temperature, receive less springtime precipitation, and have more frequent and severe droughts [15]. Overall, climate-induced trends will impact the persistence and dispersal of foodborne pathogens in myriad ways, especially for environmentally ubiquitous microorganisms [14]. Drought can result in decreased bacterial diversity [41], proliferation of pathogenic species [42], and increased survival of an introduced E. coli O157:H7 derivative [43,44]. Along these lines, drought may inhibit some populations of native microflora while allowing for some of the more resistant groups of pathogenic bacteria, such as Staphylococcus, Clostridium, and Bacillus to survive [14]. These authors noted that most foodborne pathogens are likely to be negatively affected by drought due to their dependence on water activity for survival and growth. Climatic factors such as temperature and season [45][46][47] may play an important role in shaping community structure.

Population of potential pathogenic bacterial sequences on leaf surfaces
Using Greengenes database (http://greengenes.lbl.gov) to identify potential pathogenic bacteria sequences at the genus level (Table 6), we identify sequences belonging to 18 potential human pathogens and two plant pathogens (Erwinia and Rastonia) among others. These two plants pathogens were included because they have been shown to degrade leaf materials and release carbon sources for pathogens such E. coli O157 to use and proliferate on leaf surfaces [48]. The relative abundance of the 20-genus identified during the three seasons showed, Bacillus, Pseudomonas, Flavobacterium, Halomonas, and Erwinia as indicated in Table 6, at the highest abundance levels. Both Pseudomonas and Erwinia occurred at the highest concentrations in May while Bacillus in June. However, most of the sequences associated with potential human pathogens such as Clostridium, Mycobacteria, Legionella, Nocardia, Burkholderia, Corynebacterium, and Enterococcus were least abundant in June than in March and May.
Halomonas sequences were detected in higher numbers following a warming trend from March to June (Table 6). Halomonas is a newly recognized human pathogen causing infections and contamination in a dialysis center [49] and can grow in a wide range of salt concentrations [38]. The detection of sequences for Arcobacter and Acinetobacter were also of increasing magnitude from March to June (Table 6). Furthermore, A. butzleri displays microbiological and clinical features like those of Campylobacter jejuni; however, A. butzleri is more frequently associated with watery diarrhea [50]. Acinetobacter may also be ubiquitous in the environment and this pathogen is a serious infectious agent of importance to hospitals worldwide [51]. It also can accumulate diverse mechanisms of resistance, leading to the emergence of strains that are resistant to many antibiotics [52]. It should be noted that some of these are human or plant pathogens but many are rather beneficial commensals. and Sphingomonas as dominant genera [47,53]. It should be noted that bacteria establishing colonies at the phyllosphere are limited by various factors including both biotic and abiotic factors. Abiotic factors such as the available nutrient [1], seasonal variation, rainfall, temperature, plant immunity, and competitor microbes [54] may influence the survival of microbes in the phyllosphere. As previously reported, environmental factors can drastically influence the microbiome changes on phyllosphere. This is common to epiphytic microorganisms exposed to heavy stress during the seasonal cycle, the day/night cycle, and the growth, age, and anatomical dynamics of the plant. As previously reported drought condition can influence the epiphytic microbial community on Holm oak [55]. At the same time, during high temperature, bacterial endophytic communities are altered in lower leaves of paddy, but not in the epiphytes [56], which was in contrast to epiphytic fungal community that responded well during warmer seasons [44,57,58]. Pathogens survival in host tissues can use hemibiotrophs and necrotrophs mode of life. Some chemicals of plant tissues inhibit the microbial association either through salicylic acid or jasmonic acid type and the reactive O 2 species may have an inhibitory effect on the pathogens [59]. Plants use jasmonic acid, methyl jasmonate, ethylene, flavonoid, 12-oxo-phytodienoic acid, and salicylic acid-mediated signals for quenching pathogens on its surface. Bacteria need carbon, nitrogen, inorganic, and organic energy sources. However, the phyllosphere is a very harsh environment, and in the absence of some of these nutrients, the phyllosphere is still colonized by a large number of bacteria (10 5 -10 7 CFU/g of the leaf) in the presence of high relative humidity and free water at suitable environmental conditions [60,61]. This is due to the release of nutrients or leaf exudates which adequately can support microbial growth. There are varieties of molecules leached from the plant leaves such as sugar, amino acids, organic acids, minerals, etc. that will also support bacterial growth [62][63][64]. These leaching materials may differ with plant species and the environmental condition [62,63,65]. Therefore, it is not surprising that in a study like ours, we were still able to identify more than 5000 OTUs using next generation sequencing, and this is in agreement with reports that on the average the phyllosphere may contain large amount of bacteria (10 5 -10 7 CFU/g of the leaf) [61,62].

Pathways and genes of phyllosphere community
Shotgun sequencing is the direct method for the analysis of functional capabilities of environmental microbiome. However, changes in metabolic pathways and functional genes based on Piphillin algorithm predictions [21] to generate inferred metagenomes of sample type communities based on a 97% identity cutoff has recently been introduced [21]. Approximately 25% of the sample type microbial communities were assigned to genomes by the Piphillin algorithm. This is very low, however other studies have captured the same low correlation between the distance metrics based on phylogenetic/taxonomic and predicted the metagenomic composition of their samples [66,67].
The top pathways on leaf samples during the three growing seasons were located to biosynthesis of amino acids, ABC transporters, ribosome, aminoacyl-tRNA biosynthesis, two-component system, carbon metabolism, purine metabolism, and pyrimidine metabolism (Fig 4A  and 4B). There were 161 significantly different functional pathway features detected out of 308 tested. Only the 37 functional pathways with an absolute log-2-fold change greater than 3 and not equal to Inf or -Inf are shown. The photosynthesis antenna proteins pathway was significantly enriched in June leaf samples when compared to March (Fig 4C) based on KEGG pathway functional prediction feature selection. Features selections were considered significant if their FDR-corrected P-value was less than or equal to 0.05, and the absolute value of the log-2-fold change was greater than or equal to 1. A total of 157 pathways were enriched in March leaf samples and 4 pathways were enriched in June leaf samples, while 121 pathways were significantly different between the May and June leaf samples. Two pathways were enriched in June samples and 119 pathways were enriched in May samples (Fig 4D). Only the 29 pathways with an absolute log-2-fold change greater than 4 are shown.
At the gene level, the most abundant genes for March-Junes leaf samples were mcp (methyl-accepting chemotaxis protein), tRNAs genes which included tRNA-alanine, tRNAarginine, tRNA-glycine, tRNA-leucine, tRNA-methionine, tRNA-serine, and tRNA-valine (Fig 5A), while methyl-accepting chemotaxis protein (mcp), RNA polymerase sigma-70 factor (ECF subfamily; SIG3.2, Sig3.2, rpoE), iron complex outer membrane receptor protein (TC. FEV.OM), polar amino acid transport system permease protein (ABC.PA.P), multiple subunits of the peptide/nickel transport system permease protein (ABC.PE.P1;), and tRNA Leucine were associated with the most abundant genes in May versus June (Fig 5B). There were 3,497 KEGG ortholog significantly different functional pathway prediction features for March versus June that were detected out of 7,127 tested ( Fig 5C). Overall, 2,697 genes were enriched in March leaf samples and 800 were enriched in June leaf samples. Only genes with an absolute log-2-fold change greater than 6 and not equal to Inf or -Inf are shown. For May versus June leaf samples, there were 3,033 significantly different features detected out of 7,132 tested. A  Fig 4A); May vs. June (Fig 4B). KEGG pathway prediction feature selection for March vs. June leaf (Fig 4C). There were 161 significantly different features detected out of 308 tested. Only the 37 pathways with an absolute log-2-fold change greater than 3 and not equal to Inf or -Inf are shown. KEGG pathway prediction feature selection for May vs. June leaf (Fig 4D). Features were considered significant if their FDR-corrected P-value was less than or equal to 0.05, and the absolute value of the log-2-fold change was greater than or equal to 1. There were 121 significantly different features detected out of 308 tested. Only the 29 pathways with an absolute log-2-fold change greater than 4 are shown.
https://doi.org/10.1371/journal.pone.0252242.g004 total of 291 genes were enriched in June samples and 2,742 genes were enriched in May samples. Only 29 genes with an absolute log-2-fold change greater than 6.5 and not equal to Inf or -Inf are shown (Fig 5D). Several genes related to toxin co-regulated pilus biosynthesis proteins were also significantly enriched in June leaf samples, when compared to March and May leaf samples. Overall, several pathways and genes were enriched by different bacterial OTUs with increasing salinity. This may be a function of microbial consortia as individual isolate may not perform all the functions from breaking down the main components and their metabolites in the saline environment [68]. In this study, Halomonas was significantly detected in higher numbers following warning tread from March to June in the high saline environment. This suggests that the genus was able to degrade nutrients in the environment for growth. According to Wang et al. [69], the genus Halomonas was among a group of bacteria that contributed to phenanthrene degradation intermediates through catechol 1, 2-dioxygenase (C120).
There are a lot of flexible metabolic adaptations for microbes in the phyllosphere, and this is why they survive in the harsh microenvironment. Plant releases a lot of compounds for  (Fig 5A). Plot shows the most abundant functional genes. Proportional abundance of the top inferred genes for May vs. June leaf (Fig 5B). Plot shows the most abundant functional genes. KEGG ortholog prediction feature selection for March vs. June leaf (Fig 5C). There were 3,497 significantly different features detected out of 7,127 tested. Only genes with an absolute log-2-fold change greater than 6 and not equal to Inf or -Inf are shown. KEGG ortholog prediction feature selection for May vs. June leaf (Fig 5D). There were 3,033 significantly different features detected out of 7,132 tested. Only 29 genes with an absolute log-2 fold change greater than 6.5 and not equal to Inf or -Inf are shown. Features were considered significant if their FDR-corrected p-value was less than or equal to 0.05, and the absolute value of the log-2 fold change was greater than or equal to 1.
https://doi.org/10.1371/journal.pone.0252242.g005 metabolic functions such as carbohydrates, polyols, amino acids, amines, isoprenoids, halogenated compounds, alcohols, water and salts, and these are the main sources of nutrients for phyllosphere microorganisms [70] as well as saline or alkaline pH which generates stress in phyllosphere microbes. Also, phyllosphere microbes can develop multiple mode of adaptation to survive in phyllosphere such as tolerance, antimicrobial, extracellular polysaccharides, phytohormonal, and immunity compounds against a microbial competitor [70] (Trouvelot et al. 2014). The phyllosphere microbiome acts as a vital role for leaf surface environment and their surrounding ecosystem functions [71]. It has been shown that plants release a variety of volatile organic compounds (VOCs) and its precursors on the surface of leaves such as terpenes, monoterpenes, flavones, methanol, methane, and halogenated methane [72], and these could regulate the microorganisms in response to changes in the environment. Many of the phyllosphere microbial communities in this study share the common metabolic properties of the soil and rhizosphere microbes from our previous study [73]. Phyllosphere bacterial communities such as Halomonas, Pleomorphomonas, Prosthecobacter, Geodermatophilus, Pirellula, Blastococcus, Opitutus, Allochromatium, Desulfosporosinus, Burkholderiaas, etc., were also identified in soil and rhizosphere. These bacteria also have carbohydrate metabolizing genes involved in utilization of starch, hemicellulose, pectin, and cellulose, rich in humus materials [74][75][76]. Therefore, there seems to be some form of synergistic activities involving the above and below ground plant microbiomes.

Conclusion
In this study, seasonal variations significantly affected bacterial composition on leaf surfaces. Seasonal differences were major factor driving most of the sequences associated with potential human pathogens such as Clostridium, Mycobacteria, Legionella, Nocardia, Burkholderia, Corynebacterium, and Enterococcus which were less abundant in June than in March and May. Also, photosynthesis antenna proteins pathway was significantly enriched in June leaf samples in comparison to March, and the high enrichment may be associated with stress conditions such as high light, high salinity, elevated temperature, and nutrient limitation. Therefore, the application of functional pathway analysis has provided another evidence of spinach response to different stress conditions examined in this study during the different seasons and salinity.