Introduction

F1000Research

2046-1402

F1000 Research Limited

London, UK

10.12688/f1000research.12490.1

Research Article

Articles

Agriculture & Biotechnology

Community Ecology & Biodiversity

Plant-Biotic Interactions

Plant-Environment Interactions

Maize phyllosphere microbial community niche development across stages of host leaf growth

[version 1; peer review: 1 approved with reservations]

Manching

Heather C.

Formal Analysis Investigation Methodology Supervision Writing – Review & Editing 1 2 Carlson

Kara

Formal Analysis Investigation Methodology Writing – Review & Editing 2 3 Kosowsky

Sean

Formal Analysis Software Writing – Review & Editing 2 Smitherman

C. Tyler

Data Curation Formal Analysis Investigation Methodology Resources Writing – Review & Editing 2 4 Stapleton

Ann E.

Conceptualization Data Curation Methodology Project Administration Supervision Writing – Original Draft Preparation Writing – Review & Editing https://orcid.org/0000-0001-7242-2673 a 2 1Department of Plant and Soil Sciences, University of Delaware, Newark, DE, USA 2Department of Biology and Marine Biology, University of North Carolina Wilmington, Wilmington, NC, USA 3North Carolina Museum of Science, Raleigh, NC, USA 4Nestlé, Inc., Danville, NC, USA

a stapletona@uncw.edu

Competing interests: No competing interests were disclosed.

18 9 2017

2017

1698

19 6 2026

2017

This is an open access article distributed under the terms of the Creative Commons Attribution Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Background The phyllosphere hosts a variety of microorganisms, including bacteria, which can play a positive role in the success of the host plant. Bacterial communities in the phylloplane are influenced by both biotic and abiotic factors, including host plant surface topography and chemistry, which change in concert with microbial communities as the plant leaves develop and age. Methods We examined how the Zea mays L. leaf microbial community structure changed with plant age. Ribosomal spacer length and scanning electron microscopic imaging strategies were used to assess microbial community composition across maize plant ages, using a novel staggered experimental design. Results Significant changes in community composition were observed for both molecular and imaging analyses, and the two analysis methods provided complementary information about bacterial community structure within each leaf developmental stage. Conclusions Both taxonomic and cell-size trait patterns provided evidence for niche-based contributions to microbial community development on leaves.

epiphytes plant host community assembly Zea mays L. ecological niche scanning electron microscopy

The author(s) declared that no grants were involved in supporting this work.

Introduction

The phyllosphere, the microbial community inhabiting the surface of plant leaves ¹, can play important positive roles for the plant partner. For example, disease resistance is increased when supportive bacteria are present ^{2,
3} and increased host plant growth in the presence of specific community members has been documented ^{4–
6}, along with a positive effect of overall increased microbial diversity ⁷. Leaf habitats are also extremely important for the microbial community, as this habitat is quite large in surface area globally. In temperate regions seasonal growth and annual agricultural plantings lead to repeated opportunities for colonization by adapted microbes ^{8,
9}.

Determinants of community assembly/succession in microbial systems include biotic and abiotic factors, often classified into different community assembly rules ¹⁰. Typically random or neutral stochastic models are tested against abiotic effect models, which are called environmental, or, more recently, habitat filtering rules. Environmental filtering can produce different patterns depending on scale ¹¹. As recommended by Kraft et al. (2015) we will refer to environmental filtering as habitat filtering when there is no separate test for abiotic and biotic effects. A second class of assembly rules addresses the interactions of organisms, by either facilitation or competitive exclusion. These niche interactions are often modeled separately from abiotic effects, or abiotic effects are assumed to happen first.

Community assembly and succession processes are typically inferred from diversity patterns in community samples ¹². The number of different taxa (alpha-diversity), the difference in the composition of taxa between replicates (beta-diversity) and the proportion of co-occurring taxa all provide evidence for these different community-generating mechanisms. Niche construction is inferred from a pattern of decreasing beta-diversity, increasing co-occurrence, and increasing richness, while habitat filtering is characterized by decreasing beta-diversity with no trends in co-occurrence or richness ^{13,
14}. In contrast, neutral community assembly would exhibit increased richness (provided immigration rates are maintained), with unchanged beta-diversity and co-occurrence levels. Once we have measurements of the various components of diversity, we can compare observed trends to these expectations. The combination of increasing/decreasing trends would then support or fail to support a particular mechanism of community succession. Both trait-based (morphological, in our study) and taxonomic measures of diversity can be compared to the expected trends ^{15,
16}. For microbial systems, cell size is one of the few traits that is measurable at the microscopic scale; other functional traits, such as metabolic activity, are more often examined at the molecular scale, by whole-community DNA sequencing, proteome analysis, or metabolite profiling.

