Crown closure affects endophytic leaf mycobiome compositional dynamics over time in Pseudotsuga menziesii var. menziesii

Old-growth


Introduction
Endophytic fungi live a significant part of their life cycles within their hosts without producing symptoms of disease (Saikkonen et al., 1998). Those that infect plant leaves are generally horizontally transferred between host plants, comprise a diverse taxonomic and phylogenetic assemblage of fungi, and can modulate foliar disease severity (Busby et al., 2016), stomatal conductance (Arnold and Engelbrecht, 2007), and nutrient uptake (Christian et al., 2019). Given the functional importance of foliar fungal communities for regulating plant responses to biotic and abiotic stresses, elucidating the factors controlling their composition will have important implications for both basic and applied plant sciences.
Communities of foliar endophytic fungi vary along environmental gradients at landscape (Zimmerman and Vitousek, 2012) and regional (Barge et al., 2019) scales, yet there are few studies examining patterns of fungal community variation along the vertical axis of individual host plants (Osono and Mori, 2004;Harrison et al., 2016;Izuno et al., 2016;Oono et al., 2017). While this omission may be inconsequential for annuals and plants of modest stature, the crowns of tall, old-growth conifer species like Pseudotsuga menziesii var. menziesii support vertically stratified ecological communities in which compositions change with height in crown (McCune, 1993;Sillett and Rambo, 2000). And although P. menziesii crowns were among the first rigorously sampled for endophytic needle fungi (Sherwood and Carroll, 1974), further characterization of these endophyte communities with modern sequencing techniques is needed.
Tree crowns span multiple microclimatic gradients. UV-B and photosynthetically active radiation are abruptly attenuated by foliage in the mid-crown (Parker, 1997), while the upper crown is only intermittently wet, with generally warmer and less stable air temperatures than the lower portions of the crown that are buffered by vegetation (Heffernan, 2017). Further, unshaded foliage can become nearly 7 • C warmer than surrounding air temperatures, while leaf temperatures of shaded foliage are relatively stable (Doughty and Goulden, 2008). These and other microclimatic factors can also affect fungal colonization potential in other contexts (Carisse et al., 2000;Unterseher and Tal, 2006;Norros et al., 2015), and may also influence foliar mycobiome composition by acting as environmental filters dependent on the degree of canopy openness (Gilbert et al., 2007). Additionally, foliage at different heights in the crown possesses different morphological, physiological, and spectral properties (Ishii et al., 2002;Gebauer et al., 2015;Schweiger et al., 2020), potentially allowing host-mediated, within-crown biotic filtering of needle communities. Because P. menziesii canopy foliage can be long-lived, community assembly outcomes for the canopy mycobiome at a point in time may have lasting impacts on future community structure and assembly as foliage is exposed to infection over time (Bernstein and Carroll, 1977;Arnold and Herre, 2003;Osono, 2008).
The potential for light detection and ranging (LiDAR) point cloud data to account for variation in the diversity and composition of terrestrial fungal communities (Peura et al., 2016;Thers et al., 2017;Valdez et al., 2021) and to characterize the habitat preferences of crown-dwelling organisms (Barnes et al., 2016) has recently been demonstrated. These and other approaches rely on LiDAR technologies, which associate several emitted and detected laser pulse echoes ("returns") with accurately recorded positions. Further processing yields three-dimensional "point clouds'', with each point in space representing a position where the laser pulse reflects off vegetation or other surfaces.
Researchers can then obtain accurate estimates of canopy height, cover, and structure to inform their analyses of communities (Lefsky et al., 2002). LiDAR-derived metrics also correlate with specific aspects of the microclimate such as near-surface air temperatures in old-growth P. menziesii forests (George et al., 2015;Frey et al., 2016;Davis et al., 2019). Thus, combining mycobiome data with a co-localized measurement of the crown microenvironment may help to elucidate how endophyte communities vary along tree stems and among trees. Here, point cloud data derived from airborne LiDAR and high-throughput amplicon sequencing data were leveraged to test the hypothesis that the composition and diversity of fungal needle endophyte communities are correlated with crown exposure. In addition to examining the height above ground as a compositionally relevant factor, our definition of exposure considered an airborne LiDAR-derived measure of crown closure, which quantifies the amount of intervening foliage (measured as "points" in a point cloud) between the sky and a specified height in the tree. Altogether, positions in the crown that are further from the  Spies, 2016. Relationship between height in crown and crown closure (B). Grey lines show data from point clouds associated with 5828 trees randomly sampled across HJA (each with randomly sampled crown heights), while colored lines and points represent point clouds associated with trees that were climbed and sampled as part of this study. Point cloud cylinders are centered on the treetop of a focal tree, but include all point returns within a 10 m radius (C). Although the same vertical height can be sampled in two different trees, these heights can be associated with different crown closure indices. ground and less obscured by foliage are more exposed; positions closer to the ground and more obscured by foliage (more closed) are less exposed. We further explored the relevance of exposure as a concept and hypothesized that more-or less-exposed crowns undergo different shifts in community structure with each age class of needle, potentially reflecting differences in needle endophyte community assembly. We report that only closure uniquely contributed to an ecologically informed metric of exposure, accounting for variation in community composition and diversity, and exposure groups showed differing patterns in community structure over the needle ages sampled.