Previous analyses of controlled inoculation phyllosphere community development indicate that habitat filtering is important, with niche effects less often detected across time, though all previous studies have used seasonal or through-time designs that confound exposure time and plant development. Pre-colonization reduces immigrant success ^{8,
17}, suggesting biotic competition is key, though high beta-diversity in replicate samples from both pre-colonized and newly inoculated communities suggested that there are many options for species that arrive first to succeed. Several crop species have been tested for seasonal patterns of phyllosphere community development. Microbial populations increased more rapidly on younger leaves of lettuce ¹⁸, and alpha-diversity increased over the first few weeks then plateaued ¹⁹, which provides evidence that habitat filtering was a contributor to community succession later in the season. In apple flowers, species richness increased early in the season and beta-diversity dropped, then rose ²⁰. About half the apple flower microbial species had significant co-occurrence relationships and the numbers of co-occurrences increased and decreased over time. This pattern supports niche construction as a key mechanism for apple flower colonization. As in apple flowers, across-season dicot crop phyllosphere samples became less diverse (alpha-diversity went up then declined). In the crop samples, more shared and plant-specific taxa were found later in the season ²¹, again supporting development of specific community structure. In contrast, Arabidopsis epiphyte beta-diversity increased with time ²², providing stronger support for stochastic effects in this system. The history of annual opportunities for new colonization is reflected in the processes for dispersal and succession on plant leaves. Colonization of new growth in the presence of older leaves may proceed differently than colonization from distant source materials. The contributions of abiotic changes and host leaf age are confounded when sampling occurs through time, as in seasonal sampling, thus necessitating alternative experimental designs and/or additional phylogenetic or functional trait information to examine the roles of these assembly processes. Sampling serially across time, as all prior studies have done, could confound short and long-term trends and does not enable differences in plant host to be distinguished from inoculation processes ^{19,
20}. Therefore, for this study host plant seeds were planted at different times and sampled once.

How can leaf microbial communities be measured? Molecular methods for detection of species using ribosomal fragments or sequences have become increasingly common ^{20,
21,
23}. Microscopy can be used to detect bacteria (and fungi) on corn leaf surfaces ²⁴, and surface bacteria can be followed with reporter genes ²⁵ as well as by direct enumeration. Microscopic or particle-based enumeration methods are especially useful for microbial communities where we have limited information on traits and horizontal gene transfer is likely ²⁶. For this study we used both a molecular approach to measure relatedness-based (taxonomic) microbial community composition, and direct enumeration of scanning electron micrograph (SEM) images to measure cell-size trait changes.

Our specific questions concerned the pattern of community development – how community succession varies with plant age, and whether the microbial community was primarily shaped by microbe-host interactions (environmental or more generally, habitat filtering) or by microbemicrobe interactions (niche construction). We found complex patterns in the maize leaf microbial community structure even across the short series we sampled, with niche construction and stage-specific increasing-decreasing diversity trends observed.

Methods Experimental design and sample collection