Sampling and surface sterilization
Six trees (BIRD348, CC, LATREEN, PC02, PC12, and PC17) at the H. J. Andrews Experimental Forest (HJA, Fig. 1) were sampled in late winter 2016-2017, and two other trees (DISCOVERY and DT NEIGHBOR) were sampled in winter 2017-2018. HJA is located on the west side of the Cascade Range of central Oregon; its 6400 ha range from 410 to 1630 m in elevation, and canopies can exceed 93 m in height. Using single-rope climbing techniques (Video S1), where possible, accessible branches from opposing aspects were sampled every 1-5 m in each tree, and the height of each sample collection above ground level was noted. Collection of branches began with the first available branch closest to the ground, and continued as high as could be safely climbed. Due to safety concerns, climbing did not reach the tops of tree crowns. Branches were transported in coolers with ice blocks and then later refrigerated at 4 • C until they could be processed (within 72 h). Needles that experienced one (A1), two (A2), three (A3), and four (A4) growing seasons were sampled from each branch and represented four unique needle age classes. To target endophytic fungi, needles were surface sterilized via immersion in 70% ethanol for 1 min, 4 min in 5% sodium hypochlorite, 1 min again in 70% ethanol, followed by a sterile distilled water rinse before being frozen at − 80 • C (modified from Thomas et al., 2016). Needles from different aspects at a given height on each tree were pooled prior to 72 h of lyophilization. Each sample consisted of a unique tree-by-height-by-age combination (71 heights x 4 ages = 284 samples).

DNA extraction, Illumina MiSeq sequencing, and sequence processing
DNA was extracted using OPS Diagnostics' 96 Well Synergy™ Plant Homogenization plates (Lebanon, New Jersey, USA). Each well received five needles of a single age class, sourced from each sampled height; empty wells were processed per the kit instructions and served as extraction controls. To reduce PCR inhibition issues, DNA extractions were cleaned with 1.2X volumes of Mag-Bind® TotalPure NGS magnetic beads (Omega Bio-Tek, Norcross, Georgia, USA). The ITS2 region was amplified in a primary PCR using the 5.8S-Fun and ITS4-Fun fungalspecific primers (Taylor et al., 2016) Each primer consisted of sequences that amplified the ITS2 region, followed by 3-6 bp length heterogeneity spacers, and subsequently terminated with an Illumina adapter sequence. 25 μl PCRs were carried out with MyTaq™ HS Red Mix (Meridian Bioscience, Inc., Cincinnati, Ohio, USA), 0.4 μM of each primer, 5 μl of template, and applying the following cycling parameters: 3 min denaturation at 95 • C, 28 cycles of 95 • C (30 s), 58 • C (30 s), 72 • C (30 s), followed by a 2 min extension at 72 • C. Using a secondary PCR, sample barcodes and Illumina flow cell adapters were then incorporated into these PCR products after amplification was visually confirmed via agarose gel electrophoresis. 25 μl PCRs were carried out with MyTaq™ HS Red Mix, 0.5 μM of each barcode primer, 1 μl of template, and with the following cycling parameters: 3 min at 95 • C, eight cycles of 95 • C (30 s), 55 • C (30 s), 72 • C (30 s), followed by a 2 min extension at 72 • C. This PCR product was purified and normalized (to 2.5-3.0 ng μl − 1 ) using Charm Biotech's Just-a-Plate™ 96 Well Purification and Normalization plates (Cape Girardeau, Missouri, USA) and submitted the library to the Center for Genome Research and Biocomputing (Oregon State University, Corvallis, Oregon) for 2 × 300 paired end sequencing on the Illumina MiSeq platform (Reagent Kit v3, San Diego, California, USA).
Reads were demultiplexed with Pheniqs (Galanti et al., 2017) and adapters were trimmed using Cutadapt (Martin, 2011) and SeqPurge (Sturm et al., 2016). Reads were then denoised, filtered of chimeras, and merged using the DADA2 R package (Callahan et al., 2016), producing amplicon sequence variants (ASVs) using the pseudo-pooling method. ITS2 sequences were extracted with ITSx (Bengtsson-Palme et al., 2013) and were aligned to the 2020-02-04 UNITE eukaryote release (singletons included, Abarenkov et al., 2020a) with VSEARCH (Rognes et al., 2016). If at least half of each ASV's global alignments (at > 0.5 identity) were fungal, the ASV was retained as a putative fungal sequence. The LULU R package, implementing the LULU algorithm, was used to collapse potentially infragenomic or artifactual ASVs into OTUs (Frøslev et al., 2017), with minimum match = 97% similarity, and minimum relative occurrence = 95%. Although not all ASVs were affected by this approach, putative taxa will be referred to as OTUs. After this step, OTUs for which >1% of non-control reads were found in the extraction and PCR controls were identified as contaminants and removed. Taxonomy was then assigned to the remaining subset of OTUs using the fungal version of the 2020-02-04 UNITE release (without singletons, Abarenkov et al., 2020b) and the RDP classifier (Wang et al., 2007), implemented with DADA2. This release was supplemented with full ITS sequences belonging to fungal taxa known to associate with P. menziesii or other conifers, but not included in the UNITE release. The full ITS sequence of P. menziesii was also included in the release, allowing for host amplification to be identified and removed. The PERFect R package was used to remove spurious OTUs by identifying those OTUs which made insignificant contributions to the total covariance of the dataset (Smirnova et al., 2019); the permutation method with default significance thresholds (retaining OTUS with P < 0.1, with P values calculated using 1000 permutations) was applied. Altogether, LULU conditional OTU curation and PERFect filtering were employed as a means to retain as much sequence resolution as possible while also attempting to minimize the effects of infragenomic taxon splitting and artifactual sequence generation.