Zea maysL. B73 inbred seeds obtained from the Maize Genetics Cooperation Stock Center ( http://maizecoop. cropsci.uiuc.edu/) were planted using a hand held planter between April 26th and June 15th 2009 at the Central Crops Research Station in Clayton, NC ( http://maizecoop.cropsci.uiuc.edu/). All necessary permits were obtained for the described study, which complied with all relevant regulations; permits were managed by the research station staff. Approximately twenty seeds were planted each date, with plantings between four and eleven days apart in a single field. All plants within a developmental stage were planted in the same row, bordered on either side by older and younger stages, then by additional crops, grass borders, and tree-shrub hedgerows.

All ten replicate leaf samples were collected at once on July 15, 2009. Both scanning electron microscope (SEM) and ribotype (ARISA) samples were obtained from tips of separate central leaves from ten plants at each developmental stage. SEM samples were approximately five centimeters square and were cut along the leaf mid-rib near the tip ( Figure 1). Sections were immediately fixed, by immersion in a solution of 2.5 percent glutaraldehyde in 0.1 M Sorenson’s phosphate buffer at a pH of 7.1.

Figure 1. Conceptual diagram of the leaf sampling procedure.

Each plant in the drawing represents 10 plants from one of seven developmental stages; the experimental design had n=10. Small squares represent SEM sampling areas and large squares represent ARISA sampling areas. All samples were taken from center leaves within each developmental group, and both samples per plant were taken from leaves of comparable position based on height.

Leaf wash samples for ribosomal DNA analysis were then collected from the same ten plants, using adjacent leaf numbers and a similar portion of the leaf tip ( Figure 1). Approximately 10 cm of leaf was cut, submerged in buffer (0.2mM Tris pH 7.5, 0.02 mM EDTA, 0.00012X Triton X100), agitated vigorously and removed using the method described in 27. Leaf washes were stored on dry ice while in the field then placed in a -70 degrees C freezer for long-term storage. One mock control (where no leaves were placed in the tubes) and one wash control (SEM sample from washed leaf to check for washing effectiveness) was taken for each stage.

SEM sample processing

Excess glutaraldehyde was removed from leaf SEM samples by rinsing for 30 minutes in 0.1 M Sorenson’s phosphate buffer followed by an additional 30 minutes in distilled water. Sections were then dehydrated in an acetone series (50, 70, 90, 100 percent), each lasting 30 minutes, followed by a 30-minute treatment in hexamethyldisilazane. Sections were then allowed to air dry before being attached by the lower leaf surface to stubs using double-sided sticky tape. One leaf section chosen arbitrarily from each plant was sputter coated with 10 nm of platinum/palladium and observed in the SEM at 5 kV with a spot size of 3. The SEM procedure was tested prior to sampling using young (approximately 25 cm tall), greenhouse grown maize plants and Pseudomonas aureginosa ²⁸.

SEM image acquisition

Using the SEM stage control, five evenly dispersed random images (quadrants) per plant were taken at 1000x and 2000x for 35 sections (35 plants) each to give two viewpoints by which to count bacteria. All images are available from the image repository BisQue via the iPlant/CyVerse Data Commons, at doi: 10.7946/P2WC77.

Ribosomal spacer amplification and separation (ARISA)

For the DNA analysis, leaf washes were syringe filtered using previously developed methods ²⁹ (videos of the method are available upon request) and microbial DNA was extracted using an UltraClean Microbial kit according to the manufacturer’s instructions (MoBio Laboratories, Carlsbad, CA). An equal volume of each DNA sample was amplified to ensure fair representation of the epiphytes collected.

Two separate, nested PCR reactions were performed on experimental samples as described by Cardinale et al. ³⁰, along with four positive controls from a mixture of DNAs from cultured bacteria (Carlson, 2015) and four negative (nuclease-free water) controls. For the initial PCR step, a master mix was made to include: Q5 master mix (New England Biolabs, Ipswich, MA), ITSF primer (5’-GTCGTAACAAGGTAGCCGTA-3’), LD primer (5’-CCGGGTTTCCCCATTCGG-3’), and nuclease-free ultrapure water. Ten microliters of master mix was added to sterile PCR tubes followed by experimental and sample controls. The samples and controls were then run on an iCycler (BioRad, Hercules, CA) PCR machine with the following settings: 1X 95C for 30 seconds, 35X (95C for 30 seconds; 61.5C for 30 seconds; 72C for 1 minute, 30 seconds), 1X 72C for 10 minutes. For the second PCR step, a second master mix was made with the following: Q5 master mix, ITSreub primer (5’-GCCAAGGCATCCACC-3’), which contained a fluorescent HEX marker for fragment analysis in a Genetic Analyzer, SD primer (5’-TGCGGCTGGATCCCCTCCTT-3’), and nuclease-free water.

Samples were run on the Applied Biosystems Genetic Analyzer (ABI, Inc, Carlsbad, CA) and analysis of PCR product length was carried out using the instrument software. For each sample, two microliters of PCR product and 13 microliters of GS500-Rox plus HiDi (ABI, Inc, Carlsbad, CA) were combined. The samples were denatured at 94C for ten minutes, placed on ice for two minutes, then placed into the analyzer and the program started.

SEM image analysis

For each SEM quadrant image, Adobe Photoshop was used to manually outline each bacterial organism present (examples are shown in Figure 2). All single bacteria present in the image were outlined; bacteria were identified by size and smoothness of circumference, as originally described by Davis ²⁴. Bacterial areas were measured using Image Pro Plus (Media Cybernetics, Inc., Rockville, MD, USA) and recorded for each quadrant picture that was taken. Image Pro Plus sized bacteria that were touching other bacteria as one large organism, therefore as many bacteria as possible were outlined in the clusters without outlining touching bacteria. Each image was individually measured using calibration settings based on the magnification used to take the image. Each image analyst completed analysis of a training image set and the coefficient of variation of their analyses was checked to ensure that it was less than 5 percent.

Bacterial sizes acquired from all images were binned using methods described by Shimazaki et al. ³¹. Bin sizes were assigned to the microbial samples as size-class units (SCU) and a microbial community profile was created for the samples using the RiboSort package in R ³².

Figure 2. Example SEM (left panel) and bacterial overlay images (right panel).

A) Typical image after 30 days of growth in the field. B) 72-day image, with leaf hair and stomata clearly visible. C) Typical image after 80 days of growth in the field.

ARISA data processing

Files from the Genetic Analyzer were opened in Peak Scanner 2 software ( http://www.appliedbiosystems.com) using Analysis Method=Sizing Standard-PP, and Size standard=GS500. Only the green (hex-primer-labeled) peaks were used to determine level of bacterial abundance. The size and heights of the green peaks were calculated using the ROX GS500 standard. The individual output files were assembled into a single microbial profile file using the RiboSort package ³² using the script, RiboSort(data,dataformat="standard",dye="G", output="proportions",zerorows=FALSE,repeats=1, mergerepeats="presentinall"). This program automatically bins the peaks to generate abundances; we used the default threshold of 50 fluorescence units as recommended by the program authors.

Statistical analysis of presence-absence matrices

The OTU (operational taxonomic unit) abundance matrices output from Ribosort were analyzed by Permdisp and Permanova ^{33–
35} using Primer-E software v6 ³⁶. Overall square root transformation and Bray Curtis similarity resemblance were applied to the initial abundance data within Primer-E before Permanova and Permdisp tests were run. The P-values from the program were Bonferroni-corrected. Primer-E SIMPER analysis was used to compute the number of OTU contributing to differences between stages. Primer-E functions were used to compute diversity measures and taxa numbers for each sample (Supplementary Dataset).

The evenness (Pielou’s J’) of OTU and SCU for each sample was computed using Primer-E v6 ³⁶. The values of J’ for OTU and SCU were compared across plant stages by rank-transformed analysis of means using the routines in JMP Pro v11 (SAS Inc, Cary, NC) with P-values adjusted for multiple comparisons. For the rank-transformed analysis of means the P-value threshold was set at the Bonferroni-corrected value of P less than 0.007.

This rank-transformed analysis of means is not sensitive to heteroschedasticity and is more powerful than other nonparametric tests when the data distribution is heavy-tailed ³⁷. Nonparametric pairwise comparisons of evenness between stages were done using the Steel-Dwass multiple-comparison method. The distributions of SCU areas in each stage of plant growth were also computed with JMP Pro v11.

Co-occurrence data analysis

Pairwise analysis was completed using the R package Co-Occur ³⁸; this package allows determination of the number of OTU/SCU that occur together in the same sample day at a probability above random, along with the number of taxa that are not found together (negatively correlated taxa) above the random-pair threshold. The analysis code, input and output files are available at https://github.com/skosowsky/Cooccur. We created network visualizations of the co-occurring OTU and SCU with Cytoscape (v. 3.4) and compared networks using the Cytoscape plug-in app DyNet (v1.0) ³⁹. Intensity of node connection change across the seven stages is reflected in the intensity of the color of the DyNet-analyzed nodes.

Statistical analysis of presence-only data

The number of OTU and SCU present within each sample was computed using Primer-E v6 (abbreviated as S, also known as Chao1). The number of OTU and SCU compared to the overall mean across plant stages was analyzed using rank-transformed analysis of means via routines in JMP Pro v11 (SAS Inc, Cary, NC), with P-values adjusted for multiple comparisons. For the ranktransformed analysis of means the P-value threshold was set at the Bonferroni-corrected value of P=0.007. This rank-transformed analysis of means is not sensitive to heteroschedasticity and is more powerful than other nonparametric tests when the data distribution is heavytailed. Nonparametric pairwise comparisons between stages were carried out with the Steel-Dwass multiple-comparison method within JMP.

Results

Our experimental design allowed for a single microbial sampling that was not confounded by differential exposure to recent weather events, with sampling for molecular and microscopic analyses in each plot ( Figure 1). Molecular analysis of community ribotype composition (ARISA) and image analysis of microbial community size cell assemblage size classes (by scanning electron microscopy, SEM) provided two complementary views of community assembly processes.

Presence-absence matrix comparisons

ARISA community composition and beta-diversity. Microbial community composition changed over time in bacterial ribosomal taxa (Permanova P-value of 0.026), as measured by differences in operational taxonomic units (OTU). The Permanova pairwise results for bacterial ribosomal OTU summarized in Table1 and with additional details in Supplementary Table 1 and Supplementary Table 2 indicate that, specifically, day 30 and day 80 were significantly different from each other in community composition. Thus we conclude there was a difference in either the richness or arrangement of the communities’ ribotypes in the earliest and latest stages of community development. To determine which aspect of community structure contributed to the difference detected by Permanova, we examined the dispersion as a measure of beta-diversity, i.e. the amount of range in the set of species present across replicates within each time. There was no significant change in the dispersion in the ribosomal OTU datatype across stages (Permdisp P=0.445, details in Supplementary Results). To further analyze the ARISA community structure, we compared the evenness of species distributions across stages.

Table 1. ARISA Permanova Community Structure Comparisons.

Significance groups with different letters are significantly different by pairwise test with a Bonferroni multiple-test-correction threshold of P=0.001315.

Stage (day after planting)	Significance group
30	A
41	A,B
48	A,B
58	A,B
62	A,B
72	A,B
80	B

Evenness is a measure of the range of OTU abundances in a sample unit; we used the Pielou’s J’ measure of sample evenness for comparison of stage evenness values. Evenness values were not normally distributed and were heavy-tailed (with kurtosis less than 3), so nonparametric tests were used for analysis of evenness differences across stages. Evenness was significantly different (P=0.0354) between stages 30 and 80 in the ARISA datatype in nonparametric pairwise tests (see Supplementary Results section 1.b for details), with the evenness value lower in the latest stage (i.e. a few OTU are more abundant in the stage 80 samples, compared to a more similar set of abundances in the first samples), though overall the evenness averages are not significantly different from the overall mean across all the sample days ( Supplementary Results section b). This pairwise-significant evenness difference in important OTU abundance ranges in the earliest and latest sample dates is also visible in a SIMPER analysis of ARISA OTU, with the top contributions to the day 30 OTU more evenly spread (for day 30: OTUs with 31, 19, 11, 11, 8, 6, 4 percent) as compared to day 80 (with OTUs of 50, 19, 10, 9, 5 percent, Supplementary Table 4).

SEM community composition and beta-diversity. In contrast to the ARISA datatype, in the SEM-based size-class unit (SCU) datatype the beta-diversity significantly decreased in the day 41 samples relative to the day 48 samples (Permdisp P=0.004, Table 2, details in Supplementary Table 3). This beta-diversity difference confounds Permanova tests, so we cannot further compare SEM community compositions across stages with that method. The significant SEM beta-diversity community change tells us that there is an important successional process at this specific stage relative to the earlier and later samples. Comparison of the Pielou’s J’ evenness values in the SCU datatype using nonparametric tests indicated that day 72 had significantly lower evenness of size-classes relative to the overall mean ( Supplementary Results section 2), with pairwise tests showing that day 72 was significantly lower in evenness than day 80 ( Supplementary Results 2.b). Pairwise comparisons among all other days were not significantly different.

Table 2. SEM Permdisp Community Structure Comparisons.

Significance groups with different letters are significantly different by pairwise test with a Bonferroni multiple-test-correction threshold of P=0.001315.

Stage (day after planting)	Average centroid	Significance group
30	61.789	A,B
41	55.974	A
48	64.629	B
58	63.854	A,B
62	62.28	A,B
72	61.929	A,B
80	61.892	A,B

Figure 3. Non-parametric analysis of means of microbial OTU number (richness).

The pale blue area is the confidence interval (UDL and LDL) above and below the overall mean. Rank-means for each stage are points at the tip of the drop-down bars. Significant means (adjusted P<0.05) have red line tips, while non-significant means within the confidence interval have green tips. A) ARISA mean differences with transformed ranks. B) SEM mean difference of transformed ranks of taxa.