Airborne LiDAR data-processing
Crown closure values associated with each sample height were obtained by processing light detection and ranging point cloud data that had been collected for a larger survey of the McKenzie River study area. Between June 07, 2016 and June 21, 2016, Quantum Spatial (Portland, Oregon, now NV5 Geospatial, Hollywood, Florida, USA), commissioned by the Oregon LiDAR Consortium, obtained discrete-return airborne laser scans of the HJA area from a Leica ALS80 (Leica Geosystems, St. Gallen, Switzerland) mounted on a Cessna Grand Caravan (Textron Aviation Inc., Wichita, Kansas, USA), flying 1500 m above ground level. Across the entire 2016 McKenzie River survey area, an average pulse density of 12.64 m -2 (SD = 4.867, min = 6.70 pulses m − 2 , max = 65.66 pulses m − 2 ) was verified by the Oregon Department of Geology and Mineral Industries. Using 1,723 ground survey points, it was also found that the aerial survey had a horizontal accuracy of 0.37 m root mean square error (RMSE) and a vertical accuracy of 0.053 m RMSE. Point cloud data from this survey (OCM Partners, 2020) were downloaded from the NOAA Data Access Viewer as a projection to UTM Zone 10 (horizontal datum = NAD83 (2011), vertical datum = NAVD88, geoid = GEOID18), with meters as horizontal and vertical units.
Point clouds were processed with the lidR package (Roussel et al., 2020). Point clouds were purged of point returns with GPS time errors and intensity values of zero. For each tree, cylindrical volumes with radius 10 m (and centered on tree positions) were clipped from these point clouds. This was the minimum radius that contained the horizontal breadth of each crown for each sampled tree. Prior to clipping, tree positions were refined, starting with GPS coordinates on record for each tree. The expertise of the HJA forest director (Dr. Mark Schulze) was used to verify final tree positions, which were centered on each treetop. For each cylinder containing the point cloud of each tree, the bare earth elevation under each tree was estimated via k-nearest neighbor inverse-distance weighting (k-nearest neighbor points = 10) using all available ground points. This bare earth elevation was added to the heights at which needle material was sampled in each tree. From here, a crown closure index was defined for each sampled tree height as the number of point returns above the sampled height, divided by the total number of point returns in each cylinder (Lovell et al., 2003). These closure values were then multiplied by the median height of elevation-normalized, non-ground point returns (Z med , Fig. 1C) to account for the fact that a point cloud cylinder with a greater amount of tall vegetation will generally yield more, intervening point returns than a cylinder containing shorter vegetation (assuming equal scan coverage).