Presence-only comparisons

We compared the number of OTU present in each sample to the overall mean using nonparametric comparisons suitable for non-normal, heavy-tailed distributions ( Figure 3). The number of OTU – the richness – increases significantly in the ARISA datatype at day 62 (this mean is above the confidence interval and is thus indicated with a red dot in Figure 3 a), with earlier and later samples showing no significant difference from the overall mean richness (means indicated with green dots) and no pairwise differences in means compared to each other with nonparametric tests ( Supplementary Results section 1a). In contrast to the ARISA datatype pattern, the number of different SCU present in the size-class datatype is significantly lower at day 48 ( Figure 3b), while all other stages are within the confidence limit for the overall mean richness. Nonparametric pairwise comparisons of SCU indicated that day 48 was significantly lower than all later-stage samples and that days 72 and 80 had significantly higher SCU means present than day 30 ( Supplementary Results section 2a).

Core community members

For the ARISA datatype, there were three OTU found in all stages of community development; these were OTU9, OTU22, and OTU24, which comprised 1.5 percent of the total number of OTU. For the SEM datatype, there were 13 SCU (5.7 percent of the total); a visual display of the proportions of all 13 across the seven sample dates is shown in Supplementary Figure 3.

Co-occurrence pattern

To further examine community assembly processes we looked at trends in the co-occurrence of pairs of species in both phyllosphere datatypes. We examined these patterns in our data to identify ways that community composition could differ. Positive co-occurrence is defined as species occurring together more often than random (and thus having potential mutualistic interactions), while negative co-occurrence values indicate that the species pair is found together less often than the random expectation (potential antagonistic interactions). Co-occurrence patterns changed across plant stages in both ARISA ( Figure 4) and SEM ( Figure 5). Participation of certain OTU and SCU in varying interactions across multiple stages is indicated by the red color intensity and number of overall edge interactions in the DyNet networks ( Figure 4b and Figure 5b). In the ARISA datatype, highly variable co-occurring OTU (intense red color in Figure 4b) had both positive and negative interactions with other community members. In contrast, the SEM datatype exhibited more positive interactions connecting highly dynamic nodes (intense red nodes in Figure 5b). There were more positive than negative significant co-occurrences in both the ARISA and SEM datatypes ( Figure 4 and Figure 5, edge arrows with points). The number of co-occurrences increased with increasing plant growth time in each datatype, though not at the same stage; as seen in Figure 4, ARISA had the largest number of network nodes at day62 and for SEM ( Figure 5) the largest number of nodes was day72. The overall proportion of co-occurring pairs was similar in each datatype ( Table 3), with a minority of microbial types found to have significant co-occurrence partners. Specific OTU and SCU identifiers are individually displayed in stacked plots that show details about which OTU occurred in more than one pair within Supplementary Figure 1 and Supplementary Figure 2.

Figure 4. Network visualization of significantly co-occurring ribotype OTU pairs at each plant stage, with the identifier number for each OTU shown with the label “a” plus the number.

If there were no significant co-occurring OTU, then that sample stage was not shown. Negative interactions between OTU nodes are shown as edges with a barred arrow tip and positive correlations between OTU are indicated with pointed arrows. a) Significant co-occurring OTU from each sampled stage. b) Results of DyNet analysis of co-occurring OTU most changed across stages. Only significant co-occurring OTU (as shown in part a) were analyzed using DyNet. Intensity of the red color of nodes indicates amount of change across stages, with darker color indicating the most ‘rewiring’ across stages.

Figure 5. Network visualization of significantly co-occurring SCU pairs at each plant stage, with the identifier number for each SCU shown with the node label “s” plus the number.

If there were no significant co-occurring SCU, then that day was not shown. Negative interactions between SCU nodes are shown as edges with a barred arrow tip and positive correlations between SCU are indicated with pointed arrows. a) Significant co-occurring SCU from each sampled stage. b) Results of DyNet analysis of co-occurring SCU most changed across stages. Only significant co-occurring SCU (as shown in part a) were analyzed using DyNet. Intensity of the red color of nodes indicates amount of change across stages, with darker color indicating the most ‘rewiring’.