Analysis
To determine how needle community composition varied across trees, needle ages, and crown properties, features for each sample were relativized by sample sequencing depth (i.e., converted to relative abundance), and a generalized log transformation was applied to dampen the compositional influence of highly abundant OTUs. Bray-Curtis dissimilarity matrices were generated for use in downstream analyses. All subsequent analyses and ordinations were performed using functions from the vegan R package (Oksanen et al., 2020) unless stated otherwise. Dissimilarities associated with the full dataset were ordinated with non-metric multidimensional scalings (NMDS), using a two-dimensional solution. NMDS significance was evaluated by performing a Monte Carlo test (with 999 permutations) on the results using R scripts (calling vegan functions) modified from those originally created by Dr. Kevin McGarigal (pers. comm.). Four samples which prevented the NMDS ordination from converging were purged from the first needle age class. In trimming each needle age subset down to only the samples that shared a tree height with every other subset prior to analysis, samples were removed from the other three needle age classes if one needle age class was missing. This resulted in the retention of 64 samples for each needle age class subset (256 samples total). Permutation analyses of variance (PERMANOVA; Anderson, 2001) were performed with the adonis2 function on the full dataset to determine whether needle compositions differed by tree and needle age, and the homogeneity of variances among levels of each factor was tested using the betadisper (implementing PERMDISP2; Anderson, 2006) and permutest functions to assess whether differences in composition were influenced by differences in group dispersion. For these tests, 999 random permutations were used to generate null pseudo-F ratio distributions, designating source tree ("tree") as a permutation stratum. Analyses were then followed by a PERMANOVA on the full dataset testing the effects of height and closure after accounting for needle age class, employing 999 permutations to generate P values and with "tree" as permutation stratum. Distance-based redundancy analyses (dbRDA) were constrained on crown variables, and PERMANOVAs were performed on each needle age subset to determine whether compositional variation accounted for by both crown variables changed with each needle age class.
Crown terms (i.e., among height and closure) accounting for unique compositional variation in the marginal test on the full dataset were retained to generate exposure groups. To ensure that group assignment procedure was not influenced by the small sampling of trees in our dataset, the point clouds of 5828 randomly sampled trees across HJA were also characterized by calculating closure at 10 random heights in each tree, using the method described above (Fig. 1B). Random sampling was restricted to trees attaining the minimum treetop height (35 m) and Z med (22.67 m) values observed among the climbed trees in our dataset, and random height sampling was also restricted to the range observed for climbed trees (lowest height = 10 m, greatest height = 55 m). Closure values calculated at each sampled height were clustered into two groups using the Hartigan-Wong k-means clustering algorithm (Hartigan and Wong, 1979) with the maximum number of iterations and the number of random starts both set to 999. Samples from heights associated with the lower mean cluster were designated as belonging to the "open" exposure group, and samples associated with the higher mean cluster were designated as belonging to the "closed" exposure group. This was followed by indicator species analysis on the full dataset (Dufrene and Legendre, 1998; as implemented with the labdsv package by Roberts, 2019), which was used to identify taxa that occurred more frequently and with greater relative abundance than expected by chance in either open or closed exposure groups. Only OTUs occurring in at least 10% of samples were subjected to indicator species analysis. P-values were adjusted to control for the false discovery rate of indicator taxa using a Benjamini-Hochberg correction (Benjamini and Hochberg, 1995).
To determine whether alpha diversity varied with tree, needle age, and closure, and height richness (Hill number q = 0) and the exponent of Shannon entropy (q = 1) were estimated for each sample using the iNEXT package (Hsieh et al., 2016). Both metrics were estimated at a depth of 1000 reads, and variation in richness and log 2 -transformed Shannon entropy (i.e, the Shannon index of diversity) was tested among source tree and needle age (and their interaction) using analysis of variance (ANOVA) F-tests (Fox and Weisberg 2019). Assumptions of normality and homoscedasticity were investigated. In the event of putative heteroscedasticity, a White standard error adjustment was applied in an attempt to correct for its influence in the ANOVA (Long and Ervin, 2000). Results of F-tests were investigated further with Tukey honest significant differences post-hoc tests (alpha = 0.05). Alpha diversity means were bootstrapped and confidence intervals were generated with the bca function of the coxed R package (Kropko and Harden, 2020). Linear mixed-effect models were created using the nlme package (Pinheiro et al., 2021) to test the correlations of alpha diversity metrics with closure and height, setting source tree as a random effect. Reduced models containing either height or closure were compared to each other, a full model containing both crown variables, and a base model that only contained the tree random effect according to their corrected Akaike information criterion (AICc) and Bayesian information criterion (BIC) values using the AICcmodavg R package (Mazerolle, 2020). Reduced model P-values were obtained by performing a likelihood-ratio test comparing the reduced model to the base model. Models included needle age class as a fixed effect term if the term was significant in earlier ANOVAs. Pseudo-R 2 values of fixed effects were calculated using the piecewiseSEM package (Lefcheck, 2016).
Open and closed exposure groups were examined for changes in community structure over time by performing Mantel tests at each transition between needle age classes (A1 to A2, A2 to A3, A3 to A4). Rows and columns associated with samples belonging to either the open or closed exposure groups were extracted from the Bray-Curtis distance matrix associated with each needle age class and used as input matrices for Mantel tests, each using 999 permutations, setting source tree as a permutation stratum. Spearman ranked Mantel statistics were also bootstrapped with 999 permutations. With each permutation, tree heights shared by the needle age classes being compared were randomly sampled with replacement and the Bray-Curtis dissimilarity matrix was generated among the corresponding samples belonging to each needle age class. As above, samples belonging to either the open or closed exposure groups were extracted and the ranked Spearman correlation between matrices was found. Means were calculated from these bootstrapped correlation coefficients, and confidence intervals were generated as above. To compare the relative abundances of dominant OTUs across needle age transitions, paired Wilcoxon signed rank tests were performed in R. When the null hypothesis was rejected for a pair of age classes (alpha = 0.05), a two-sample Wilcoxon ranked sum test was performed comparing OTU relative abundances between open and closed exposure groups for the older needle age class. Exact P values (as opposed to normally approximated P values) were calculated when possible.