Table 3. Proportion of Co-occurring Pairs within the ARISA and SEM Community Measures.

Datatype	Total number of OTU or SCU	Percent that significantly co-occur at least once	Percent with no interaction
ARISA OTU	195	14	93
SEM OTU	229	28	88

ARISA and SEM taxonomic abundance data

Detailed description of column headers, row labels and cell content for Dataset 1 and Supplementary materials

Discussion

Our experimental design was selected to avoid confounding by differential exposure to recent weather events. This design, which has recently been recommended by Shade et al.(2013) and by Williams et al.(2013), allowed us to examine the contributions to community assembly due to dispersal, filtering and niche-related processes. Comparisons of presence-absence patterns to presence-only richness within a datatype are complementary analyses, as the presence-absence matrix allows for overall tests of community difference (Permanova-based, 33) and tests for beta-diversity difference as difference in homogeneity of variance (Permdisp, 34). Presence-only lists of taxa or size-classes allow examination of the richness dimension of the microbial communities, abstracted away from other structural differences.

Are two datatypes better than one? We measured a morphological trait (microbe size) and a relatedness (ribosomal length difference) metric, and the two datatypes were not especially similar at specific time points, making them complementary rather than redundant; we infer that predictive power from one datatype to the other would be low. The differences between datatypes suggest that size-classes may include different taxa, in other words our size-classes may not have a phylogenetic signal. Some recent work in phytoplankton and soil microbial communities contrasts with our results, showing that traits can be conserved across taxonomic groupings ^{26,
40}; generally conservation of traits and evolutionary relatedness is likely to depend on the genetic architecture of the phenotypic trait ⁴¹. The relationship between metabolic and functional traits and phylogenetic signal is as yet unclear for microbial communities that have substantial levels of gene transfer ^{42,
43}. Complementary information provided by two (or more) datatypes is thus recommended for future phyllosphere experiments.

Misclassification of ribosomal spacer lengths to different taxa when they are in fact due to population differences, and misclassification of large bacterial cells in a size-class separate from their smaller offspring, would have consequences for trends across samples if the misclassification probability were different at the different stages of community development. In our study, the underlying assumption is that our probability of misclassification is the same across stages. For size-classes, misclassification could increase with time due to growth of cells, resulting in a continuous increase in the number of SCU over time. We do not see this trend in our data, so from our study there is no evidence of confounding misclassification. We cannot rule out misclassification that occurs only at specific stages using our experimental design. We suggest that sample processing for future experiments include inline size separation so that size-classes, ribosomal data, and molecular/metabolic data are collected as true multivariate data in each sample. Image information such as our scanning electron micrographs or fluorescence-based microbial images could also be analyzed for patterns that could illuminate resource use traits and specific cell-shape-based mutualisms ⁴⁴.

Trends in richness, co-occurrence trends, and beta-diversity changes together define overall community assembly mechanistic contributions, with habitat filtering characterized by unchanged richness, unchanged co-occurrence frequency, and decreasing beta-diversity. Niche development would be observed as a combination of increasing richness, increasing co-occurrence, and lowered beta-diversity. Finally, a neutral pattern – with no change in immigration probability over time – would exhibit trends toward increasing richness, unchanged co-occurrence, and unchanged beta-diversity. Comparison of recent phyllosphere analyses and our work using this framework emphasizes the role of niche development as a mechanism of community assembly.

First we compare our results with recent work in apple flower microbial communities, where beta-diversity was initially at relatively high levels, decreased, then increased over the remainder of the season; richness (measured as phylogenetic diversity in the apple flower work) increased then leveled off, and finally co-occurrences showed early increase then decline and rebound toward higher levels over the season ²⁰. To summarize, early in the apple flowering season the diversity patterns support niche development mechanisms, and later in the season both habitat filtering and niche development patterns were seen. In our study of maize leaf microbial epiphytes we also found that there was no overall temporal trend. Our leaf ARISA pattern of diversity change was not especially similar to the apple flower changes; for example, we found no significant change in beta-diversity, while beta-diversity declined then rose over the later part of the season in apple flowers. Our ARISA richness and the apple flower phylogenetic diversity also differed, with apple flower phylogenetic diversity increasing then leveling off, while on maize leaves richness increased only in one later growth stage. The evenness in both our ARISA and apple flower phylogenetic diversity varied; in apple flower the ‘mid’ stage had a lower evenness while in our ARISA dataset the day 80 evenness was lower than in day 30. Co-occurrences in our study increased in mid-season, as did the apple flower communities in the early flower stage. Overall, the early season patterns in the across-season apple flower microbial community analysis were the most similar in community assembly to our maize leaf staggered-stage analysis. In addition to the apple flower study, in a recent study of three different crop leaf phyllosphere seasonal patterns, Copeland et al. ²¹ found that alpha-diversity increased then decreased across the season, with more shared species later in the season. In all three crops, beta-diversity (as measured as Unifrac distance for leaf-specific microbes) decreased across season, though not evenly; the beta-diversity decrease was more prominent later in season and at one time-point was higher during the later phase. Co-occurrences were not analyzed in this study, so these crop study results are consistent with niche and/or habitat mechanisms.

Core microbial communities in crops are those that are either recruited early and retained, or continually input across the season. These are often found by presence ⁴⁵, as in our analysis of ARISA and SEM core microbial community members. We found a relatively small proportion of core OTU and SCU; other phyllosphere studies such as on apple flowers ²⁰ found more (for apple flowers the generalist category was 14 percent, for example). Our core set should be considered as a first description, as we did not consider absences or compare our core set to air or soil samples. Instead, our core set focused on microbial community members that colonized leaves of all developmental stages.

Conclusions

Our ARISA taxonomic diversity patterns and microbial cell size trait datatypes exhibited different specific patterns of diversity, but there was similarity in trend overall, as both datatypes support niche process contributions to community assembly. Niche (deterministic) processes have recently been found to be important in soil microbial communities ^{46,
47}. In lakes, microbial generalists found at multiple spatially separated sites had a better fit to neutral stochastic models, with specialists exhibiting a more deterministic pattern ⁴⁸. As leaf surface communities are quite different than air or soil communities ⁹, and the number of core species in our inventory is low, we suggest that many of the maize leaf microbial community member that we detect are specialists. Our staggered experimental design highlighted the role of microbe-microbe interactions in leaf microbial community development. We thus suggest that future experiments use staggered experimental designs and multivariate metabolic, genic, and size-trait data collection.

Data availability

The data referenced by this article are under copyright with the following copyright statement: Copyright: ï¿½ 2017 Manching HC et al.

Data associated with the article are available under the terms of the Creative Commons Zero "No rights reserved" data waiver (CC0 1.0 Public domain dedication).

Dataset 1: ARISA and SEM taxonomic abundance data. DOI, 10.5256/f1000research.12490.d177516 ⁴⁹

Dataset 2: Detailed description of column headers, row labels and cell content for Dataset 1 and Supplementary materials. DOI, 10.5256/f1000research.12490.d177517 ⁵⁰

Files in the CyVerse/iPlant Data Commons repository are available from http://datacommons.cyverse.org/browse/iplant/home/shared/commons_repo/curated/Stapleton_MaizeLeafMicrobesByLeafAge_2015, with DOI, 10.7946/P2WC77. This repository includes the following individual files:

SupplementarySEMfiles Folder

binnedSEMSizes.csv

readmeSEMfolderfilemetadata.txt

SEMImages Folder and Subfolders

SEMoverlays Folder

SEMribosortinputfiles Folder

SEMRiboSortoutputfiles Folder

SupplementaryARISAfiles Folder

ARISAfilenamemetadata.csv

ARISARibosortinput Folder

ARISARibosortoutput Folder

ARISAruns Folder and SubFolders

readmeARISAfolderandfilemetadata.txt

Source data files formatted for CytoScape visualization are available at https://github.com/AEStapleton/plantstage.

Acknowledgements

We are grateful to A. Ulloa-Bacillio, UNCW and J. Geyer, UNCW for assistance with SEM image analysis, to the UNCW microscopy facility for training and equipment support, and to the staff of the Central Crops research station for experimental field management.

Supplementary material

Supplementary Figure S1: a-n Arisa Plots of Signifcant Co-occurring OTU Pairs.