Reproducibility
All statistical analyses were conducted in the R v4.0.3 statistical computing environment (R Core Team 2020). The functions associated with the tidyverse R package (Wickham et al., 2019) featured heavily in many scripts and facilitated the manipulation of data structures. Handling of metabarcoding data was accomplished with the phyloseq package (McMurdie and Holmes, 2013). Novel code and metadata associated with bioinformatics, point cloud processing, statistical analyses, and figures have been publicly archived (https://doi.org/10. 5281/zenodo.6360767). Demultiplexed sequencing data are deposited in the NCBI Sequence Read Archive (BioProject ID PRJNA748821, BioSample accession SAMN20345134). LiDAR point cloud data are available for download from the NOAA data portal (https://www.fisheri es.noaa.gov/inport/item/49949).

Overview
Before the removal of contaminants and non-fungal sequences, ITS

Fig. 2.
Taxonomic assignments of OTUs present in the top 95% of reads, pooled at the tree-level (A). OTUs are colored by genus. OTU abundance distribution across the eight trees sampled in this study (B). Occupancy-abundance curve, plotting the mean relative abundance of each OTU against the number of trees it was found to occupy (C). The y-axis of (C) was log10-transformed. Each point represents a different OTU. extraction, chimera removal, and LULU clustering, 1786 ASVs were identified across 284 non-control samples. After ITS extraction, there were 1316 ASVs and 1292 ASVs after chimera filtering. 138 ASVs were attributed to Pseudotsuga ITS sequences, totalling 561,671 reads. 1115 putatively fungal ASVs were retained after the VSEARCH filtering process. After 97% LULU clustering and contaminant removal, 899 OTUs were delineated, represented by 6,123,130 sequences over 282 samples. This was reduced to 218 OTUs, 5,001,904 sequences, and 256 samples after PERFect OTU filtering, contaminant removal, and the removal of samples that did form a full complement of needle age classes. OTU accumulation curves, generated with functions from the BiodiversityR package (Kindt and Coe 2005), indicated that each tree was adequately sampled (Fig. S1).

Table 1
PERMANOVA results displaying the amount of compositional variation accounted for by needle age and tree, and the marginal variation of height and crown closure. P values are the results of permutation tests of pseudo F-ratios using 999 iterations, setting "tree" as a stratum.   (Table S3).

Structure
Fungal community composition from heights in the open and closed exposure groups were differentially structured over needle age class transitions (Fig. 5). For closed exposures, community structure became increasingly correlated over the A1-A4 transitions. No significant structural correlation was detected at the A1 to A2 transition (P = 0.762), nor for the A2 to A3 transition (P = 0.204) but community structures were strongly correlated with each other at the A3 to A4 (P = 0.034, Spearman correlation ρ = 0.782). Conversely, the community structures of open exposures were not correlated at any transition stage (P ≥ 0.220). The relative abundance of N. gaeumannii (OTU.1) in all exposures only increased at the A2 to A3 transition (Wilcoxon signed rank test, V = 208, P < 0.001). OTU.1 relative abundances in closed exposures, however, remained lower than relative abundances in open exposures at A3 (Wilcoxon rank sum test, W = 815.5, P < 0.001) with slightly elevated proportions in A3 and A4 (Fig. 6A). Conversely, R. parkeri OTUs (OTU.2, OTU.3, and OTU.6) across all exposures increased in relative abundance from A1 to A2 (V = 406, P < 0.001), but larger relative abundances were achieved in closed exposures at A2 (W = 238, P < 0.001, Fig. 6B). These R. parkeri OTUs decreased in relative abundance at the A2 to A3 transition (V = 1649, P < 0.001), but closed exposures had higher relative abundances than open exposures (W = 206, P < 0.001) at A3.

Discussion
This study is the first of its kind to combine airborne LiDAR-derived metrics with fine-scale, within-crown foliage sampling of known needle age classes, and high-throughput amplicon sequencing. Additionally, forest-level LiDAR data were exploited to more objectively classify the crowns of climbed trees in the dataset. Altogether, these approaches showed that the diversity, composition, and structure of the endophytic needle mycobiome are correlated with the surrounding crown microenvironment in a way that depends on the age classes of needles. This study demonstrates the unique utility that remote sensing technologies like LiDAR can provide to studies of "hard-to-reach" ecological communities, while also allowing researchers to frame the sometimes limited environmental variation among their sampled sites in the context of variation observed at larger scales. Further, the high throughput amplicon sequencing approach taken here resulted in identifications of the same dominant taxa described in earlier, established work. Namely, Nothophaecryptopus gaeumannii (Stone et al., 2008), Rhabdocline parkeri (Sherwood-Pike et al., 1986), and Phyllosticta Fig. 4. Estimates of OTU richness or the Shannon index of diversity associated with needle age (A), crown closure index (B), and tree (C). Estimates were obtained for each sample via interpolation or extrapolation to a depth of 1000 reads. Violin plots illustrate the distribution of richness and Shannon index diversity estimates for each grouping. Letters for (A) and (C) indicate significant difference groups from Tukey HSD post-hoc comparisons (alpha = 0.05) among levels of each factor after performing ANOVA F-tests. Horizontal lines delineate the bootstrapped means of estimates, while vertical lines correspond to bootstrapped 95% confidence intervals (n = 999).

Crown closure, interacting with needle age, accounted for more community variation than height
Source tree accounted for the greatest amount of compositional variation among samples (Table 1, Fig. S2), likely due to the large number of OTUs specific to only one or a few trees (Fig. 2). These results are similar to those of Harrison et al., (2016) and Oono et al., (2017), where, despite spatial ranges differing between the three studies by an order of magnitude (~10 km for the current study vs ~100 km for Oono et al. vs ~700 km for Harrison et al.), source tree or site also accounted for a large portion of variation in the composition of the foliar endophytic mycobiome of conifers. Further, these and other (Izuno et al., 2016), earlier studies find compositional differences between sampled crown positions. Analogous to crown position, height measured in the present study was unable to account for unique compositional variation (Table 1), perhaps suggesting that previously reported vertical stratification effects are actually unmeasured crown closure effects, given that crown closure decreases monotonically with increasing height in crown. This is supported by culture-based and visual survey studies which showed that certain foliar fungal taxa appeared to demonstrate preferences for either qualitatively defined shaded or exposed foliage (Osono and Mori, 2004;Gilbert et al., 2007;Unterseher et al., 2007;Scholtysik et al., 2013), but this effect is not consistently observed (Taudière et al., 2018). However, other work examining terrestrial macrofungal communities has established correlations between airborne LiDAR-derived measurements of canopy cover, maximum canopy height aboveground, and compositional turnover (Thers et al., 2017). In any case, quantifying exposure beyond canopy position or height should facilitate comparisons between studies of crown-associated foliar fungi (Arnold and Herre, 2003).
Analyzing each needle age class separately showed that compositional signal attributable to crown closure was detected at all needle ages, but peaked at A2 and A3 (Table S2, Fig. 3). The closure signal was stronger overall, but constrained to one needle age class when metabarcoding data were analyzed according to relative abundances that were not log-transformed. This result is indicative of the extent to which endophytic fungal communities were dominated by a few common taxa, a condition also reported by Harrison et al., (2016). Regardless of whether rarer or more common taxa were given more weight, closure-mediated compositional signal continued to collapse after needle age class A3 (Fig. 3). Thus, this shift in community composition for older needle age classes marks when smaller-scale, intra-crown community dynamics might give way to larger-scale or non-closure-mediated assembly dynamics. Needle age class has also been found to account for compositional variation of foliar fungal communities in both metabarcoding (Würth et al., 2019) and culture-based studies (Osono, 2008), and this is likely due to higher incidences of infection observed for older needle age classes (Sherwood and Carroll, 1974;Bernstein and Carroll, 1977;Petrini and Carroll, 1981) resulting from longer exposure to inoculum sources (Arnold and Herre, 2003). However, differences in needle physiology and morphology also accrue with increasing needle age class (Porte and Loustau, 1998;Bernier et al., 2001;Apple et al., 2002;Ishii et al., 2002;Yan et al., 2012;Eimil-Fraga et al., 2015), potentially accounting for its influence on endophyte community composition.
The relevance of crown closure is also supported by the finding that OTU richness and diversity increased with increasing closure (Fig. 4B,  Fig. S3), with the model selection process explicitly excluding height. Again, these findings are similar to those reported from foliar fungal surveys and airborne LiDAR-supported macrofungal studies (Osono and Mori, 2004;Gilbert et al., 2007;Unterseher et al., 2007;Scholtysik et al., 2013;Thers et al., 2017), but evidence from metabarcoding studies is mixed and complicated by the inconsistent use of canopy position or light exposure groups between studies (Harrison et al., 2016;Izuno et al., 2016;Oono et al., 2017;Taudière et al., 2018). If compared with findings that composition and richness vary along larger-scale temperature and precipitation gradients (Giauque and Hawkes, 2016;Bowman and Arnold, 2021;Oita et al., 2021), the results observed in the current study might suggest the more immediate action of microclimate-mediated environmental filtering mechanisms that operate at smaller scales. However, this interpretation relies on the strength of the relationship between within-crown LiDAR and microclimatic metrics, which has yet to be sufficiently demonstrated. Further, physiological and morphological needle traits vary within the crowns of conifers in response to these same microclimatic gradients (Ishii et al., 2002(Ishii et al., , 2008Gebauer et al., 2015;Schweiger et al., 2020), and this variation may have the potential to act with even more immediacy in structuring conifer endophyte communities. In fact, correlations of crown closure with composition and diversity can still suggest a role for closure-mediated dispersal limitation in structuring these communities (Gilbert and Reynolds, 2005), modifying the effect of closure-mediated environmental filtering.

Community structure diverges via interactions between needle age and needle exposure
Further exploration of needle age class and exposure group showed that inter-age differences in community structure differed between open and closed exposure groups. Although no structural correlations were detected at the A1-A2 and A2-A3 transitions for both exposure groups, strongly correlated structure was detected for closed exposures at the A3-A4 transition, indicating that turnover proceeded in a way that preserved some of the structure associated with the previous needle age class (Fig. 5). Conversely, open exposures never demonstrated significant structural correlations, signifying dynamic shifts in one or many compositional elements between needle age classes. Although the relative abundance of open exposure indicator OTU.1 (identified as Nothophaeocryptopus gaeumannii; Table S3) increased at the A2 to A3 transition, this increase was less pronounced for needles sampled from closed exposures (Figs. 3 and 6A), suggestive of altered community assembly dynamics relative to needles from open exposures. However, without sampling even older needle age classes, it is not apparent that (1) endophyte communities in needles from closed exposures would ever fully converge on a composition defined by N. gaeumannii dominance, (2) that the community structure of needles from open exposures ever stabilizes, or (3) that closed exposures retain age-class-to-age-class community structure.
Overall, observations from this study align with what is known about the infection cycle of N. gaeumannii and the general pathology of the Swiss needle cast (SNC) disease it can cause in P. menziesii. Conidium production has never been observed. Instead, needles are initially infected solely via ascospores (Stone et al., 2008) dispersed by rain. Ascospore release generally coincides with P. menziesii bud burst, and older needles appear to resist ascospore infection to some degree, leading pathologists to conclude that infection primarily occurs on newly emerged foliage during the spring and summer. After entering via stomata, N. gaeumannii proceeds to extensively colonize intercellular spaces over the lifetime of the needle, perhaps assisted by periodic epiphytic growth and stomatal reentry. Thus, continued colonization after early infection likely accounted for the majority of the observed increases in N. gaeumannii relative abundance at each needle age class transition (Fig. 6). The extent of this advantage was apparent given that the increase in OTU richness from A1 to A4 (Fig. 4A) was not met with an increase in the Shannon index of diversity over these needle age classes (Table S4). This suggests that while other taxa are able to at least infect a needle over its lifetime, they are generally unable to colonize the needle as effectively as N. gaeumannii over this period. The greater relative abundance and incidence of N. gaeumannii OTU.1 in open exposures is corroborated by experimental findings that shade-treated P. menziesii seedlings showed reduced SNC severity relative to the full-sun treatments (Manter et al., 2005). This measure of severity was also correlated with higher winter temperatures measured over the course of the experiment. Results from this earlier experiment are also in accordance with the common observation that SNC disease severity is negatively correlated with depth in crown in natural and managed stands of P. menziesii (Hansen et al., 2000;Lan et al., 2019Lan et al., , 2022Ritokova et al., 2020). A greater abundance of N. gaeumannii in warm, open exposures could reflect morphological and/or physiological adaptations to this niche. For example, previous work has reported that cultures of N. gaeumannii can deposit red pigments into the growth medium, suggesting that N. gaeumannii, like Cercospora spp. also in the Mycosphaerellaceae, may produce cercosporin or cercosporin-like perylenequinones (Winton et al., 2007). Cercosporin, which is red in its active, oxidized state, was originally isolated from Cercospora kikuchii and since been shown to disrupt plant cell membranes via the Fig. 6. Relative abundances of OTU.1 (identified to Nothophaeocryptopus gaeumannii) (A) and OTU.2, OTU.3, and OTU.6 (identified to Rhabdocline parkeri) (B) in each exposure group and needle age class. photoactivation of singlet oxygen (Daub and Ehrenshaft, 2000). Similar to SNC patterns, the incidence and progression of several diseases caused by perylenequinone-producing fungi (Cercospora spp., Alternaria alternata) are also positively associated with increased light exposure (Daub et al., 2013). Further investigation of the most recent assembly of the N. gaeumannii genome (BioProject PRJNA212511, BioSample SAMN02254965) with the antiSMASH BGC prediction program (fun-giSMASH default parameters, Blin et al., 2021) revealed the presence of a putative cercosporin BGC resembling one found in Cercospora beticola with 31% similarity (Table S5). Additionally, a T1PKS melanin BGC was identified, showing 100% similarity with a cluster found in Bipolaris oryzae. These findings, taken alongside increasing relative abundance of N. gaeumannii with decreasing crown closure, suggest that N. gaeumannii may be adapted to a niche characterized by warm temperatures and/or pronounced light exposure.
In contrast to N. gaeumannii, dominant Rhabdocline parkeri OTUs (OTU.2, OTU.3) were identified as taxa associated with closed exposures (Table S3), which is also in accordance with the known biology of this endophyte. Conidia of R. parkeri infect needles of P. menziesii and remain within a single epidermal cell until the needle approaches senescence. Rain-dispersed conidia produced on abscised needles can infect living needles of all needle age classes during the wet winter months (Sherwood- Pike et al., 1986). Needles are soon overtaken by generalist saprotroph communities once abscised needles reach terrestrial leaf litter (Stone, 1987), but R. parkeri has been shown to dominate the endosphere of abscised needles shaken directly from ten year-old P. menziesii crowns (Gonen, 2020). These abscised needles lodged in the crown then serve as inoculum sources for new infections in the crown. To account for the results of our indicator species analysis, it is hypothesized that the crowns above closed exposures intercept and/or retain more abscised needles on branches compared to portions above open exposures. P. menziesii canopy soils are derived from bark and decomposing needles deposited on branches (Pike et al., 1975), and these soils might then serve as reservoirs for R. parkeri and potentially other endophytic fungi (Looby et al., 2020). The capacity for a tree crown to retain needle litter reservoirs close to living foliage might also be adaptive if fungi sourced from this litter can reduce pest or pathogen damage, as has been described for R. parkeri (Carroll, 1988). This reservoir effect could also account for the finding that OTU richness and diversity increased with increasing crown closure ( Fig. 4B; Fig. S3). However, closed exposures could also be associated with more moist, thermostable, and otherwise accommodating microhabitats that are favorable to infection by R. parkeri and/or fungi in general, potentially accounting for or contributing to the indicator status of R. parkeri OTUs and the greater taxonomic diversity associated with closed exposures (Unterseher and Tal, 2006).

Declaration of competing interest
The authors declare no conflicts of interest.