Aponte_Bolivar_mim2_manuscript.qmd

---
title: "Foliar fungal symbionts in sympatric yellow monkeyflowers along elevation gradients in the Sierra Nevada"
date: last-modified
date-format: "MMMM D, YYYY"
bibliography: references.yaml
#csl: apa.csl
csl: american-journal-of-botany.csl

author:
  - name: Bolívar Aponte Rolón
    orcid: 0000-0002-2544-4551
    email: baponterolon@tulane.edu
    corresponding: true
    affiliations:
      - ref: tu
        name: Tulane University
        department: Ecology and Evolutionary Biology
        city: New Orleans
        state: LA
        country: US
        url: www.tulane.edu
  - name: Kathleen G. Ferris
    email: kferris@tulane.edu
    corresponding: false
    affiliations:
      - ref: tu
  - name: Sunshine A. Van Bael
    orcid: 0000-0001-7317-3533
    email: svanbael@tulane.edu
    corresponding: false
    affiliation:
      - ref: tu

affiliations:
  - id: tu
    name: Department of Ecology and Evolutionary Biology, Tulane University, 6823 St. Charles Avenue, New Orleans, LA 70118
    
filters:
  - authors-block

format:
  docx:
    reference-doc: custom-reference-doc.docx
    number-sections: true
    tbl-colwidths: true
    citation-package: biblatex
    
crossref:
  custom:
    - kind: float
      key: suppfig
      latex-env: suppfig
      reference-prefix: Fig. S
      space-before-numbering: false
      latex-list-of-description: Supplementary Figure

editor: source
execute:
  echo: false
warning: false
---

```{r}
#| label: setup
#| echo: false
#| eval: true
#| tidy: true

# Core packages
library("tidyverse")
library("data.table")
library("ggpubr")
library("ggfortify")
library("rstatix")
library("broom")
library("readxl")

# Linear Models
library("car")
library("nlme")
library("MASS")
library("MuMIn")
library("sjPlot")
library("class")
library("caret")

# Community Diversity
library("vegan")
library("hillR")
library("geosphere")
library("indicspecies")

# ggplot2 extensions
library("ggtext")
library("ggpmisc")
library("MetBrewer")

# Tables
library("huxtable")
library("flextable")
library("broom.mixed")
library("officer")
library("knitr")

# Maps
library("ggmap")
library("ggspatial")

# Phyloseq friendly
#Some of these are redundant and have conflicts, but have functions that I like better than others or compliment each other.
library("phyloseq")
library("microeco") # New package for microbial analyses
library("file2meco") # File to microeco (phyloseq friendly)
library("metagMisc") # Miscellaneous functions
library("microbiome") # Expands phyloseq 
library("mirlyn")

# Parallel processing
library("parallelly")

# Misc
library("conflicted")
conflict_prefer("select", "dplyr")
conflict_prefer("desc", "dplyr")
conflict_prefer("filter", "dplyr")
conflict_prefer("alpha", "ggplot2")

# File path

path <- "../mim2_bioinformatics" # Change accordingly

# Main data sets
# The data cleaning and wrangling of `leaf_traits` and `plant_traits` data sets are in the `mim2_statistics.qmd` notebook, chunk `cleaning_chaping_data`. The data sets are saved as RDS files and CSV files in the `statistics` folder. The data cleaning and wrangling for phyloseq objects are in the `mim2_bioinformatics.qmd` notebook. The phyloseq objects are saved as RDS files and CSV files in the `taxonomy` folder.

# The FEF community data is only available for 160/320 samples. For the downstream analyses, we eliminated *M. bicolor* or species "B" in the data set. The removal of this is due to a small sample size for leaf trait measurements and even smaller for FEF community data (*n* = 3). The removal is not explicit in `plant_traits` it is only explicit in the subset `leaf_traits_noB` and latter on in the `asv_matrix` that gives way for community analyses in [Community Diversity] in the mim2_statistics.qmd. We include a custom function for p-value formatting.

# We will load the data sets and phyloseq objects here.


#### Load leaf traits master files ####
traits <- read.csv(file.path(path, "field_data/mim2_leaf_traits_field_survey.csv")) #Raw file just in case

leaf_traits <- readRDS(file.path(path, "clean_data/statistics/leaf_traits.rds"))

plant_traits <- readRDS(file.path(path, "clean_data/statistics/plant_traits.rds"))

plant_traits_MB <- readRDS(file.path(path, "clean_data/statistics/plant_traits_MB.rds")) #Substitues ps_melt df

#Removing species B
leaf_traits_noB <- leaf_traits |>
  filter(!Species == "M. bicolor")

#### Phyloseq objects ####
#Phyloseq joins various objects that we have already prepare: taxonomic table, ASV table and our sample data.
#These are the data frames resulting from the phyloseq section in mim2_bioinformatics.qmd notebook.

#Load Phyloseq object cleaned of singletons (231 ASVs)
ps_clean_3 <- readRDS(file.path(path, "clean_data/taxonomy/02-TAXA_8450_phyloseq_nonsingletons_noB.rds"))


pseq_rrfb <- readRDS(file.path(path, "clean_data/statistics/pseq_rrfb.rds")) #Mirlyn object from Method B

rarefied_phyloseq <- readRDS(file.path(path, "clean_data/statistics/rarefied_phyloseq.rds"))

# Load ps_melt
#ps_melt <- readRDS(file.path(path, "clean_data/statistics/ps_melt.rds"))

#### Matrices ####
# Initial matrices, before rarefaction
# ASV abundance matrix
# Select randomized samples for beta diversity rarefaction
# A total of 84 samples result from method B
final_names_methodB <- sample_names(pseq_rrfb[[1]]) # Sample names

# # Matrix
# Simple Matrix of 84 samples from Method B pipeline
asv_matrix <- otu_table(ps_clean_3) |> # ASV matrix
  as.data.frame() |>
  select(contains(final_names_methodB)) |> # Samples from method B randomization
  as.matrix()

# Load
rrfy_hell_matrix <- readRDS("clean_data/statistics/rrfy_hell_matrix.rds")

# Mean Hellinger transformation of rarefied data sets
# For Mantel tests
# set.seed(433)
# asv_avgdist <- avgdist(t(asv_matrix),
#                        sample = 750,
#                        distfun = vegdist,
#                        dmethod = "hellinger", #Transformation, later Bray_Curtis
#                        meanfun = mean,
#                        iterations = 50) |>
#   as.matrix()
# 
# # Save
# saveRDS(asv_avgdist, file = "clean_data/statistics/asv_avgdist.rds")

# Load  
asv_avgdist <- readRDS("clean_data/statistics/asv_avgdist.rds") 

# # Adding column and rownames
# rownames(asv_avgdist) <- final_names_methodB
# colnames(asv_avgdist) <- final_names_methodB
# 
# # Geographical matrix
# geo_distm <- distm(plant_traits |>
#                       filter(Unique_ID %in% final_names_methodB) |>
#                       column_to_rownames(var = "Unique_ID") |>
#                       select(Longitude, Latitude), fun = distVincentyEllipsoid)
#                       
# # Adding column and rownames
#rownames(geo_distm) <- final_names_methodB 
#colnames(geo_distm) <- final_names_methodB  

# Save
#saveRDS(geo_distm, file = "clean_data/statistics/geo_distm.rds")
# Load
geo_distm <- readRDS(file.path(path, "clean_data/statistics/geo_distm.rds"))

# Traits for dbRDA
dbrda_hell_matrix <- rrfy_hell_matrix |>
  column_to_rownames(var = "X") |>
  select(-c(1, 233:256))

dbrda_traits <- rrfy_hell_matrix |>
  select(
    X,
    Unique_ID,
    Site,
    Species,
    Elevation_m,
    Elevation_cat,
    Sample_Date,
    logACI,
    logLT,
    logLPS,
    logLMA,
    logLBI
  ) |>
  unite(Species_Elevation,
        Species,
        Elevation_cat,
        sep = "_",
        remove = FALSE) # To model interactions visually.

# Dataset with missing values imputed
dbrda_traits_imp <- dbrda_traits |>
  mutate(across(c(, 8:13), ~ ifelse(is.na(.), median(., na.rm = TRUE), .)))

#### Summary stats about the data ####

# Leaf level traits
#sum(!is.na(leaf_traits_noB$ACI)) # No missing values
#sum(!is.na(leaf_traits_noB$LT)) # No missing values
#sum(!is.na(leaf_traits_noB$LPS)) # No missing values
#sum(!is.na(leaf_traits_noB$LMA)) # No missing values
#sum(!is.na(leaf_traits_noB$LBI)) # No missing values

# Plant level FEF community
# Excludes species "B"
#sum(!is.na(plant_traits$Observed[plant_traits$Species != "B"])) # No missing values
#sum(!is.na(plant_traits$hill_1[plant_traits$Species != "B"])) # No missing values
#sum(!is.na(plant_traits$hill_2[plant_traits$Species != "B"])) # No missing values

# Functions
source("src/mim2_function.R")

#Comparisons
species_comparisons <- list(c("M. nasutus", "M. laciniatus"), c("M. nasutus", "M. guttatus"), c("M. laciniatus", "M. guttatus"))

elevation_comparisons <- list(c("LOW", "MID"), c("LOW", "HIGH"), c("MID", "HIGH"))
```

## Keywords

# Abstract

The *Mimulus guttatus* species complex is a phenotypically and ecologically diverse group of flowering plants with leaf functional trait variation that provides key advantages for adapting to harsh environments. The potential for microbial symbionts to contribute to their ecological and evolutionary success, however, has been overlooked. We explored how leaf functional traits and foliar endophytic fungi (FEF) abundance, richness, and community composition in three sympatric Monkeyflowers change across elevation in the Sierra Nevada, CA, USA. We asked: Q1) Are there differences in leaf functional traits and FEF communities *among* sympatric *Mimulus* spp. populations along at similar elevations? Q2) How do traits and FEF communities change across an elevation gradient within species? Q3) Are FEF richness, diversity and community composition correlated with leaf functional traits and/or elevation? Q4) How does FEF community composition differ with geographic distance within each species? We collected *M. guttatus*, *M. nasutus*, and *M.laciniatus* individuals from natural populations across the Sierra Nevada, measured leaf functional traits and used ITS sequencing to describe the leaf endophyte community. We found  significant associations of FEF community composition with 
<!-- leaf lobing,  -->
host species, and elevation, suggesting that these factors may influence fungal community composition. Furthermore, FEF community dissimilarity was correlated with geographical distance indicating isolation by distance and limited dispersal of fungal endophytes. We detected the prevalence of *Vishniacozyma victoriae*, an endophyte found most commonly in Antarctica, across all Sierran *Mimulus* populations. The presence of *V. victoriae* could play a role in the adaptation of *Mimulus* to cold, high elevation environments.  Our findings offer novel insights into the intricate interactions between fungal endophyte communities, plant traits, and elevation.


# Introduction

Microbial dispersal is challenging to trace. The general assumption is that most, if not all, microbes are ubiquitous. Nevertheless, the dispersal of microbial communities is not random as it is the result of neutral and selective processes influenced by environmental conditions, geographic distance, and host species [@rebolledagomez2019; CITE MORE]. The parallel evolution and symbiotic interactions with plants microbes, such as fungi, aided plants’ colonization of land across multiple varying climates [reviewed in @peay2016; @remy1994; @wang2006; @field2015] and bacteria [@soltis1995; @adams2002; @delaux2015]. However, little is known about how plant microbiomes influence or are influenced by local adaptation of their host species. Geographic location and environmental conditions are known to influence community composition of symbiotic fungi, with varying effects across different fungal guilds [@bowman2021; CITE]. Previous studies have suggested that foliar endophytic fungi (FEF), which live inside plant leaves, can alter plant traits under stressful conditions such as drought [@song2016] and therefore may contribute to local adaptation. Thus, the resulting FEF community is influenced by the interplay of environmental conditions, geographic distance and host species. Understanding the factors that shape FEF communities is essential for understanding the ecological and evolutionary dynamics of plant-microbe interactions.

Host traits and host identity can influence FEF communities in a variety of ways. Leaf functional traits, such as leaf thickness, leaf mass per area, and leaf toughness are known to influence plant's response to abiotic factors like light, temperature, and hydraulic constraints [@nicotra2011; @oguchi2018]. These traits serve as host imposed filters that play role structuring FEF communities [@saunders2010; @tellezTraits2022]. Jointly, host identity has been shown to structure the foliar microbiomes of plants. For example, a study on the classical plant model, *Arabidopsis thaliana*, revealed that bacterial and fungal communities in leaves are shaped by host genotype, at least for the most common operational taxonomic units (OTUs) [@horton2014]. Recently, Karasov et al. [-@karasov2024], found robust trends in support of this hypothesis, showing that the bacterial leaf microbiome in diverse European populations of *A. thaliana* is influenced by host genotype. At the local scale, *A. thaliana* can impose strong selection on the composition of it leaf microbiome, and at a continental scale these effects are ubiquitous, but geography and abiotic factors contribute strongly to the patterns observed [@karasov2024]. The results and patterns observed in this model system encourage the idea to test these hypotheses in other plant species to determine their broader applicability. The *Mimulus guttatus* species complex lends itself be an ideal model system to explore questions about FEF communities and their relationship with host traits and environmental conditions due its closely related species that are interfertile, considerable morphological variation, and wide range of habitats [@vickery1978; @wu2008].

Multiple studies have examined various tissue types for microbial symbionts in *Mimulus* species, yet little is known about their foliar symbionts. Aboveground, Beslile et al. [-@belisle2012] reported on the distribution of diverse fungal communities in the flower nectar of *Mimulus aurantiacus*. The authors considered flowers as islands in a metapopulation system and found that the frequency of micro-fungi (i.e., yeasts) was significantly correlated to the location of the plant and marginally correlated with the density of the flowers in the plant [@belisle2012]. Another study by Rebolleda et al [-@rebolledagomez2019], examined epiphytic bacterial communities of different flower parts in *Mimulus guttatus* and that pollinator selection drives strong microenvironmental selection of bacterial communities, overwhelming the effects of geographic distance. Belowground, arbuscular mycorrhizal fungi (AMF) are known to associate with *M. guttatus*. In an experimental study, McIntosh et al. [-@mcintosh2024], observed that coastal perennial and montane annual ecotypes of *M. guttatus* harbored distinct fungal communities, with the coastal perennial ecotype being more colonized by AMF than the montane annual ecotype. The authors suggested that population divergence in life history and related traits shapes plant responsiveness to soil microbes and the specificity of their associations with root endophytic fungi, such as AMF [@mcintosh2024]. 

Foliar symbionts can increase plant fitness under stressful conditions [@giauque2019; @aimone2023; @shankarnaik2019; @karasov2024]. Considering the wide range of habitats the *M. guttatus* species complex occupies, ranging from granite outcrops to lush alpine meadows, as well as it is crucial to examine the composition of FEF communities across *Mimulus* species to determine if they align with broader ecological patterns. Altitudinal gradients provide excellent natural laboratories to explore environmental effects on plant-microbe interactions. For example, a study by Kazenel et al. [-@kazenel2019], found that the composition of root and leaf endophytes in cool-season, perennial grasses, varied differently with altitude and warming. Others have found that elevation and site are the strongest predictors of fungal taxa [@cordier2012]. To our knowledge, no previous studies have considered symbionts in leaf tissue of *Mimulus* or how host species and elevational gradients influence symbiont community composition. 

In this study we use amplicon sequencing to explore how FEF communities vary among sympatric populations of *M. guttatus*, *M. laciniatus* and *M. nasutus* collected along elevation gradients throughout the central Sierra Nevada, CA ([Fig. @fig-figure1]). We test the following questions: Q1) Are there differences in leaf functional traits and FEF communities among sympatric *Mimulus* species populations at similar elevations? Q2) How do traits and FEF communities change across elevation within species? Q3) Are FEF richness, diversity and community composition correlated with leaf functional traits and/or elevation? Q4) How does FEF community composition differ with geographic distance within each species? We predict that abundance, diversity,and richness would decline with increased elevation, and community composition of FEF would be more similar among the same sites (alpha diversity) than between sites (beta diversity) regardless of host species. With these questions we gauged whether FEF community responses to altitude were consistent across host species. Our results consider replicated altitudinal gradients and provide insights into the ecological dynamics of plant-microbe interactions in the *Mimulus* species complex.

<!-- In addition to the central role of photosynthesis for explaining leaf shape, emphasizing specific leaf functional traits such as leaf thickness (LT), leaf mass per area (LMA), and leaf toughness can illustrate how mechanical properties impact plants' response to their abiotic and biotic environments [@oguchi2018]. Such characteristics are part of a leaf economic spectrum, with leaves at one extreme having short lifespans, high nitrogen content, low LMA, thin leaf blades, and thin cell walls. At the opposite end are long-lived leaves with low nitrogen content, high LMA, thick leaf blades and thick cell walls [@mason2015; @wright2004]. These traits not only determine the chemical, structural, and longevity attributes of leaves but influence their associated microbial communities [@saunders2010; @tellezTraits2022].  -->

<!-- Leaves are essential structures for plants, playing a crucial role in various physiological processes such as photosynthesis, transpiration, and nutrient uptake [@oguchi2018; @fritz2018; @tsukaya2018]. They display various shapes driven by both genetic and environmental factors [@tsukaya2018; @tsukaya2005; @nicotra2011; @fritz2018]. Since photosynthetic leaves are crucial for plant growth and survival, it has been proposed that natural selection has shaped leaf functional traits in response to environmental conditions and ecological demands [@jump2005; @oguchi2018; @nicotra2010]. The leaf functional traits involved in the resulting leaf shape respond to three main stimuli: light, temperature modulation, and hydraulic constraints. Therefore, proposed explanations for the evolution of different leaf shapes have included: 1) thermoregulation to arid and hot environments, 2) hydraulic constraints and biomechanical constraints, 3) adaptations against plant enemies, and 4) adaptations to optimize light interception [@nicotra2011; @oguchi2018].  -->
<!-- These are not mutually exclusive, and the expected patterns are not always clear, for example, the optimal leaf width, the distance from mid-vein to the leaf margin, is driven evolutionarily by the environmental conditions experienced by species [@tsukaya2018].   -->

 <!-- For example, Wang et al. [-@wang2022] provide evidence of the favored genetic pathways for leaf shape differentiation in a comparative study  *Arabidopsis thaliana* which has narrow entire leaves, and close relative *Cardamine hirsute* which has deeply lobed leaves.  -->

# Materials and Methods

## Sample Collection

We collected plant specimens during April-July 2021 and 2022 from populations of *M. guttatus*, *M. laciniatus*, and *M. nasutus* (syn. *Erythranthe guttata*, *Erythranthe laciniata*, and *Erythranthe nasuta*) across Stanislaus National Forest (SNF), Sierra National Forest (SINF) and Yosemite National Park (YNP), CA, USA. <!-- You should create a table of each location and description of every population that you sampled for this paper. That table should be cited here. -->  We haphazardly selected sites close to the main roads based on the presence of a viable population with at least \~ 50 individuals per species. Samples collected from YNP were collected from non-wilderness areas on the side of the road. We determined population viability ensuring that they had individuals flowering or close to flowering stage. We collected between 5 - 20 individuals per species from a total of 25 sites (Table 1). We selected individuals that possessed healthy looking leaves and no visible signs of pathogen damage or senescence. At sites where two species were present, we collected individuals that were at least \~ 25 meters apart. We collected sample specimens by carefully uprooting the plant and placing into individual plastic bags (e.g., Ziploc®) and preserving in an ice chest until return to the field laboratory at the UC Merced Yosemite Field Station, YNP, CA, USA. Plant specimens were processed within 8 hrs of collection.

## Leaf traits measurement

From each plant, we measured leaf traits: leaf thickness (LT), leaf punch strength (LPS), leaf mass per area (LMA), anthocyanin content index (ACI) which are known to be associated with the structure of FEF communities [@tellezTraits2022] as well as leaf lobe index (LBI) [@ferris2015]. We cleaned plants with tap water to remove all soil and debris remnants from the leaves and roots. We removed all healthy leaves (\~ 5 - 10) from the stems and took three measurements per trait from three haphazardly selected leaves from individual plants, except for LBI, only one leaf per plant. We used a transparency film to hold the leaf in place and flatten, after which we took a digital photograph for further analysis. To calculate the LBI, we followed Ferris et al. [-@ferris2015]. Briefly, leaf lobing is calculated as the convex hull area minus the true leaf area divided by convex hull area of a digital photograph of a leaf in ImageJ [v1.52r; @schneider2012]. We measured ACI content with ACM-200plus (Opti-Sciences Inc. Hudson, New Hampshire, U.S.A.) on haphazardly selected locations of the leaf surface (working from the petiole out to the leaf tip) [@tellezTraits2022]. The ACM-200 calculates an ACI value from the ratio of % transmittance at 931 nm/% transmittance at 525 nm [@opti-sciencesinc], effectively accounting for leaf thickness. We measured LT (μm) with a Mitutoyo 7327 Micrometer Gauge (Mitutoyo, Takatsu-ku, Kawasaki, Japan) on haphazardly selected locations of the leaf lamina, taking care to avoid major and secondary veins. We used an Imada DST-11a digital force gauge (Imada Inc., Northbrook, IL, United States) to measure LPS, a measure of leaf toughness, on the lamina of each leaf selected, avoiding minor leaf veins when possible [@tellezTraits2022]. It functions by conducting punch-and-die tests with a sharp-edged cylindrical steel punch (2.0 mm diameter) and a steel die with a sharp-edged aperture of small clearance (0.05 mm). Once LPS was measured, we used a 4 mm diameter punch hole to puncture disks for LMA measurements. We collected one disk per leaf (see Supplementary material for details). The disk punches dried were shipped to Tulane University, New Orleans, LA, USA to dry at 60 ℃ for 48-72 hours before being weighed.

## Molecular Work

### Tissue preservation

Upon completion of the leaf traits measurements, we prepared and preserved samples at the UC Merced Yosemite Field Station. We started by removing the main vein and margins from photosynthetic tissue. The leaf lamina was haphazardly cut with a sterile blade into 2 mm wide strip in parallel to the main vein [@arnold2003; @tellezTraits2022; @higgins2014]. Leaf strips were then sterilized with sequential washes in 95% EtOH (10 s), 0.5% sodium hypochlorite (NaOCl) (60 s), and 70% EtOH (60 s) and air dried under sterile conditions. Due to the small size of monkeyflower plants, the maximum amount of leaf lamina was preserved in sterile 15 mL tubes with \~ 10 mL CTAB solution (1 M Tris–HCl pH 8, 5 M NaCl, 0.5 M EDTA, and 20 g CTAB). The leaf tissue in CTAB solution was used for amplicon sequencing (described in detail below). All leaf tissue handling was performed in a sterile environment with an alcohol burner lamp inside a portable biosafety cabinet. All surfaces were previously sterilized sequentially with 0.5% NaOCl, 95% EtOH, and 70% EtOH. We surface sterilized surfaces and instruments in between sample handling to prevent cross contamination.

### Amplicon sequencing

We stored leaf tissue in CTAB solution for 2 months at room temperature before extracting DNA at Tulane University. To prepare for sample DNA extraction procedure, we decontaminated all instruments, materials, and surfaces in biosafety cabinet with 0.5 % NaOCl, 70 % EtOH, and 95% EtOH, and subsequently treated with UV light for 30 minutes. We subsampled 0.2 - 0.3 g of leaf tissue from each sample and placed into a sterile 2 mL tubes containing an assortment of beads: 3.2 mm stainless steel beads (Next Advance, Cat# SSB32), 100 µL stainless steel bead blend, 0.9-2.0mm (NextAdvance, Cat# SSB14B) and 2-3 of the autoclaved 2 mm zirconium oxide beads (Next Advance, Cat# ZRoB20). The 2 mL tubes with beads were previously prepared. We then proceeded to lyophilize samples for 72 hours to fully remove CTAB content from tissue. After, we submerged the sample tubes in liquid nitrogen for 30 s and homogenized samples at 30 Hz for 3 minutes in a TissueLyser LT (QIAGEN, Valencia, CA, USA). We stored samples in 20 ℃ until DNA extraction procedure.

We used a DNA extraction protocol for high-molecular weight DNA extraction adapted from Russo et al., [-@russo2022]. Briefly, it is a CTAB:chloroform:isoamyl extraction combined with a solid-phase reversible immobilization (SPRI) bead step [@russo2022; @rohland2012; @liu2023]. Protocol modifications allowed us to optimize extractions for fungal DNA from preserved leaf tissue [see details in @aponterolon2023]. After all genomic DNA was extracted, we quantified the DNA using Quant-iT® dsDNA HS Assay kit with Qubit Flourometer (Thermo Fisher Scientific, Waltham, MA, USA., Cat# Q33120) and followed a two-step amplification approach described by Sarmiento et al. [-@sarmiento2017] and U'Ren & Arnold [-@uren2017]. We used standard primers ITS1F [@gardes1993] and ITS2 [@white1990_inbook] modified with the Illumina TruSeq adaptor (see Supplementary Material). Every sample was amplified in three parallel reactions at the annealing temperatures 52 ℃, 54 ℃, and 56 ℃ to amplify a wide range of fungal taxa and reduce amplification bias for short ITS sequences [@uren2017; @lumibao2018]. Each PCR (PCR1) reaction contained 2 µL of sample DNA template. We visualized PCR1 reactions with SYBR™ Safe DNA Gel Stain (Thermo Fisher Scientific, Waltham, MA, USA., Cat# S33102) on 2% agarose gel [@oita2021]. We combined 5 µL of amplicon product from parallel reactions into a single tube per sample and purified using Sera-Mag™ SpeedBead Carboxylate-Modified Magnetic Particles (Hydrophobic) (Thermo Fisher Scientific, Waltham, MA, USA., Cat#09-981-123) prepared as per [-@rohland2012 and -@liu2023] and used a ratio of 1.2x:1 with 80% EtOH following manufacturer's instructions. We used 3 µL of PCR1 product from samples, DNA extraction controls, and PCR1 negative controls for a second PCR (PCR2) with barcoded adapters (IDT, Coralville, Iowa, USA). Each PCR2 reaction (total 30 µL) contained 1X Phusion Flash High Fidelity PCR Master Mix (Thermo Fisher Scientific, Waltham, MA, USA., Cat# F548L), 0.075 µM of barcoded primers (forward and reverse pooled at an initial concentration of 2 µM) and 0.20mg/mL of BSA (Thermo Fisher Scientific, Waltham, MA, USA., Cat# B14) following [-@sarmiento2017 and -@uren2017]. Before final pooling for sequencing, we purified and concentrated amplicons using SPRI beads to a total volume of 20 µL. We quantified PCR2 product with Quant-iT™ PicoGreen™ dsDNA Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA., Cat# P7589) with the BioTek Synergy LX plate reader (Agilent, Santa Clara, CA) and combined equimolar amounts of libraries, including DNA extraction controls, PCR1, and PCR2 negative controls into a 10nM library pool. We did not detect any contamination visually or fluorometrically. Libraries were sequenced on the Illumina MiSeq platform with Reagent Kit v3 (2 \u00D7 300 bp) at Duke Genome Sequencing and Analysis Core Facility (Durham, NC, USA). Throughout all these steps, we used a separate set of sterile pipettes, tips, and equipment to reduce contamination in a designated PCR area to restrict contact with pre-PCR materials [@oita2021].

### Bioinformatic analyses

We assessed the quality of the reads using FastQc v0.12.1 [ v0.12.1; @andrews2010] and MultiQC [@ewels2016] tools. A total of 60,696,808 total ITS1 reads yielded from 343 (including 27 controls) libraries sequenced in two separate sequencing runs. The first sequencing event yielded 32,117,684 and the second 28,579,124 ITS1 reads. We tailored the open-source DADA2 bioinformatic pipeline for our data set [@callahan2016]. Based on our initial quality assessment, both forward and reverse reads were of low quality, with base calls deteriorating after 100 bp. We filtered our reads for ambiguous calls before removing the adapters by using `filterAndTrim` function and argument `maxN = 0` from the `dada2` package [v1.28.0; @callahan2016]. We removed forward and reverse primer adapters (and their reverse compliments) and eliminated reads shorter than 20 bp using the `cutadapt` tool [v4.6, @martin2011]. After removing ambiguous calls and forward and reverse primers, we assessed the quality of the reads again and saw slight improvement. We proceeded to apply stringent filter and truncation parameters to ensure quality of reads when assigning taxonomy. We filtered and truncated reads based on maximum expected errors (maxEE) rather than read length as it provides a reliable quality filtering [@edgar2015]. For this we set the arguments `trunQ = 2`, `maxEE = c(2,2)` for forward and reverse reads, and minimum read length of 50 bp with `minLen = 50` in the used the `filterAndTrim` function [@callahan2016]. These parameters eliminated 151 samples from our data set, all from our second sequencing event. After this filter, we dereplicated reads with the `derepFastq` function and merged pairs using `mergePairs` functions with an overlap of 20 bp, minimum. We then inferred composition of the samples with `dada` function, which applies the DADA algorithm [@callahan2016; @rosen2012]. We removed chimeras via the "consensus"method with the `removeBimeraDenovo` function and ultimately we used the `assignTaxonomy` function to assign taxonomy the amplicon sequence variants (ASV) referenced against the UNITE database [@abarenkov2023a]. After taxonomy assignment we used the `phyloseq` package [@mcmurdie2013] to create a phyloseq object for downstream analyses.

We used the `decontam` package [v1.20.0; @davis2018] to statistically determine which ASVs are likely contaminants based on their frequency in our samples and remove them using `prune_taxa` function from the `phyloseq` package [v1.44.0; @mcmurdie2013]. After which, we calculated the average read count found in DNA and PCR extraction controls, considered to be laboratory contaminants, and subtracted that from the samples' read counts. We then used custom scripts to remove any ASV that represented less than 0.1% of the abundance per sample on the assumption that it originates from contamination throughout handling of samples in the DNA and PCR processes. We removed singleton ASVs with the `prune_taxa` function [@mcmurdie2013]. We identified core taxa members at a 1% detection threshold and 50% prevalence in samples using the `core_members` function from the `microbiome` package [v.1.22.0; @lahti2012]. All post-quality bioinformatic steps were performed in *R* [v.4.4.1; @rcoreteam2024].

## Statistical Analyses

### Community Diversity

To account for uneven sampling effort and over-representation of sequences, we normalized libraries by repeated rarefying, following Cameron et al. [-@cameron2021]. We determined a sequence depth of 750 by calculating Good's coverage and qualitative evaluation of libraries to determine a balanced coverage and breadth of samples [@schloss2024]. This approach allowed for a proportionate representation of observed sequences from host species and a robust characterization of random variation inherent in rarefaction of small libraries [@cameron2021; @schloss2024]. First, we randomly selected 136 samples out of the 157 that resulted from the bioinformatic pipeline (refer to the Results section below). The sample pool was reduced to 136 to match the theoretical reduction of one sample per site (*n* = 21). We generated 50 rarefied abundance matrices without replacement by using the `mirl` function from the `mirlyn` package [@cameron2021]. After which, we calculated alpha diversity per sample as Hill orders, the observed richness (𝑞 = 0), the exponent of Shannon's entropy (𝑞 = 1), and the Inverse Simpson's Diversity (𝑞 = 2) by applying a modified version of the function `alphadivDF` [@cameron2021], which wraps the common diversity indices in `vegan` [@oksanen2022] (see custom script in Supplementary material). Finally, for beta diversity analyses, we performed a Hellinger transformation on the rarefied abundance matrices and calculated a Bray-Curtis dissimilarity matrix for each.

 We calculated simple linear regressions understand how different measures of alpha diversity changed in response to elevation. To answer how FEF communities differ *among* host species (Q1) and *between* sites (Q2) and facilitate our understanding of the effects of elevation on FEF communities, we categorized elevation as "low" (\< 1219 m.a.s.l), "mid" (1220 - 1828 m.a.s.l.) and "high" (\> 1829 m.a.s.l.). Additionally, we applied a distance-based Redundancy Analysis (dbRDA) on the Bray-Curtis dissimilarities to statistically compare the FEF community similarities within each host species per site (Q3). Its visualizations effectively portray underlying patterns of compositional differences [@anderson2017; @legendre1999; @mcardle2001], akin to permutational analysis of variance. We calculated Spearman's rho for all leaf traits (see below and Fig. S1) and informed our selection of leaf traits with results from the PCA (see below and Fig. 1), and selected those with the lowest correlation coefficient per host species: logLBI. Our initial dbRDA model consisted of terms logLBI, sampling date, site, and elevation (m). The leaf functional trait data, as well as elevation, was not randomized or subsampled to match rarefied dataset, the same values applied to all 50 rarefied matrices. For our initial model, we determined the variance inflation factor (VIF) of each term with function `vif.cca` and eliminated redundant ones: site and sampling date. We performed manual model selection by evaluating the marginal significance of constraining variables after 999 permutations with function `anova.cca` and argument `by` set to "terms" to assess significance for each term [@legendre2011; @legendre2012; @oksanen2022]. We corroborated homogeneous dispersion of variances with a permutational analysis of multivariate dispersion (PERMDISP) using the `betadisper` with parameter `type = "median"`, and `permutest` functions from `vegan`, the latter with 999 permutations [@oksanen2022]. We used a post-hoc Tukey's test to compare the differences in the dispersion of the FEF communities among sites and elevation categories. We used the `simper` function from `vegan` to discriminate which species contribute the most to compositional differences between groups [@oksanen2022].
<!-- To determine if there is a core FEF community associated with host species or elevation categories, we conducted an indicators species analyses (ISA) using the `multipatt` function from the `indicspecies` package [@decaceres2009].  -->
We used generalized linear mixed models (GLMM) to determine the effect of host species, elevation, and leaf functional traits on the mean Bray-Curtis dissimilarity (β-diversity) of FEF communities. Leaf functional traits with missing values were imputed with the median value for the trait. We used the `lme` function from the `nlme` package [v3.1-165; @pinheiro2024; @pinheiro2000] to fit the GLMM models. We established leaf functional traits, elevation and host species as fixed terms and site as a random effect. We modeled variance structure for site with the `varIdent` argument. We manually compared and selected models based on Akaike Information Criterion (AIC) with a penalty of 2 degrees of freedom (ΔAIC) with the `model.sel` function from the `MuMIn` package [@barton2023]. We selected the best-fit model based on the lowest value obtained. All models are based on restricted maximum likelihood estimation (REML).
 
To assess correlations between pairwise FEF community dissimilarity and the geographical distance matrix per host species (Q4), we computed Mantel tests with Spearman's rho and 999 permutations using the `mantel` function [@oksanen2022]. For this test, we opted for a less taxing computational approach and rarefied sequences with the same parameters as before and calculated Hellinger transformations with the `avgdist` function [@oksanen2022]. We then calculated the Bray-Curtis dissimilarity with `vegdist` [@oksanen2022]. For the geographical distances between sites, we used `distm` function with the Vicenty (ellipsoid) method from the `geosphere` package [@hijmans2022].

### Leaf traits

We checked for normality and homoscedasticity of the leaf traits measured. We used Shapiro-Wilk and Fligner-Killen tests from the `stats` package [@rcoreteam2024] to check for normality and homoscedasticity, respectively. We established that the leaf functional trait data was not normally distributed and not homoscedastic. We then used non-parametric tests, the Wilcoxon Rank Sum test, to compare leaf functional trait means among species and sites to answer the first portion of Q1 and Q2. We used the `compare_means` and `stat_compare_means`functions from the `ggpubr` package [@kassambara2023] to perform these tests and properly visualize them. We adjusted *p* values to account for false discovery rates in multiple comparisons by using “BH” method [@benjamini1995]. We performed Principal Component Analysis (PCA) to understand patterns and relationships among leaf traits of host species. We used the `prcomp` function from the `vegan` package [@oksanen2022] compute the PCA analysis with variables ACI, LT, LPS, LMA, and LBI, all log-transformed. We used only complete raw leaf functional traits measurements to compute the PCA analysis (*n* = 504), samples with missing values were omitted.  All statistical analyses were performed in *R* programming language [v.4.4.1; @rcoreteam2024]. We present the log-transformed leaf functional trait data at the leaf level: ACI (*n* = `r sum(!is.na(leaf_traits_noB$logACI))`), LT (*n* = `r sum(!is.na(leaf_traits_noB$logLT))`), LPS (*n* = `r sum(!is.na(leaf_traits_noB$logLPS))`), LMA (*n* = `r sum(!is.na(leaf_traits_noB$logLMA))`), LBI (*n* = `r sum(!is.na(leaf_traits_noB$logLBI))`). The FEF community data is presented at the plant/sample level (*n* = `r sum(!is.na(plant_traits$Observed))`).


# Results

Our PCA analysis showed how leaf functional traits were related (Q1). We plotted leaf functional traits according to species groups on the PCA axes to show how the variance in the complete data set is explained by PC1 (42.57%) and PC2 (23.27%) ([Fig. @fig-figure2]). The PCA analyses showed correlations between ACI, LT, LPS, and LMA as loadings that tracked along PC1 towards positive values ([Fig. @fig-figure2]). We observed that LBI loading was orthogonal along PC2 to the other traits, indicative of low correlation. We noted distinct groupings by species along PC2 such that *M. laciniatus*, the most lobed species, was distinct from *M.guttatus* and *M. nasutus*. The latter two overlapped along PC1 and PC2 ([Fig. @fig-figure2]).

We found that leaf functional traits differed among and within *Mimulus* species across elevation (Fig. S5 - S9). For LMA, when we observed all host species, there was a statistically significant positive correlation between LMA and elevation (*R^2^~adj~* = .041, *p* \< .001, Fig. S5A). When we observed this relationship with categorical elevation, we found that species did not differ in LMA values at low elevations (Fig. S5B). At mid elevations, *M. laciniatus* and *M. guttatus* had significant differences in LMA (*p* \< .001), while at high elevation this difference dissipated. At high elevation we saw that *M. nasutus* differed significantly from both *M. laciniatus* and *M. guttatus* (*p* \< .0001, Fig. S5B). For ACI, when we observed all species, we found no correlation with elevation (*R^2^~adj~* \< -.000, *p* \< .0001, Fig. S6A). But when we observed the elevation categories, we saw that *M. laciniatus* had statistically significant differences (*p* \< .0001) from *M. nasutus* and *M. guttatus*, while the latter two did not differ (Fig. S6B). At mid elevation we saw that *M. laciniatus* and *M. guttatus* had statistically significant differences (*p* \< .0001), while at high elevations *M. laciniatus* had significantly lower levels (*p* \< .0001) of ACI than *M. nasutus*, and *M. guttatus* (*p* \< .01, Fig. S6B). We also saw statistically significant differences between *M. nasutus* and *M. guttatus* at high elevations (*p* \< .0001) (Fig. S6B). The LBI trait is a measure of leaf "lobeness" and it is confounded with species, since not all host species exhibit the trait plasticity with elevation change. Regardless, our comparisons show that LBI is significantly correlated with elevation (Fig. S7A). Our categorical comparison shows that *M. laciniatus* had statistically significant greater values of LBI compared to *M. nasutus* and *M. guttatus* at all elevation categories (Fig. S7B). At low elevations, *M. nasutus* and *M. guttatus* did not differ, but they did have significant differences at high elevations (p \< .0001, Fig. S7B). Our measure of leaf toughness, LPS, was significantly positively correlated with elevation (*R^2^~adj~* = .004, *p* = .032, Fig. S8A). At low elevations we saw a statistically significant difference in LPS between *M. nasutus* and *M. guttatus* (*p* \< .01), while at mid elevations *M. laciniatus* and *M. guttatus* showed a similar pattern (*p* \< .0001, Fig. S8B). At high elevations we only saw significant differences between *M. laciniatus* and *M. nasutus* (*p* \< .01, Fig. S8B). Finally, we saw a positive correlation between LT and elevation (*R^2^~adj~* = .013, *p* \< .0001), when we observed all species together (Fig. S9A). We only saw statistically significant differences between *M. laciniatus* and *M. nasutus* at low elevations (*p* \< .05), while at mid we see significant differences between *M. laciniatus* and *M. guttatus* (*p* \< .0001) and at high elevations as well (*p* \< .0001, Fig. S9B). At high elevations, *M. guttatus* and *M. nasutus* also showed significant differences (*p* \< .0001, Fig. S9B).

We obtained 5,082,229 reads representing 726 ASVs from 174 samples after processing samples through the DADA2 pipeline. The raw reads obtained were composed of 26.81% Ascomycota, 71.53% Basidiomycota, \<0.05% Chytridiomycota, \<0.5% Mortierellomycota, \<0.03% Olpidiomycota, 0.01% Rozellomycota, \<0.001% Aphelidiomycota, and 1.19% unidentified. After decontamination, and removal of singletons we retained 4,856,220 reads representing 231 ASVs from 157 samples composed of 26% Ascomycota, 73% Basidiomycota, 0.01% Chytridiomycota, \<0.1% Mortierellomycota, 0.03% Olpidiomycota, \<0.002% Rozellomycota and \<1.0% unknown reads ([Fig. @fig-figure3] and Fig. S2). After rarefaction of sequences, we were left with 84 samples in which we found the most prevalent taxa: *Vishniacozyma victoriae* (Basidiomycota, ASV_1), *Cladosporium herbanum* (Ascomycota, ASV_2) and *Cladosporium* spp. (Ascomycota, ASV_7), *Dyszogia patagonica* (ASV_3), *Filobasidium chernovii* (Basidiomycota, ASV_5), and *Alternaria tenuissima* (Ascomycota, ASV_8) (Fig. S2).

In general, we observed a decrease in FEF endophyte diversity with increasing elevation (Q1). At the genus level, we observed a negative correlation between elevation and observed richness (𝑞 = 0) (*R^2^~adj~* \< .01, *F*~1,~ ~4198~= 16.69, *p* \< .001), the exponent of Shannon's entropy (𝑞 = 1) (*R^2^~adj~* \< .01, *F*~1,~ ~4198~= 37.91, *p* \< .001), and the Inverse Simpson's Diversity (𝑞 = 2) (*R^2^~adj~* \< .01, *F*~1,~ ~4198~= 23.73, *p* \< .001). At the species level, we observed an increase in alpha diversity for *M. laciniatus* ([Fig. @fig-figure4] A - C) while *M. nasutus*' declined with elevation ([Fig. @fig-figure4] A - C). *M. guttatus*' alpha diversity increased for Hill number 0 (𝑞 = 0) but did not change for the other Hill numbers ([Fig. @fig-figure4] A - C). The alpha diversity measures showed significant differences among all host species, except *M. laciniatus* and *M. nasutus* in Hill number 0 (𝑞 = 0, *p* \> .05, Fig. S4A). We observed a similar pattern in beta diversity among elevation sites (Fig. S4D -S4F). We saw no differences in beta diversity between low and mid elevation sites for Hill order 2 (𝑞 = 2, *p* \> .05, Fig. S4F).

The FEF community composition varies by species and elevation category (Q2). The best-fit model revealed that 22% of the overall variance in FEF communities was accounted for by log-transformed LBI, and the interaction of host species and elevation, the constraining variables. Each of these variables was significantly correlated  (*p* \< .001) with FEF communities in their host species ([Fig. @fig-figure5] and Table S3). We saw that the first axis (dbRDA1) explained 51.7% and the second axis (dbRDA2) explained 15.9% of the constrained variance ([Fig. @fig-figure5]). We observed high overlap in the groupings of FEF communities by host species at LOW elevations and greater clustering at HIGH elevation sites ([Fig. @fig-figure5]). We found that FEF communities were not homogeneous across all species (*F*~2,~ ~4197~= 320.3, *p* < .001). The post-hoc Tukey’s test revealed that all species comparisons were significantly different at ⍺ = 0.05. We detected significant differences in the dispersion of FEF communities by elevation category (*F*~2,~ ~4197~= 228.1, *p* \< .001) (Tukey’s test ⍺ = 0.05). The interaction between log LBI and elevation category also showed differences in group dispersion (*F*~2,~ ~4192~= 268.1, *p* < .001). The only comparisons that were not different were; *M. guttatus* across all elevations, *M laciniatus* at LOW and MID elevations compared to *M. guttatus* at HIGH elevations, and *M. nasutus* at MID and HIGH elevations compared to *M. laciniatus* at MID elevations (Tukey’s test ⍺ = 0.05). Due to the significant differences between host species and elevation categories in PERMDISP tests we cannot rule out that the observed differences FEF communities are due to dispersion. 

The GLMMs best-fit model, Model 4 (Table 1), included only leaf functional traits as significant predictors of mean β-diversity. The ACI levels were a significant predictor with the greatest effect size (estimate = -0.078, *t*(4200) = -32.48, *p* \< .001). LT was a significant predictor (estimate = -0.051, *t*(4200) = -9.51, *p* \< .001) and LMA as well (estimate = 0.013, *t*(4200) = 3.74, *p* \< .001). Elevation was not a significant predictor of mean β-diversity across all models examined (Table 1). In accordance, simple linear regressions showed that mean β-diversity in *M. laciniatus* did not change significantly with elevation (*R^2^~adj~* \< .01, *F*~1,~ ~1598~= 1.84, *p* = .176, [Fig. @fig-figure6]). At the same time, *M. nasutus* (*R^2^~adj~* = .46, *F*~1,~ ~748~= 642.3, *p* \< .001) and *M. guttatus*  (*R^2^~adj~* \< .01, *F*~1,~ ~1848~= 38.3, *p* \< .001) showed a significant decline in mean β-diversity with elevation.

Finally, we observed a significant correlation between FEF community and geographical distance (Q3). This was only true, however, for *M. laciniatus* (*r* = 0.27, *p* \= .01, [Fig. @fig-figure7]A) and *M. nasutus* (*r* = .29, *p* \< .01, [Fig. @fig-figure7]B). For *M. guttatus* we saw no significant correlation between FEF community composition and geographic distance (*r* = -.004, *p* = .48, [Fig. @fig-figure7]B]).


# Discussion

 <!-- Q1) Are there differences in leaf functional traits and FEF communities *among* sympatric *Mimulus* spp. populations along an elevation gradient?  -->
 <!-- Q2) Are FEF richness, diversity and community composition correlated with leaf functional traits and/or elevation gradient?  -->
 <!-- Q3) How does FEF community composition differ in response to geographic distance? We expected the abundance, diversity, richness, and community composition of FEF would be more similar among the same sites (alpha diversity) than between sites (beta diversity) regardless of host species. -->

<!-- For slender foliage, this could be an adaptation to dry climates and areas near riverbeds prone to flooding. Plants in arid environments face the need to lower leaf temperature and minimize evapo-transpiration, while plants in flood-prone areas must balance photosynthetic capacity with enhanced structural resilience against the shearing forces brought on by floods. [@tsukaya2018; @vansteenis1981].    -->


 <!-- Auxin-related pathways have been well documented as one of the main mechanisms for the expansion of leaf lamina [reviewed in @xiong2019].   -->
Our results show how leaf functional traits differ among and within *Mimulus* species across elevation, point to a significant decline in FEF alpha diversity with increased elevation, and show differences in FEF community composition due to elevation change for all host species. Specifically, LMA, LBI, LPS and LT, were positively correlated with elevation at the genus level, reflecting a conserved plastic response to environmental changes along elevation. When we conducted Principal Component Analysis (PCA) ([Fig. @fig-figure2]), it illustrated correlations between leaf functional traits while demonstrating the distinct groupings of *M. laciniatus*, *M. guttatus*, and *M. nasutus* based on LBI differences compared to the rest of the traits: ACI, LMA, LPS and LT.  *Mimulus laciniatus* is considered the most lobed species ([Fig. @fig-figure2] and Fig. S7) [@ferris2015; @ferris2016; @tataru2023; @sexton2013a], and thrives at higher elevation. For serrated/lobed leaves, it has been observed that physiological adaptations in the leaf margin are a response to the average air temperature, but this is not a simple relationship [@tsukaya2018; @nicotra2011]. Lobed leaves are expected to have increased hydraulic efficiency and have an advantage in arid environments by influencing the leaf boundary layer, an area of still air adjacent the leaf's abaxial surface [@nobel2009; @oguchi2018]. The size of the boundary layer impacts a leaf's gas and heat exchange [@schuepp1993; @nobel2009] affecting its response to changes in mean daily temperature and water availability. The observed increase in LBI with elevation could be a response to the need to increase hydraulic efficiency and reduce heat stress in *M. laciniatus* and *M. nasutus*, while the increase in LMA, LT and LPS with elevation could demonstrate the need to increase leaf structural strength and reduce water loss at higher elevations. An experimental approach with common gardens is needed to fully answer these questions. Experimental manipulation with *Mimulus* is relatively easily due to its short generation time, high fecundity, self-compatibility, and ease of greenhouse propagation [@society2019; @wu2008]

The evolution of leaf shape for at least two of the focal species is controlled by overlapping genetic regions [@ferris2015]. Ferris et al., [-@ferris2015] points to multiple events of leaf shape evolution in the *M. guttatus* species complex and their associated habitats. Recent work has narrowed the genetic architecture at play in incomplete reproductively isolated sympatric populations of *M. laciniatus* and *M. guttatus*, shedding light on the location of quantitative trait loci (QTL) for the first flower node, flowering time, corolla width, corolla length, and leaf shape and suggesting that large-effect loci underlie these traits [@ferris2017]. Future experiments could attempt to align this genetic information with FEF diversity or community composition.

The observed decline in FEF alpha diversity with increased elevation (Fig. 4) echoes macroecological patterns of biodiversity [@sabatini2018; @kraft2011; @villacampa2019; @jimenez-hernandez2020]. Nonetheless, it is filled with nuances for each species. The observed decline in alpha diversity appears to be driven by *M. nasutus* decline and *M. guttatus*' unchanged diversity levels with increased elevation. Our dbRDA results (Fig. 4) indicate that LBI, host species, and elevation, account for 17% of the variance in FEF community composition. Subsequent PERMDISP analyses also support significant differences in FEF community composition between host species at low, mid, and high elevations. In contrast, our best-fit GLMM results suggest that a subset of leaf functional traits are the best predictors of mean β-diversity, and elevation is not. We saw that mean β-diversity declined with increased elevation ([Fig. @fig-figure6]). These conflicting results could be explained by the fact that the GLMM models allowed to model the variance structure of the site origin of host species, which could be confounded with elevation in the dbRDA analyses and not accounted for. Further insight from our Mantel tests indicated a significant positive correlation between FEF community dissimilarity and geographical distance for *M. laciniatus* and *M. nasutus*, providing evidence in support of spatially driven differences in FEF community composition for *Mimulus*. It is interesting to note that the most dissimilar FEF communities due to geographical distance are in the two species that have self-fertilizing mating systems, *M. laciniatus* [@ferris2014] and *M. nasutus* [@brandvain2014]. This could indicate that mating systems and host species' genetic structure could play a role in structuring FEF communities, but further evidence is needed to support this. Our findings overall support the idea that distinct FEF communities are structured by the interplay of host specific leaf functional traits and geographic distance. 

Patterns of diversity and community composition in microbial ecology are often constrained by both biotic and abiotic factors. For example, in an experimental setup, Kivlin et al., [-@kivlin2022] reported that host species (alpine grasses) was a stronger predictor than elevation for of alpha diversity and community composition of leaf endophytes. In contrast, root endophyte communities responded to both host species and elevation [@kivlin2022]. It is possible that our results differ due to the different phenologies and tissue types of herbaceous and gramineous plants. Similarly, Kezenel et al. [-@kazenel2019] found greater change in leaf endophytes due to altitude and warming when compared to root colonizing fungi, but the direction and magnitude of responses varied among host species and fungal functional groups. A major difference in this study is the low ASV count and the use of multiple rarefied data sets compared to Kazenel et al. [-@kazenel2019]. A study by Cordier et al. [-@cordier2012] focused on the fungal phyllosphere in European beech along an elevation gradient, and found that climatic variables, especially temperature, showed the strongest correlations with fungal community dissimilarities. While the phyllosphere of beech varies widely, Cordier et al. [-@cordier2012] also found a strong affinity of fungal taxa to elevation and site, supporting regional spatial structure. An important distinction is that they focused on the outer and inner phyllosphere, hence observing patterns that reflected the outer leaf dynamics. These may have been more susceptible to climatic factors, as opposed to inner leaf dynamics that we explored in our study. Other key differences are the host species phenology and functional leaf traits of European beech compared to *Mimulus* spp.

We saw a consistent prevalence of *Vishniacozyma victoriae*, *Cladosporium herbanum*, and *Cladosporium* in the rarefied data across all sites and samples. *Cladosporium* spp. is a well-documented saprophytic genus that occurs in senescing or dead leaves of herbaceous and woody plants [@samson2004; @schubert2007]. The basidiomycetous yeast *V. victoriae* (formerly *Cryptococcus victoriae*) is also a well-known environmentally abundant fungus capable of causing respiratory issues in humans [@rush2023]. It was first isolated in the Antarctic [@montes1999] and has since been detected worldwide [@demenezes2019]. Despite potential respiratory health detriments, *V. victoriae* has been utilized in agricultural settings for the post-harvest control of fruit diseases [@lutz2020]. It thrives at low temperatures (15 °C), but it is known to tolerate a variety of environmental conditions, and lacks a  polysaccharide capsule, which is thought to contribute to its lack of pathogenicity [@rush2023]. Its applications in wheat agriculture suggest that kernel weight is influenced by *V. victoriae*'s coexistence with other plant acquired endophytic fungi [@vujanovic2021]. Its presence might serve as an indicator of wheat's kernel resistance to pathobiota [@lutz2020]. It is proposed that through the production of various bio-active compounds, it can contribute to plant growth and ecological adaptation to cold environments [@vujanovic2021; @buzzini2018]. According to Vujanovic [-@vujanovic2021] and  Ogaki et al. [-@ogaki2020], no antagonism has been detected between *V. victoriae* with other yeasts, or plant pathogens. We need further quantitative studies to confirm the existence of cold-adapted microbial taxa and their associated hosts [@marian2022]. The presence of *V. victoriae* in our samples might be indicative of its potential role in the local adaptation of *Mimulus* to cold and high elevation environments. 

To optimize microbial studies in *Mimulus*, it is crucial to explore alternative sampling methods such as using fresh tissue or rapidly preserved tissue in liquid nitrogen to improve DNA extraction yields. We urge future investigations to expand their sampling efforts per species and populations while prioritizing prompt DNA extraction to enhance FEF capture. We also acknowledge that experimental manipulations are needed to confirm the causal relationships between leaf functional traits, FEF communities, and elevation. Future studies should consider the role of plant genotype and genetic loci in shaping FEF communities and how these communities might contribute to the host's adaptation to cold environments.

# Conclusions

The *Mimulus guttatus* species complex serves as a robust ecological and evolutionary model system. The identification of FEF communities in *Mimulus* spp. leaf tissue represents a substantial contribution to this field of study, opening up new avenues of inquiry. Our study uncovers potential beneficial FEFs that may contribute to the species complex's adaptation to cold environments. Future research should focus on exploring the interactions of FEF communities and *Mimulus* host genotypes that contribute to the expanded phenotype. We further need to understand how highly prevalent FEF taxa respond to seasonal and temporal changes and contribute to overall plant fitness. An experimental approach taking into consideration populations’ phenotype and genotypes can help disentangle the effect of site and host species on FEF communities. We can understand how plants and their symbionts might respond to climate change – and how symbionts may alleviate plant stress as the planet warms. 

# Author Contributions
B.A.R, K.G.F. and S.A.V. designed research study. B.A.R. conducted field surveys, laboratory work, and collected data. B.A.R. processed the data and performed statistical analyses. B.A.R. wrote the manuscript with input from all authors. All authors gave final approval for publication. 


# Acknowledgements

This work was supported by the National Institute of General Medicinal Sciences of the National Institute of Health (NIH) under award number R35GM138224 to KGF. We thank the Yosemite National Park Service for permit support (YOSE-2021-SCI-0033, YOSE-2022-SCI-0051), and Breeanne Jackson at Yosemite Field Station (doi:10.21973/n3v36c) for logistical support and providing accommodation. We thank the USDA Forest Service for providing access to the Sierra National Forest (MISSING PERMIT NUMBER) and Stanislaus National Forest (Botanical Collection Permit: 031528). We thank Caroline, M. Dong, Fidel J. Machado Perez, Elizabeth MacDougal, Rachael Dennis, and Lissette Montano for field assistance. This research was supported in part using high performance computing (HPC) resources and services provided by Information Technology at Tulane University, New Orleans, LA. Lastly, we would like to acknowledge the original inhabitants of the unceded land on which our research was conducted, the Southern Sierra Miwuk Nation, Bishop Paiute Tribe, Bridgeport Indian Colony, Mono Lake Kutzadikaa, North Fork Rancheria of Mono Indians of California, Picayune Rancheria of the Chukchansi Indians and the Tuolumne Band of Me-Wuk Indians. Their culture and stewardship remain an integral part of the land. As scientists we strive to take responsibility for the impacts of colonialism in our field and move forward with respect and support of indigenous movements and knowledge.

# Conflict of Interest Statement

The authors declare no competing interests.

# Data Availability Statement
The genomic data that support the findings of this study will be submitted to NCBI-GenBank upon acceptance of this manuscript. The manuscript and code for reproducibility is available from corresponding author's GitHub:[https://github.com/jibarozzo/endophyte_mimulus.git](https://github.com/jibarozzo/endophyte_-leaf-traits_mimulus.git)

# References

::: {#refs}
:::

# Figures

```{r, themes}
#| echo: false
#| eval: true
#| tidy: true
#| warning: false

#ggplot themes
# Theme for leaf trait and diversity plots: discrete
theme_lfspp_discrete <-      theme_classic(base_size = 12) +
  theme(legend.title = element_text(face = "bold"),
        strip.text = element_text(face = "bold"),
        panel.border = element_rect(linetype = "blank", fill = NA),
        legend.text = element_text(face = "italic"), #Make legend text italic
        axis.text.x = element_blank())

# Theme for leaf trait plots: continuous
theme_lfspp_continuous <- theme_classic(base_size = 12) +
  theme(legend.title = element_text(face = "bold"),
        panel.border = element_rect(linetype = "blank", fill = NA),
        legend.text = element_text(face = "italic"), #Make legend text italic
        axis.text.x = element_text(size = 10),
        legend.position = "bottom")

# Theme for species by diversity  plots: discrete
theme_sppdiv_discrete <-    theme_classic(base_size = 12) +
   theme(legend.title = element_text(face = "bold"),
         strip.text = element_text(face = "italic", size = 10),
         panel.border = element_rect(linetype = "blank", fill = NA),
         legend.text = element_text(face = "italic", size = 12), #Make legend text italic
         axis.text.x.bottom = element_blank(),
         axis.title.x = element_blank())

# Theme for elevation plots: discrete
theme_elevdiv_discrete <-    theme_classic(base_size = 12) +
   theme(legend.title = element_text(face = "bold"),
         strip.text = element_text(face = "italic", size = 10),
         panel.border = element_rect(linetype = "blank", fill = NA),
         axis.text.x.bottom = element_blank(),
         axis.title.x = element_blank())

# Theme for leaf trait and diversity  plots: continuous
theme_lfdiv_continuous <- theme_classic(base_size = 12) +
   theme(legend.title = element_text(face = "bold"),
         strip.text = element_text(face = "bold.italic", size = 10),
         panel.border = element_rect(linetype = "blank", fill = NA),
        legend.text = element_text(face = "italic"), #Make legend text italic
         axis.text.x = element_text(size = 10))
```


## Figure 1

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 6
#| label: fig-figure1
#| fig-cap: "Map of *Mimulus* spp. sampling sites in the Sierra Nevada, California, USA. A) Map of the Sierra Nevada region showing the location of the three species sampled. The color gradient represents the elevation gradient from low (yellow) to high (purple)."

# Map of sampling sites

#register_google(key ='SECRET_KEY', write = TRUE)

map_latlong <- read_excel("field_data/Mimulus_CH2_Field_Survey.xlsx", sheet = "clean_site_lat_long")

map_latlong <- plant_traits |>
  select(Site, Species, Longitude, Latitude, Elevation_m) |>
  filter(!Site == "FVWB" & !Site == "LYFB" & !Site == "TRTB" & !Species == "M. bicolor")

mappops <- get_map(
  location = c(
    lon =
      mean(map_latlong$Longitude),
    lat =
      mean(map_latlong$Latitude)
  ),
  zoom = 9,
  maptype = "terrain",
  scale = 4
)
map <- ggmap(mappops, extent = "device", darken = 0) +
  geom_point(
    data = map_latlong,
    aes(
      x
      = Longitude,
      y = Latitude,
      size = Elevation_m,
      fill = Species
    ),
    shape = 21,
    alpha = 0.1
  ) +
  scale_fill_manual(
    labels = c("*M. laciniatus*", "*M. nasutus*", "*M. guttatus*"),
    values = met.brewer(
      "Isfahan2",
      n = 4,
      type = "discrete",
      direction = c(-1)
    )
  ) +
  scale_color_manual(
    labels = c("*M. laciniatus*", "*M. nasutus*", "*M. guttatus*"),
    values = met.brewer(
      "Isfahan2",
      n = 4,
      type = "discrete",
      direction = c(-1)
    )
  ) +
  theme_classic(base_size = 12) +
  theme(
    legend.title = element_text(face = "bold"),
    legend.text = element_markdown(size = 10),
    panel.border = element_rect(linetype = "blank", fill = NA),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10)
  ) +
  ylab("Latitude") +
  xlab("Longitude") +
  annotation_north_arrow(
    location = "tr",
    which_north = "true",
    height = unit(0.7, "cm"),
    width = unit(0.7, "cm"),
  ) +
  guides(fill = guide_legend(title = "Species", override.aes = list(alpha = 0.7)),
         size = guide_legend(title = "Elevation (m a.s.l.)", override.aes = list(alpha = 1)))

map
ggsave(
  "figures/figure1.jpg",
  plot = map,
  units = "in",
  width = 6,
  height = 6,
  dpi = 600
)
```

## Figure 2

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 6
#| label: fig-figure2
#| fig-cap: "The Principal Component Analysis (PCA) shows how leaf functional traits are correlated. The PCA analyses shows correlations between ACI, LT, LPS, and LMA as loadings track along PC1 towards positive values, while leaf lobing correlates positively with M. laciniatus. We observed distinct groupings by species along PC2 such that *M. laciniatus*, the most lobed species, is distinct from *M. guttatus* and *M. nasutus*. The latter two have greater overlap along PC1 and PC2."

# PCA using covariates to explain species richness/abundance ####
data.pca <- leaf_traits_noB |>
  na.omit() |>
  select(c(20:24)) #Selecting only the log-transformed leaf trait columns. We use the log transformed data due to our findings of non-normality in the data. See this post for a hot take: https://stats.stackexchange.com/questions/164381/why-log-transforming-the-data-before-performing-principal-component-analysis

leaf_traits.pca <- as.data.frame(na.omit(leaf_traits_noB))


###Run this to create pca with prcomp function
pca <- prcomp(data.pca, scale = TRUE)
pca$rotation=-pca$rotation
pca$x=-pca$x

#Checking the PCA
#plot(pca,type = "lines")
#biplot(pca) # Base type PCA

# PCA using autoplot() and prcomp()and modifying with ggplot syntax ####
pca_auto <- autoplot(pca, data =leaf_traits.pca, 
         alpha = 0, #Setting alpha to zero render the automatic circle point null.Manipulate shapes with geom_point(). 
         loadings = TRUE,
         loadings.colour = "black",
         loadings.label = TRUE, 
         loadings.label.colour="black",
         loadings.label.size = 4, 
         loadings.label.vjust = 0, 
         loadings.label.hjust = -0.15) + 
  geom_point(aes(fill = Species, color = Species), alpha = 0.5, size = 3) +
  geom_hline(yintercept = 0, colour = "gray45") +
  geom_vline(xintercept = 0, colour = "gray45") +
  scale_fill_manual(labels=c("*M. laciniatus*", "*M. nasutus*", "*M. guttatus*"), 
                    values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  scale_color_manual(labels = c("*M. laciniatus*", "*M. nasutus*", "*M. guttatus*"),
                     values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  stat_ellipse(aes(color = Species), 
               geom = "path", 
               linewidth = 1.3, 
               position = "identity", 
               type = "t", 
               linetype = 1, 
               level = 0.95, 
               segments = 51, na.rm = TRUE, show.legend = NA, inherit.aes = TRUE) +
  #scale_x_continuous(expand = c(0.0, -0.09)) +
  theme_classic(base_size = 12) +
  theme(legend.text = element_markdown(),
        legend.position = "bottom",
        legend.title = element_text(face = "bold")) 
  # labs(caption = expression("ACI = anthocyanins, LPS = leaf punch strength, LMA = leaf mass per area, LT = leaf thickness. LBI = Leaf lobe index. All leaf replicates per species ("~ italic("n") ~ "= 501)."))


pca_auto$layers <- c(pca_auto$layers, pca_auto$layers[[2]], pca_auto$layers[[3]]) # This adds/copies layers 2-3 and overlays them. It makes the arrows be on top of the points. There must be a better ways of doing this.

pca_auto
# layer_data(pca_auto,2)
# pca_auto + geom_point(
#   data = layer_data(pca_auto, 2),
#   aes(xend, yend),
#   size = 3
# ) +
#    geom_text(vjust=-1, label=rownames(pca_labels))

ggsave(filename = file.path(path, "figures/figure2.png"), plot = pca_auto, dpi=300, units=c("mm"), width=150, height=150)
```

## Figure 3


<!-- Seems weird to have so much Basidiomycota - how many taxa does this represent? -->

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| message: false
#| fig-width: 6
#| fig-height: 8
#| label: fig-figure3
#| fig-cap: "Relative abundance of fungal phyla in rarefied sequence data by species and elevation category."


#Barplots with microeco
# Quick addition traits log transformed
names_list <- colnames(ps_clean_3@otu_table) #List of samples
ps_clean_3@sam_data <- plant_traits |>
  filter(Unique_ID %in% names_list) |>
  select(c(1:12, 15:20)) |>
  mutate(X = Unique_ID) |>
  column_to_rownames(var = "X") |>
  sample_data()


# Raw data
micro_psclean <- phyloseq2meco(ps_clean_3) # First convert phyloseq to meco


# t1 <- trans_abund$new(dataset = micro_psclean, taxrank = "Phylum", ntaxa = 8)
# t1gg <- t1$plot_bar(others_color = "grey70", facet = c("Species", "Elevation_cat"), xtext_keep = FALSE, legend_text_italic = FALSE, barwidth = 1, color_values = RColorBrewer::brewer.pal(8, "Set3"))
# t1gg <- t1gg +   labs(title = "Raw reads") +
#   theme(legend.title = element_text(face = "bold"),
#          strip.text = element_text(face = "italic", size = 10),
#          panel.border = element_rect(linetype = "blank", fill = NA),
#          axis.text.x = element_blank())
#t1gg


# Rarefied data
micro_rarefied_ps <- phyloseq2meco(rarefied_phyloseq) # First convert phyloseq to meco

r1 <- trans_abund$new(dataset = micro_rarefied_ps, taxrank = "Phylum", ntaxa = 8)
r1gg <- r1$plot_bar(others_color = "grey70", facet = c("Species", "Elevation_cat"), xtext_keep = FALSE, legend_text_italic = FALSE, barwidth = 1, color_values = RColorBrewer::brewer.pal(8, "Set3"))
r1gg <- r1gg +
  labs(title = "Rarefied reads") +
  theme(legend.title = element_text(face = "bold"),
         strip.text = element_text(face = "bold.italic", size = 10),
         panel.border = element_rect(linetype = "blank", fill = NA),
         axis.text.x = element_blank())
r1gg
ggsave(filename = file.path(path, "figures/figure3.png"), plot = r1gg, dpi=300, units=c("mm"), width=300, height=300)

# # Arrange plots
# ra_plots <- ggarrange(t1gg, r1gg,
#                           nrow = 2, ncol = 1,
#                           common.legend = TRUE, legend = "bottom",
#                           labels = c("A", "B"))
# ra_plots
#
# ggsave(file.path(path, "figures/figure2.png") , plot = ra_plots, dpi=300, units=c("mm"), width=200, height=250)
```



## Figure 4

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 7
#| label: fig-figure4
#| fig-cap: "Relationship between Hill orders by host species and elevation. A) Observed ASV richness (𝑞 = 0); B) Shannon's entropy (𝑞 = 1); and C) Inverse Simpson's index (𝑞 = 2) per host species as elevation increases. Blue filled points correspond to *M. laciniatus*, green filled to *M. nasutus*, and yellow to *M. guttatus*.  The solid line represents the linear model fit and the shaded area represents the 95% confidence interval."

# Shannon Diversity or hill number 1
# Three main faceted plots come out of this chunk

########## Labels ###########
#Legend title label
legend_title <- "Species"
# Species labels
labels_list <- c("M. laciniatus", "M. nasutus", "M. guttatus")
names(labels_list) <- c("M. laciniatus", "M. nasutus", "M. guttatus")
#############################


# Hill number 0 
# Species richness or hill number 0

#Linear model summary
#summary(lm(Observed ~ Elevation_m, data = plant_traits_MB))
#richness_p <- format.p(cor.test(plant_traits$Observed, plant_traits$Elevation_m)$p.value) # This is just to make the p-value nicer. 
# Species labels
labels_list <- c("M. laciniatus", "M. nasutus", "M. guttatus")
names(labels_list) <- c("L", "N", "G")

# Plot
hill0_diversity <- ggplot(data = plant_traits_MB, 
                          aes(x = Elevation_m, y = Observed, color = Species)) +
   geom_jitter(aes(color = Species, shape = Site), 
               #stat = "summary",
               #fun = mean,
               size = 3, alpha = 0.3)+ 
  geom_smooth(method = lm, se = F) +
  geom_smooth(method = lm, se = T, level = 0.95, na.rm = F, color = "black") +
  stat_poly_eq(formula = y ~ x, color = "black",
               aes(label = paste(#..eq.label.., 
                 ..adj.rr.label.., #"~italic(p) ==",
                 ..p.value.label.., 
                 #signif(..p.value.., digits = 4),
                 sep = "*`,`~~")), 
               rr.digits = 2,
               p.digits = 3 ,
               parse = TRUE,
               label.x = "right",
               label.y= "top",
               vstep = 0.08) +
  stat_poly_eq(
      formula = y ~ x,
      aes(
        color = Species,
        label = paste(#..eq.label..,
          ..adj.rr.label..,
          #"~italic(p) ==",
          ..p.value.label..,
          #signif(..p.value.., digits = 4),
          sep = "*`,`~~")
      ),
      parse = TRUE,
      label.x.npc = "center",
      vstep = 0.08
    ) +
# stat_poly_eq(formula = y ~ x, color = aes(color=Species),
#                aes(label = paste(#..eq.label.., 
#                  ..adj.rr.label.., #"~italic(p) ==",
#                  ..p.value.label.., 
#                  #signif(..p.value.., digits = 4),
#                  sep = "*`,`~~")), 
#                rr.digits = 2,
#                p.digits = 3 ,
#                parse = TRUE,
#                label.x = "left",
#                label.y= "top",
#                vstep = 0.08) +
   scale_fill_manual(values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) + 
   scale_color_manual(legend_title, values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) + 
 scale_shape_manual(values = c(1:22)) +
   theme_lfdiv_continuous +
   labs(y =  "Alpha diversity (\U1D492 = 0)", x = "", col = "") + 
   scale_y_continuous(expand = c(0.05, 0.05)) +
  scale_x_continuous(expand = c(0.05, 0.05)) 
  #facet_wrap(~ Species) +
  #guides(color = "none")        #To turn off color legend
  # annotate(geom = "text", x = 2100, y = 2.5, label = "p = 0.03", color = "black")
#annotate(geom="text", x=2100, y=50, label="p = 0.06", color = "black")

#hill0_diversity

#ggsave(filename = file.path(path, "figures/richness_elevation.png"), plot = hill0_diversity, dpi=300, units=c("mm"), width=250, height=150)

# Linear models per species
# Linear model summary
# M. laciniatus
#summary(lm(Observed ~ Elevation_m, data = plant_traits[plant_traits$Species == "L",]))
#Not significant

# M. nasutus
#summary(lm(Observed ~ Elevation_m, data = plant_traits[plant_traits$Species == "N",]))
#Not significant

# M. guttatus
#summary(lm(Observed~ Elevation_m, data = plant_traits[plant_traits$Species == "G",]))
#Not significant

# Hill number 1
# Shannon diversity or hill number 1
# Linear model summary
#summary(lm(hill_1 ~ Elevation_m, data = plant_traits_MB))
#elev_div_p <- format.p(cor.test(plant_traits$Shannon, plant_traits$Elevation_m)$p.value) # This is just to make the p-value nicer. 

# Plot
hill1_diversity <-ggplot(data = plant_traits_MB, 
  aes(x = Elevation_m, y = hill_2, color = Species)) +
   geom_jitter(aes(color = Species,  shape = Site),
               #stat = "summary",
               #fun = mean,
               size = 3, alpha = 0.3) + 
   geom_smooth(method = lm, se = F) + 
   geom_smooth(method = lm, se = T, level = 0.95, na.rm = F, color = "black") + 
  stat_poly_eq(formula = y ~ x, color = "black",
               aes(label = paste(#..eq.label.., 
                 ..adj.rr.label.., #"~italic(p) ==",
                 ..p.value.label.., 
                 #signif(..p.value.., digits = 4),
                 sep = "*`,`~~")), 
               rr.digits = 2,
               p.digits = 3 ,
               parse = TRUE,
               label.x = "right",
               label.y = "top",
               vstep = 0.08) +
  stat_poly_eq(
        formula = y ~ x,
        aes(
          color = Species,
          label = paste(#..eq.label..,
            ..adj.rr.label..,
            #"~italic(p) ==",
            ..p.value.label..,
            #signif(..p.value.., digits = 4),
            sep = "*`,`~~")
        ),
        parse = TRUE,
        label.x.npc = "center",
        vstep = 0.08
      ) +
   scale_color_manual(legend_title,
                     values = met.brewer("Isfahan2", n = 4, type = "discrete", 
                                          direction = c(-1))) + 
 scale_shape_manual(values = c(1:22)) +
  labs(y = "Alpha diversity (\U1D492 = 1)", x = "", col = "") +
  scale_y_continuous(expand = c(0.05, 0.05)) +
  theme_lfdiv_continuous +
  #theme(axis.title.x = element_text()) +
  theme(strip.text = element_blank()) 
  #facet_wrap(~ Species) +
  #guides(color = "none")        #To turn off color legend

#hill1_diversity

#ggsave(filename = file.path(path, "figures/shannon_elevation.png"), plot = hill1_diversity, dpi=300, units=c("mm"), width=250, height=150)

# Linear models per species
# Linear model summary
# M. laciniatus
#summary(lm(hill_1 ~ Elevation_m, data = plant_traits[plant_traits$Species == "L",]))
#Not significant

# M. nasutus
#summary(lm(hill_1 ~ Elevation_m, data = plant_traits[plant_traits$Species == "N",]))
#Not significant

# M. guttatus
#summary(lm(hill_1 ~ Elevation_m, data = plant_traits[plant_traits$Species == "G",]))
#Not significant

# Hill number 2
# Inverse Simpson diversity or hill number 2
# Linear model summary
#summary(lm(hill_2 ~ Elevation_m, data = plant_traits))
#elev_div_p <- format.p(cor.test(plant_traits$Shannon, plant_traits$Elevation_m)$p.value) # This is just to make the p-value nicer. 

# Plot
hill2_diversity <- ggplot(data = plant_traits_MB, 
  aes(x = Elevation_m, y = hill_2, color = Species)) +
   geom_jitter(aes(color = Species,  shape = Site), 
               #stat = "summary",
               #fun = mean,
               size = 3, alpha = 0.3) + 
   geom_smooth(method = lm, se = F) + 
   geom_smooth(method = lm, se = T, level = 0.95, na.rm = F, color = "black") + 
  stat_poly_eq(formula = y ~ x, color = "black",
               aes(label = paste(#..eq.label.., 
                 ..adj.rr.label.., #"~italic(p) ==",
                 ..p.value.label.., 
                 #signif(..p.value.., digits = 4),
                 sep = "*`,`~~")), 
               rr.digits = 2,
               p.digits = 3 ,
               parse = TRUE,
               label.x = "right",
               label.y = "top",
               vstep = 0.08) +
  stat_poly_eq(
        formula = y ~ x,
        aes(
          color = Species,
          label = paste(#..eq.label..,
            ..adj.rr.label..,
            #"~italic(p) ==",
            ..p.value.label..,
            #signif(..p.value.., digits = 4),
            sep = "*`,`~~")
        ),
        parse = TRUE,
        label.x.npc = "center",
        vstep = 0.08
      ) +
   #scale_fill_manual(labels = c("*M. laciniatus*","*M. nasutus*", "*M. guttatus*"), 
                     #values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) + 
   scale_color_manual(legend_title, labels = c("M. laciniatus", "M. nasutus", "M. guttatus"), 
                     values = met.brewer("Isfahan2", n = 4, type = "discrete", 
                                          direction = c(-1))) + 
 scale_shape_manual(values = c(1:22)) +
  labs(y = "Alpha diversity (\U1D492 = 2)", x = "Elevation (m.a.s.l.)", col = "") +
  theme_lfdiv_continuous +
  theme(strip.text = element_blank()) 
 # facet_wrap(~ Species) 
  #guides(color = "none")        #To turn off color legend

#hill2_diversity

#ggsave(filename = file.path(path, "figures/invsimpson_elevation.png"), plot = hill2_diversity, dpi=300, units=c("mm"), width=250, height=150)

# Linear models per species
# Linear model summary
# M. laciniatus
#summary(lm(hill_2 ~ Elevation_m, data = plant_traits[plant_traits$Species == "L",]))
#Not significant

# M. nasutus
#summary(lm(hill_2 ~ Elevation_m, data = plant_traits[plant_traits$Species == "N",]))
#Not significant

# M. guttatus
#summary(lm(hill_2 ~ Elevation_m, data = plant_traits[plant_traits$Species == "G",]))
#Not significant


###### Joined plots ######
hill_numbers <-
  ggarrange(
    hill0_diversity,
    hill1_diversity,
    hill2_diversity,
    common.legend = TRUE,
    legend = "right",
    labels = c("A", "B", "C"),
    nrow = 3,
    ncol = 1
  )

hill_numbers
ggsave(filename = file.path(path, "figures/figure4.png"), plot = hill_numbers, dpi=300, units=c("mm"), width=275, height=300)
```



## Figure 5

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 5
#| fig-height: 6
#| label: fig-figure5
#| fig-cap: "FEF community composition association with leaf functional traits and elevation by host species. Distance-based Redundancy Analysis (dbRDA) plot of rarefied FEF community and leaf functional traits by species facetted by elevation category. Each cluster of points represents a rarefied FEF community sample from one host species. Solid arrow lines represent significant associations (*p* < .01)with non-significant leaf trait associations left out for clarity. The length and direction of the arrows indicate the strength and direction of the association between the traits and the FEF community composition. Ellipses represent 95% confidence intervals. The plot is based on the Bray-Curtis dissimilarity matrix."

# dbRDA using Elevation as continuous variable
#logACI + logLT + logLPS + logLMA +
m5_hell <- readRDS("clean_data/statistics/m5_hell.rds")

# Final model
final_model <- m5_hell # Previously m4_hell was the best

B <- summary(final_model)

# Extracting scores for the plot
ordination_scores <- as.data.frame(vegan::scores(final_model)$sites)

# Extracting names for the plot
ordination_names <- rownames(ordination_scores)

# Row names to column
ordination_scores$Unique_ID <- rownames(ordination_scores)

# Subset for the plot
dbrda_df <- plant_traits_MB
  #filter(Unique_ID %in% ordination_names)

# Merge scores with trait data
merged_data <- right_join(ordination_scores, dbrda_df, by = c("Unique_ID" = "X"))

# Extracting scores for arrows
arrow_scores <- as.data.frame(vegan::scores(final_model, display = "bp"))

# Define the arrow aesthetic mapping
arrow_map <- aes(
  xend = (dbRDA1 * 2),
  yend = (dbRDA2 * 2),
  x = 0,
  y = 0,
  shape = NULL,
  color = NULL,
  fill = NULL
)
label_map <- aes(
  x = 2 * dbRDA1,
  y = 2 * dbRDA2,
  label = row.names(arrow_scores),
  shape = NULL,
  color = NULL,
  fill = NULL
)
arrowhead = arrow(length = unit(0.02, "npc"))

#subset A4 for labeling
arrow_scores <- arrow_scores[order(rownames(arrow_scores)), ]
A4 <- lapply(1:6, function(i)
  arrow_scores[i, ])

# dbRDA plot based on species
p_final_model <-
  ggplot(data = merged_data, aes(x = dbRDA1, y = dbRDA2)) +
  geom_point(
    data = merged_data,
    aes(color = Species, shape  = Elevation_cat),
    size = rel(4),
    alpha = 0.2,
    stroke = 1
  ) +
  stat_ellipse(
    aes(color = Species),
    geom = "path",
    linewidth = 1.3,
    position = "identity",
    type = "t",
    linetype = 1,
    level = 0.95,
    segments = 51,
    na.rm = FALSE,
    show.legend = NA,
    inherit.aes = TRUE
  ) +
  geom_segment(
    arrow_map,
    linewidth = rel(.8),
    data = arrow_scores[1:2, ],
    color = "black",
    arrow = arrowhead
  ) +
  geom_vline(xintercept = 0,
             color = "grey30",
             linetype = "dashed") + # plot vertical line
  geom_hline(yintercept = 0,
             color = "grey30",
             linetype = "dashed") + # plot horizontal line
  geom_text(
    label_map,
    size = rel(4),
    data = A4[[1]],
    fontface = "bold",
    show.legend = FALSE,
    label = c("Elev."),
    nudge_x = 0.45
  ) +
  geom_text(
    label_map,
    size = rel(4),
    data = A4[[2]],
    fontface = "bold",
    show.legend = FALSE,
    label = c("logLBI"),
    nudge_y = 0.2
  ) +
  xlab(label = paste(
    "dbRDA1 (",
    round(B$concont$importance[2, 1] * 100, digits = 1),
    "%)",
    sep = ""
  )) +
  ylab(label = paste(
    "dbRDA2 (",
    round(B$concont$importance[2, 2] * 100, digits = 1),
    "%)",
    sep = ""
  )) +
  scale_color_manual(
    labels = c("*M. laciniatus*", "*M. nasutus*", "*M. guttatus*"),
    values = met.brewer(
      "Isfahan2",
      n = 3,
      type = "discrete",
      direction = -1
    )
  ) +
  scale_shape_manual(values = c(1:3)) +
  theme_classic(base_size = 12) +
  theme(
    legend.title = element_text(face = "bold"),
    legend.text = element_markdown(),
    legend.position = "bottom",
    plot.caption = element_markdown(),
    strip.text = element_text(face = "bold")
  ) +
  labs(caption = "PERMDISP: *F*<sub>7,4192</sub> = 268.1, *p* = .001")  #Results from beta_perm3
  #facet_wrap(~ Elevation_cat) 
     
# Rearranging the order of layers so arrows are above the points.
#p_final_model$layers
# p_final_model$layers <- c(p_final_model$layers, p_final_model$layers[[2]], p_final_model$layers[[5]], p_final_model$layers[[6]])
p_final_model

#ggsave(filename=file.path(path, "figures/dbRDA_elevation_labels.svg"), plot = p_final_model, dpi=300, units=c("mm"), width=250, height=150)


#Save the plot
ggsave(file.path(path, "figures/figure5.png"), plot = p_final_model, dpi=300, units=c("mm"), width=150, height=150)
```

## Figure 6

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 7
#| label: fig-figure6
#| fig-cap: "Relationship between β-diversity and elevation by host species. The solid line represents simple linear regression and the shaded area represents the 95% confidence interval."

# Mean-Bray-Curtis dissimilarity
mean_bray <- vegdist(dbrda_hell_matrix, method = "bray") |>
  as.matrix() |>
  as.data.frame()

mean_bray <- rownames_to_column(mean_bray, var = "Unique_ID") |>
  mutate(mean_bray = rowMeans(mean_bray)) |>
  select(Unique_ID, mean_bray)

# Joining to dbrda_traits to avoid duplicate DFs 
dbrda_traits_imp <- left_join(dbrda_traits_imp, mean_bray, by = c("X"= "Unique_ID"))

# Simple linear regressions for the stats
# summary(lm(mean_bray ~ Elevation_m, data =subset(dbrda_traits_imp, Species == "M. laciniatus")))
# summary(lm(mean_bray ~ Elevation_m, data =subset(dbrda_traits_imp, Species == "M. nasutus")))
# summary(lm(mean_bray ~ Elevation_m, data =subset(dbrda_traits_imp, Species == "M. guttatus")))

# GGplot
beta_elevation <- ggplot(data = dbrda_traits_imp, aes(x = Elevation_m, y = mean_bray, color = Species)) +
  geom_point(
    aes(color = Species),
    size = rel(4),
    alpha = 0.2,
    stroke = 1
  ) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x =  "Elevation (m)", y = "β-diversity") +
  scale_color_manual(
    labels = c("*M. laciniatus*", "*M. nasutus*", "*M. guttatus*"),
    values = met.brewer(
      "Isfahan2",
      n = 3,
      type = "discrete",
      direction = -1
    )
  ) +
  stat_poly_eq(formula = y ~ x,
               color = "black",
               aes(
                   label = paste(#..eq.label.., 
                 after_stat(adj.rr.label), 
                 #"~italic(p) ==",
                 ..p.value.label.., 
                 #signif(..p.value.., digits = 4),
                 sep = "*`,`~~")),
               #method = "loess",
               rr.digits = 2,
               p.digits = 3 ,
               parse = TRUE,
               label.x = "right",
               label.y= "top",
               vstep = 0.08) +
  theme_classic(base_size = 16) +
  theme(
    legend.title = element_text(face = "bold"),
    legend.text = element_markdown(),
    legend.position = "bottom",
    plot.caption = element_markdown(),
    strip.text = element_text(face = "bold.italic", size = 14)
  ) +
  facet_wrap( ~ Species)

beta_elevation

ggsave(file.path(path, "figures/figure6.png"), plot = beta_elevation, dpi=300, units=c("mm"), width=250, height=250)

```


## Figure 7

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6 
#| fig-height: 7
#| label: fig-figure7
#| fig-cap: "Dissimilarity of FEF communities correlates with *Mimulus* as a function of geographic distance per site. Data represent pairwise Bray-Curtis dissimilarity (β-diversity) among sites. A) FEF community dissimilarity in considering all host species. B) FEF community dissimilarity by host species. Blue filled points correspond to *M. laciniatus*, green filled to *M. nasutus*, and yellow to *M. guttatus*."

# Prepare the data for correlation plot
# Longer format
asv_vegdist <- as.matrix(vegdist(t(asv_matrix), method = "bray"))

asv_long <-  asv_vegdist |>
  as_tibble(rownames = "A") |>
  pivot_longer(-A, names_to = "B", values_to = "bray_diss")

geo_long <- geo_distm |>
  as_tibble(rownames = "A") |>
  pivot_longer(-A, names_to = "B", values_to = "geo_dist")

# Merge the two dataframes
asv_geo_matrix_long <-
  inner_join(asv_long, geo_long, by = c("A", "B")) |>
  inner_join(
    plant_traits |>
      dplyr::select(Unique_ID, Site, Species, Elevation_cat),
    by = c("A" = "Unique_ID")
  )

# Mantel test
# All data points
set.seed(859)
abund_geo  <- mantel(
  asv_vegdist,
  geo_distm,
  method = "spearman",
  parallel = 8,
  permutations = 999,
  na.rm = TRUE
)
#abund_geo

# Mantel test by species
# M. laciniatus

# Samples with species to subset/filter matrices downstream
# M. laciniatus
lac_names_rfy <- plant_traits |>
  slice(1:157) |> # Only selecting observations with ASV data
  filter(Unique_ID %in% final_names_methodB,
         Species == "M. laciniatus") |>
  pull(Unique_ID)

# Subset the data
lac_mantel  <- mantel(
  asv_long |> # Subset the data
    filter(A %in% lac_names_rfy, B %in% lac_names_rfy) |>
    pivot_wider(names_from = B, values_from = bray_diss) |>
    column_to_rownames(var = "A") |>
    as.matrix(),
  geo_long |> # Subset the data
    filter(A %in% lac_names_rfy, B %in% lac_names_rfy) |>
    pivot_wider(names_from = B, values_from = geo_dist) |>
    column_to_rownames(var = "A") |>
    as.matrix(),
  method = "spearman",
  permutations = 999,
  na.rm = TRUE
)
#lac_mantel

# M. nasutus
nas_names_rfy <- plant_traits |>
  slice(1:157) |> # Only selecting observations with ASV data
  filter(Unique_ID %in% final_names_methodB,
         Species == "M. nasutus") |>
  pull(Unique_ID)
nas_mantel  <- mantel(
  asv_long |> # Subset the data
    filter(A %in% nas_names_rfy, B %in% nas_names_rfy) |>
    pivot_wider(names_from = B, values_from = bray_diss) |>
    column_to_rownames(var = "A") |>
    as.matrix(),
  geo_long |> # Subset the data
    filter(A %in% nas_names_rfy, B %in% nas_names_rfy) |>
    pivot_wider(names_from = B, values_from = geo_dist) |>
    column_to_rownames(var = "A") |>
    as.matrix(),
  method = "spearman",
  permutations = 999,
  na.rm = TRUE
)
#nas_mantel

# M. guttatus
gut_names_rfy <- plant_traits |>
  slice(1:157) |> # Only selecting observations with ASV data
  filter(Unique_ID %in% final_names_methodB,
         Species == "M. guttatus") |>
  pull(Unique_ID)
# Subset the data
gut_mantel  <- mantel(
  asv_long |> # Subset the data
    filter(A %in% gut_names_rfy, B %in% gut_names_rfy) |>
    pivot_wider(names_from = B, values_from = bray_diss) |>
    column_to_rownames(var = "A") |>
    as.matrix(),
  geo_long |> # Subset the data
    filter(A %in% gut_names_rfy, B %in% gut_names_rfy) |>
    pivot_wider(names_from = B, values_from = geo_dist) |>
    column_to_rownames(var = "A") |>
    as.matrix(),
  method = "spearman",
  permutations = 999,
  na.rm = TRUE
)
#gut_mantel

# Visualize the correlation
# Species labels
labels_list <- c("M. laciniatus", "M. nasutus", "M. guttatus")
names(labels_list) <-
  c("M. laciniatus", "M. nasutus", "M. guttatus")

# Correlation ggplot
mantel_test <-
  ggplot(asv_geo_matrix_long, aes(y = bray_diss, x = geo_dist)) +
  geom_point(aes(color = Species), alpha = 0.5, size = 3) +
  geom_smooth(method = "lm",
              se = TRUE,
              color = "black") +
  scale_color_manual(values = met.brewer(
    "Isfahan2",
    n = 4,
    type = "discrete",
    direction = c(-1)
  )) +
  labs(y = "β-diversity", x = "Geographic distance (m)", col = "") +
  scale_x_continuous(
    labels = function(x)
      format(x, scientific = FALSE)
  ) +
  theme_lfspp_continuous

# Adding the mantel test results
italic_r <- expression(italic("r ="))
italic_p <- expression(italic("p ="))
statistic <- round(abund_geo$statistic, digits = 2)
signif <- round(abund_geo$signif, digits = 3)


mantel_test <- mantel_test + ggplot2::annotate(
  "text",
  x = 110000,
  y = 0,
  label = paste(italic_r, statistic, italic_p, signif, sep = "~"),
  size = 5 ,
  parse = TRUE
)
#mantel_test

# By species
mantel_species <-
  ggplot(asv_geo_matrix_long, aes(y = bray_diss, x = geo_dist)) +
  geom_jitter(aes(color = Species), alpha = 0.5, size = 3) +
  geom_smooth(method = "lm",
              se = TRUE,
              color = "black") +
  scale_color_manual(legend_title, values = met.brewer(
    "Isfahan2",
    n = 4,
    type = "discrete",
    direction = c(-1)
  )) +
  labs(y = "β-diversity", x = "Geographic distance (m)", col = "") +
  scale_x_continuous(
    labels = function(x)
      format(x, scientific = FALSE)
  ) +
  theme_lfspp_continuous +
  theme_classic(base_size = 16) +
  theme(
    legend.title = element_text(face = "bold"),
    panel.border = element_rect(linetype = "blank", fill = NA),
    legend.text = element_text(face = "italic"),
    strip.text = element_text(face = "bold.italic", size = 14),
    axis.text.x = element_text(size = 10),
    legend.position = "bottom"
  ) +
  guides(guide_legend(color = "Species")) +
  facet_wrap( ~ Species)
#mantel_species

# Adding the mantel test results

italic_r <- expression(italic("r ="))
italic_p <- expression(italic("p ="))
statistic_lac  <- round(lac_mantel$statistic, digits = 2)
signif_lac <- round(lac_mantel$signif, digits = 3)
statistic_nas  <- round(nas_mantel$statistic, digits = 2)
signif_nas <- round(nas_mantel$signif, digits = 3)
statistic_gut  <- round(gut_mantel$statistic, digits = 3)
signif_gut <- round(gut_mantel$signif, digits = 3)

my_tags <-
  c(
    paste0("r = ", statistic_lac, ", p = ", signif_lac, sep = " "),
    paste0("r = ", statistic_nas, ", p = ", signif_nas, sep = " "),
    paste0("r = ", statistic_gut, ", p = ", signif_gut, sep = " ")
  )

mantel_species <- tag_fat(
  mantel_species,
  x = -Inf,
  y = -Inf,
  vjust = -1,
  hjust = -0.95,
  open = "",
  close = "",
  fontface = 3,
  size = 4,
  tag_pool = my_tags
)


# Joining plots
# mantel <- ggpubr::ggarrange(mantel_test, "", mantel_species,
#                     ncol = 1, nrow = 3,
#                     widths = c(1, 0.1, 1),
#                     heights = c(1, 0.1, 1),
#                     labels = c("A", "B"),
#                     common.legend = T, legend = "bottom")

mantel_species
ggsave(file.path(path, "figures/figure7.png") , plot = mantel_species, dpi=300, units=c("mm"), width=250, height=200)
```


## Table 1

```{r, table1}
#| echo: false
#| eval: true
#| tidy: true
#| tbl-cap: "Generalized linear mixed models (GLMMs) for predicting mean Bray-Curtis dissimilarity"

load("clean_data/statistics/glmm_models.rda")

# Huxtable to flextable
table1 <- huxreg(
  "Model 1" = B1REM,
  "Model 2" = B3REM,
  "Model 3" = B4REM,
  "Model 4" = B5REM,
  coefs = c(
    "Intercept" = "(Intercept)",
    "M. nasutus" = "SpeciesM. nasutus",
  "M. guttatus" = "SpeciesM. guttatus",
    "Elevation" = "Elevation_m",
    "ACI" = "logACI",
    "LT" = "logLT",
    "LPS" = "logLPS",
    "LMA" = "logLMA",
    "LBI" = "logLBI"
  ),
  error_pos = "right",
  error_format = "({conf.low} - {conf.high} ) [{statistic}]",
  number_format = "%.2f",
  align = ".",
  ci_level = 0.95,
  stars = c(`***` = 0.001, `**` = 0.01, `*` = 0.05),
  bold_signif = 0.05,
  borders = 0.4,
  outer_borders = 0.8,
  note = "Significance levels are represented by asterisks [p < 0.05 (*), p < 0.01 (**), p < 0.001 (***), and p < 0.0001 (****)]. T statistics in brackets.",
  statistics = c("N" = "nobs", "logLik", "AIC")
) |>
  insert_row("Predictors",
             "Estimate",
             "CI [t-statistic]",
             "Estimate",
             "CI [t-statistic]",
             "Estimate",
             "CI [t-statistic]",
             "Estimate",
             "CI [t-statistic]",
             after = 1) |>
  insert_row("Random effects", "", "", "", "", "", "", "", "", after = 11) |>
  set_all_padding(10) |>
  set_outer_padding(0) |>
  set_bold(row = 1, col = everywhere) |>
  #set_bottom_border(row = 1, col = everywhere) |>
  set_width(1.1) |>
  theme_article() |> #flextable::theme_apa() is good too.
  huxtable::as_flextable(column_to_header = TRUE) |>
  hline(12,
        part = "all",
        border = fp_border(
          color = "black",
          style = "solid",
          width = 1
        )) |> 
  hline(15,
        part = "all",
        border = fp_border(
          color = "black",
          style = "solid",
          width = 1.25
        )) |>
  hline(
    i = 1,
    j = c(2:9),
    part = "body",
    border = fp_border(
      color = "black",
      style = "solid",
      width = 1
    )
  ) |>
  flextable::align(
    i = 1:14,
    j = 2:9,
    align = "center",
    part = "body"
  ) |>
  flextable::style(i = c(1:16),
                   j = 1,
                   pr_t = fp_text(bold = FALSE)) |> 
  flextable::italic(
    i = c(4, 5, 12, 15),
    j = 1,
    italic = TRUE,
    part = "body"
  ) |>
  flextable::bold(
    i = c(2, 12),
    j = 1,
    bold = TRUE,
    part = "body"
  ) |>
  flextable::set_caption("Generalized linear mixed models (GLMMs) for predicting mean Bray-Curtis dissimilarity") |>
  flextable::font(
    i = NULL,
    j = NULL,
    fontname = "Tex Gyre Termes",
    part = "all",
    cs.family = fontname,
    hansi.family = fontname,
    eastasia.family = fontname
  )
  # add_footer_lines(
  #   value = as_paragraph("Significance levels are represented by asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]. T statistics in brackets.")
  # )
table1

#save_as_docx(table1, path = "tables/model_comparison.docx")
#save_as_image(table1, path = "tables/model_comparison.svg")
```


# Supplementary Material

## Figure S1

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| message: false
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM21
#| fig-cap: "Figure S1. Correlation matrix of log-transformed leaf functional traits by host species. The plot is based on the Spearman's rho. Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < 0.01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

t4 <- trans_env$new(dataset = micro_psclean, env_cols = c(5,14:18))

t4gg <- t4$cal_autocor(group = "Species", upper = list(continuous = GGally::wrap("cor", method= "spearman"))) + 
  theme(
    legend.title = element_text(face = "bold"),
         strip.text = element_text(face = "bold", size = 10),
         panel.border = element_rect(linetype = "blank", fill = NA)) +
  scale_fill_manual(values = met.brewer("Isfahan2", n = 3, type = "discrete", direction = c(-1))) +
  scale_color_manual(values = met.brewer("Isfahan2", n = 3, type = "discrete", direction = c(-1)))

t4gg
ggplot2::ggsave(file.path(path, "figures/supp_fig1.jpg") , plot = t4gg, dpi=300, units=c("mm"), width=300, height=300)
```

## Figure S2

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| message: false
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM22
#| fig-cap: "Figure S2. Heatmap of the top 10 most abundant FEF genera in the dataset. The heatmap is faceted by host species and elevation category."

# Plot with micro_psclean
t2 <-
  trans_abund$new(dataset = micro_psclean,
                  taxrank = "Genus",
                  ntaxa = 10)
t2gg <-
  t2$plot_heatmap(
    facet = c("Species", "Elevation_cat"),
    xtext_keep = FALSE,
    withmargin = FALSE
  ) +
  labs(title = "Un-rarefied reads") +
  theme(
    legend.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold.italic", size = 12),
    panel.border = element_rect(linetype = "blank", fill = NA),
    axis.text.x = element_blank()
  )

#t2gg



r2 <-
  trans_abund$new(dataset = micro_rarefied_ps,
                  taxrank = "Genus",
                  ntaxa = 10)
r2gg <-
  r2$plot_heatmap(
    facet = c("Species", "Elevation_cat"),
    xtext_keep = FALSE,
    withmargin = FALSE
  ) +
  labs(title = "Rarefied reads") +
  theme(
    legend.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold.italic", size = 12),
    panel.border = element_rect(linetype = "blank", fill = NA),
    axis.text.x = element_blank()
  )

r2gg
ggsave(filename = file.path(path, "figures/supp_fig2b.jpg"), plot = r2gg, dpi=300, units=c("mm"), width=300, height=250)

# Joined heatplots
heatplots <- ggarrange(t2gg, r2gg, ncol = 1, nrow = 2, 
                           common.legend = TRUE, 
                           legend = "bottom",
                           labels = c("A", "B"),
                           heights = c(1, 1))
#heatplots

# Save the plot
ggplot2::ggsave(filename = file.path(path, "figures/supp_fig2.png"), plot = heatplots, dpi=300, units=c("mm"), width=300, height=250)
```


## Figure S3
```{r}
#| echo: false
#| eval: true
#| tidy: true
#| message: false
#| fig-height: 6
#| fig-width: 6
#| label: suppfigMIM23
#| fig-cap: "Figure S3. Venn diagrams of presence and absence of ASVs among species and elevation categories un-rarified data and rarified data. A) Represents the overlap of ASVs in in host species and B) per elevation category. Number of ASVs are represented the count and in parentheses as a ratio shared between host species and elevation category"

# un-rarified data
# Per Species
dataset1 <- micro_psclean$merge_samples(use_group = "Species")
# create trans_venn object
v1 <-
  trans_venn$new(dataset1, ratio = "numratio") # Percentage represent ratio of shared ASVs
v1gg <-
  v1$plot_venn() + #https://r-charts.com/color-palettes/#discrete
  ggtitle("Un-rarified data") +
  theme(plot.title = element_text(
    face = "bold",
    size = 16,
    hjust = 0.5,
    vjust = 0.5
  ))


# Per Elevations
dataset2 <- micro_psclean$merge_samples(use_group = "Elevation_cat")
v2 <- trans_venn$new(dataset2, ratio = "numratio") 
#v2$plot_venn()

# rarified data
#micro_rarefied_ps <- phyloseq2meco(rarefied_phyloseq) # First convert phyloseq to meco
# Per Species
dataset3 <- micro_rarefied_ps$merge_samples(use_group = "Species") 
# create trans_venn object
v3 <- trans_venn$new(dataset3, ratio = "numratio")
v3gg <- v3$plot_venn() +
  #ggtitle("Rarified data") +
  theme(plot.title = element_text(
    face = "bold",
    size = 16,
    hjust = 0.5,
    vjust = 0.5
  ))

# Per Elevations
dataset4 <- micro_rarefied_ps$merge_samples(use_group = "Elevation_cat")
v4 <- trans_venn$new(dataset4, ratio = "numratio")
#v4$plot_venn()

#                                         
# # joined plots
# venn_plots <-
#   ggarrange(
#     v1gg,
#     v3gg,
#     "", # Fake plot to create space
#     "",
#     v2$plot_venn(),
#     v4$plot_venn(),
#     ncol = 2,
#     nrow = 3,
#     common.legend = TRUE,
#     heights = c(1, 0.1, 1),
#     legend = "bottom",
#     labels = c("A", "B","", "", "C", "D")
#   )

# venn_plots

rrfy_venn <-
  ggarrange(
    v3gg,
    v4$plot_venn(),
    ncol = 2,
    nrow = 1,
    common.legend = TRUE,
    legend = "bottom",
    labels = c("A", "B")
  ) + ggtitle("Rarified data")

rrfy_venn

ggsave(filename = file.path(path, "figures/supp_fig3.jpg"), plot = rrfy_venn, dpi=300, units=c("mm"), width=250, height=225)

# # Save
# ggsave(filename = file.path(path, "figures/supp_fig3.png"), plot = venn_plots, dpi=300, units=c("mm"), width=300, height=300)

```

## Figure S4

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| warning: false
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM24
#| fig-cap: "Figure S4. Alpha and beta diversity means comparisons in host species and elevation categories. Alpha diversity means comparisons; A) Observed ASV richness (𝑞 = 0); B) Shannon's entropy (𝑞 = 1); and C) Inverse Simpson's index (𝑞 = 2) per host species. A-C) Blue filled boxplots correspond to *M. laciniatus*, green filled to *M. nasutus*, and yellow to *M. guttatus*. Beta diversity mean comparisons; A) Observed ASV richness (𝑞 = 0); B) Shannon's entropy (𝑞 = 1); and C) Inverse Simpson's index (𝑞 = 2) per elevation category. D-F) Violet boxplots correspond to LOW elevation sites, pink filled to MID and light maroon to HIGH elevation sites, while squares represent *M. laciniatus*, circles *M. nasutus* and triangles *M. guttatus*.  Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

############################################
#### Alpha diversity ####
#### At the species level ####
############################################

labels_list <- c("M. laciniatus", "M. nasutus", "M. guttatus")
names(labels_list) <- c("M. laciniatus", "M. nasutus", "M. guttatus")

# Plot
hill0_species <- ggplot(data = plant_traits_MB, 
  aes(x = Species, y = Observed)) +
    geom_point(
  aes(color = Species),
  alpha = 0.2,
  size = 1.7,
  position = position_jitter(w = 0.2)
) +
  geom_boxplot(
    aes(color = Species),
    alpha = 0.5,
    outlier.alpha = 0.0) +
  scale_color_manual(name = "Species", values= met.brewer("Isfahan2", n = 4, type = "discrete", direction = -1)) +
  theme_sppdiv_discrete +
  labs(y =  "Alpha diversity (\U1D492 = 0)", x = "", col = "") + 
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test", 
                     label = "p.signif") 
  #guides(color = "none")        #To turn off color legend
  
#hill0_species
# annotate(geom = "text", x = 2100, y = 2.5, label = "p = 0.03", color = "black")
#annotate(geom="text", x=2100, y=50, label="p = 0.06", color = "black")


# Hill number 1
# Shannon diversity or hill number 1
# Plot
hill1_species <- ggplot(data = plant_traits_MB, 
  aes(x = Species, y = hill_1)) +
  geom_point(
  aes(color = Species),
  alpha = 0.2,
  size = 1.7,
  position = position_jitter(w = 0.2)
) +
  geom_boxplot(
    aes(color = Species),
    alpha = 0.5,
    outlier.alpha = 0.0) +
  scale_color_manual(name = "Species", values= met.brewer("Isfahan2", n = 4, type = "discrete", direction = -1)) +
  #scale_color_manual(aes(color = "black")) +
  theme_sppdiv_discrete +
  theme(strip.text = element_text(face = "bold", size = 14)) +
  #labs(y =  "Alpha and Beta diversity (\U1D492 = 1)", x = "", col = "") + 
  labs(y =  "Alpha diversity (\U1D492 = 1)", x = "", col = "") +
  scale_y_continuous(expand = c(0.05, 0.05)) +
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test", 
                     label = "p.signif")
 #facet_wrap(~ Elevation_cat)
 # guides(color = "none")        #To turn off color legend
  
#hill1_species


# Hill number 2
# Inverse Simpson diversity or hill number 2
hill2_species <- ggplot(data = plant_traits_MB, 
                        aes(x = Species, y = hill_2)) +
    geom_point(
  aes(color = Species),
  alpha = 0.2,
  size = 1.7,
  position = position_jitter(w = 0.2)
) +
  geom_boxplot(
    aes(color = Species),
    alpha = 0.5,
    outlier.alpha = 0.0) +
  scale_color_manual(name = "Species", values= met.brewer("Isfahan2", n = 4, type = "discrete", direction = -1)) +
  theme_sppdiv_discrete +
  labs(y =  "Alpha diversity (\U1D492 = 2)", x = "", col = "") + 
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test", 
                     label = "p.signif") 
 # facet_wrap(~ Species, labeller = labeller(Species = labels_list)) +
 #guides(color = "none")        #To turn off color legend
  
#hill2_species

# Joined plots
 hill_numbers_species <- ggarrange(hill0_species, hill1_species, hill2_species, 
                                  nrow = 1, ncol = 3,
                                  common.legend = TRUE,
                                  labels = c("A", "B", "C"),
                                  legend = "bottom")

#hill_numbers_species
#ggsave(filename = file.path(path, "figures/hill_numbers_species.png"), plot = hill_numbers_species, dpi=300, units=c("mm"), width=250, height=275)

#################################
#### Beta diversity ####
#### Across elevation ####
#################################

# Observed richness
# Hill number 0
 hill0_beta <-
   ggplot(plant_traits_MB, aes(x = Elevation_cat, y = Observed)) +
   geom_point(
     aes(color = Elevation_cat),
     alpha = 0.2,
     size = 1.7,
     position = position_jitter(w = 0.2)
   ) +
   geom_boxplot(aes(color = Elevation_cat),
                alpha = 0.5,
                outlier.alpha = 0.0) +
   scale_color_manual(name = "Elevation",
     labels = c('Low', 'Mid', 'High'),
     values = met.brewer(
       "Cassatt1",
       n = 3,
       type = "discrete",
       direction = -1
     )
   ) +
   scale_shape_manual(values = c(15:18)) +
   theme_elevdiv_discrete +
   theme(legend.text = element_text(size = 14)) +
   scale_y_continuous(expand = c(0.05, 0.05)) +
   stat_compare_means(comparisons = elevation_comparisons,
                      method = "wilcox.test",
                      label = "p.signif") +
   labs(y = "Beta diversity (\U1D492 = 0)", x = "") +
   guides(fill = guide_legend(title = "Elevation"),
          shape = guide_legend(title = "Species"))

#hill0_beta

#ggsave(file.path(path, "figures/shannon_elevation.png") , plot = bbe, dpi=300, units=c("mm"), width=200, height=150)

# Hill number 1
# Shannon diversity or hill number 1
 hill1_beta <-
   ggplot(plant_traits_MB, aes(x = Elevation_cat, y = hill_1)) +
     geom_point(
     aes(color = Elevation_cat),
     alpha = 0.2,
     size = 1.7,
     position = position_jitter(w = 0.2)
   ) +
   geom_boxplot(aes(color = Elevation_cat),
                alpha = 0.5,
                outlier.alpha = 0.0) +
   scale_color_manual(name = "Elevation",
     labels = c('Low', 'Mid', 'High'),
     values = met.brewer(
       "Cassatt1",
       n = 3,
       type = "discrete",
       direction = -1
     )
   ) +
   scale_shape_manual(values = c(15:18)) +
   theme_elevdiv_discrete +
   theme(legend.text = element_text(size = 14)) +
   scale_y_continuous(expand = c(0.05, 0.05)) +
   stat_compare_means(comparisons = elevation_comparisons,
                      method = "wilcox.test",
                      label = "p.signif") +
   labs(y = "Beta diversity (\U1D492 = 1)", x = "") +
   guides(fill = guide_legend(title = "Elevation"),
          shape = guide_legend(title = "Species"))

#hill1_beta


# Hill number 2
# Inverse Simpson diversity or hill number 2

hill2_beta <- ggplot(plant_traits_MB, aes(x = Elevation_cat, y = hill_2)) +
  
     geom_point(
     aes(color = Elevation_cat),
     alpha = 0.2,
     size = 1.7,
     position = position_jitter(w = 0.2)
   ) +
   geom_boxplot(aes(color = Elevation_cat),
                alpha = 0.5,
                outlier.alpha = 0.0) +
   scale_color_manual(name = "Elevation",
     labels = c('Low', 'Mid', 'High'),
     values = met.brewer(
       "Cassatt1",
       n = 3,
       type = "discrete",
       direction = -1
     )
   ) +
  scale_shape_manual(values=c(15:18)) +
  theme_elevdiv_discrete +
  theme(legend.text = element_text(size = 14))+
  scale_y_continuous(expand = c(0.05, 0.05)) +
  stat_compare_means(comparisons = elevation_comparisons,
                     method = "wilcox.test", 
                     label = "p.signif") +
  labs(y = "Beta diversity (\U1D492 = 2)", x = "") +
  guides(fill = guide_legend(title = "Elevation"), shape = guide_legend(title = "Species"))

#hill2_beta


# Joining plots
hill_beta <- ggarrange(hill0_beta, hill1_beta, hill2_beta,
          nrow = 1, ncol = 3,
          common.legend = TRUE,
          legend = "bottom",
          labels = c("D", "E", "F"))
# hill_beta

#### Alpha and Beta diversity ####
# Joined plots
# A joined plot of joined plots


alpha_beta_diversity <- ggarrange(hill_numbers_species, "", hill_beta,
                          nrow = 3, ncol = 1,
                          common.legend = FALSE,
                          widths = c(1, 0.1, 1),
                          heights = c(1, 0.1, 1),
                          legend = "bottom")

alpha_beta_diversity

ggsave(filename = file.path(path, "figures/supp_fig4.jpg"), plot = alpha_beta_diversity, dpi=300, units=c("mm"), width=300, height=275)

# Alpha and beta diversity plot for defense talk
# defense_alpha_beta <- ggarrange(hill1_diversity, hill1_species,
#                           nrow = 2, ncol = 1,
#                           common.legend = TRUE,
#                           legend = "bottom",
#                           labels = c("A", "B"))
# defense_alpha_beta
# 
# ggsave(filename = file.path(path, "figures/defense_alpha_beta.svg"), plot = defense_alpha_beta, dpi=300, units=c("mm"), width=300, height=300)

```


## Figure S5

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM25
#| fig-cap: "Figure S5. Comparison of log-transformed leaf mass per area (LMA) by species and elevation category. A) Compares logLMA means by species and elevation category: LOW (<1220 m), MID (1221 - 1828 m), HIGH (> 1829 m). B) Change in logLMA per species as elevation increases. The black trend line represents the fitted model for all data points, see main text for *R^2^~adj~* and *p* values. Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < 0.01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."


########## Labels ###########
#Legend title label
legend_title <- "Species"
# Species labels
labels_list <- c("M. laciniatus", "M. nasutus", "M. guttatus")
names(labels_list) <- c("L", "N", "G")
#############################

#ggpubr package used for comparisons tests
lma_compare_means <- compare_means(logLMA ~ Species,
         data = leaf_traits_noB,
         group.by = "Elevation_cat",
         method = "wilcox.test",
         p.adjust.method = "BH",
  paired = FALSE,
  alternative = "two.sided", 
  mu = 0, 
  var.equal = FALSE,
  conf.level = 0.95)

#lma_compare_means
#Report with `compare_means` because it is easier to plot with `ggpubr` package.

#Error bars for plots 
errbar_lims <- group_by(leaf_traits_noB, Site) |> 
              summarize(mean = mean(logLMA), se = sd(logLMA)/sqrt(n()), 
                        upper = mean+(2*se), lower = mean-(2*se))
#Comparisons
species_comparisons <- list(c("M. nasutus", "M. laciniatus"), c("M. nasutus", "M. guttatus"), c("M. laciniatus", "M. guttatus"))

#By Site discrete comparison
lma <- ggplot(leaf_traits_noB, aes(y = logLMA, x = Species)) +
 geom_point(
    aes(color = Species),
    alpha = 0.2,
    size = 1.7,
    position = position_jitter(w = 0.2)
  ) +
  geom_boxplot(aes(color = Species), alpha = 0.5, outlier.alpha = 0.0) +
  scale_color_manual(labels =c('M. laciniatus', 'M. nasutus', 'M. guttatus'),
                    values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  #scale_color_manual(aes(color = "black")) +
  theme_lfspp_discrete +
  labs(y = "log10[Leaf Mass per Area (LMA) (mg/mm)]", x = "") +
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test",
                     #label.y = c(0.4, 0.6, 1),
                     label = "p.signif") +
  facet_wrap(~Elevation_cat)

#lma
#ggsave(filename = file.path(path, "figures/lma_violin.png"), plot = lma, dpi=300, units=c("mm"), width=250, height=150)

# By Elevation continuous comparison
# Linear model summary
# summary(lm(logLMA ~ Elevation_m, data = leaf_traits_noB |>
#              filter(!logLMA == "-Inf"), na.action = na.omit)) # Remove -Inf values
elev_lma_p <- format.p(cor.test(leaf_traits_noB$logLMA, leaf_traits_noB$Elevation_m, na.action = na.omit)$p.value) # This is just to make the p-value nicer. 

# logLMA by Elevation
lma_cont <- ggplot(leaf_traits_noB |> # Removing infinite values
                     filter(!logLMA == "-Inf"), aes(y = logLMA, x = Elevation_m)) +
  geom_jitter(aes(color = Species), alpha = 0.6,
             size = 2.5) +
  geom_smooth(aes(color = Species), method = lm, se = FALSE, na.rm = F) +
  geom_smooth(method = lm, se = TRUE, level = 0.95, na.rm = F, color = "black") +
  stat_poly_eq(
    formula = y ~ x,
    aes(
      color = Species,
      label = paste(#..eq.label..,
        ..adj.rr.label..,
        #"~italic(p) ==",
        ..p.value.label..,
        #signif(..p.value.., digits = 4),
        sep = "*`,`~~")
    ),
    parse = TRUE,
    label.x.npc = "left",
    vstep = 0.08
  )+
  scale_color_manual(labels =c('M. laciniatus', 'M. nasutus', 'M. guttatus'),
                    values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  theme_lfspp_continuous +
  labs(y = "log10[Leaf Mass per Area (LMA)]", x = "Elevation (m.a.s.l.)") 

#lma_cont
#ggsave(filename = file.path(path, "figures/lma_elev.png"), plot = lma_cont, dpi=300, units=c("mm"), width=300, height=150)

# Joining plots
lma_elevation <- ggarrange(lma_cont, lma, ncol = 1, nrow = 2, 
                           common.legend = TRUE, 
                           legend = "bottom",
                           labels = c("A", "B"),
                           heights = c(1, 1))

lma_elevation

ggsave(filename = file.path(path, "figures/supp_fig5.jpg"), plot = lma_elevation, dpi=300, units=c("mm"), width=250, height=250)
```

## Figure S6

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 5
#| label: suppfigMIM26
#| fig-cap: "Figure S6. Change in log-transformed anthocyanin content index (ACI) by species and elevation category. A) Compares logACI means by species and elevation category: LOW (<1220 m), MID (1221 - 1828 m), HIGH (> 1829 m). The black trend line represents the fitted model for all data points, see main text for *R^2^~adj~* and *p* values. B) Change in logLPS per species as elevation increases.  Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

#ggpubr package used for comparisons tests
aci_compare_means <- compare_means(logACI ~ Species,
         data = leaf_traits_noB,
         group.by = "Elevation_cat",
         method = "wilcox.test",
         p.adjust.method = "BH",
  paired = FALSE,
  alternative = "two.sided", 
  mu = 0, 
  var.equal = FALSE,
  conf.level = 0.95)

#aci_compare_means
#Roughly the same outcome as pairwise_t_test.

#Report with `compare_means` because it is easier to plot with `ggpubr` package.
errbar_lims <- group_by(leaf_traits_noB, Species) |> 
              summarize(mean = mean(logACI), se = sd(logACI)/sqrt(n()), 
                        upper = mean + (2*se), lower = mean - (2*se))

#Comparisons
species_comparisons <- list(c("M. nasutus", "M. laciniatus"), c("M. nasutus", "M. guttatus"), c("M. laciniatus", "M. guttatus"))


aci <- ggplot(leaf_traits_noB, aes(y = logACI, x = Species)) +
   geom_point(
    aes(color = Species),
    alpha = 0.2,
    size = 1.7,
    position = position_jitter(w = 0.2)
  ) +
  geom_boxplot(aes(color = Species), alpha = 0.5, outlier.alpha = 0.0) +
  scale_color_manual(values= met.brewer("Isfahan2", n = 4, type = "discrete", direction = -1)) +
  theme_lfspp_discrete +
  theme(axis.text.y = element_text(size = 10)) +
  labs(y = paste("log[Anthocyanin Content Index (ACI)]"), x = "") +
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test",
                     #label.y = c(0.4, 0.6, 1),
                     label = "p.signif") +
  facet_wrap(~Elevation_cat)

#aci
#ggsave(filename = file.path(path, "figures/aci_violin.png"), plot = aci, dpi=300, units=c("mm"), width=300, height=150)

# By Elevation continuous comparison
# Linear model summary
#summary(lm(logACI~ Elevation_m, data = leaf_traits_noB, na.action = na.omit))
elev_aci_p <- format.p(cor.test(leaf_traits_noB$logACI, leaf_traits_noB$Elevation_m, na.action = na.omit)$p.value) # This is just to make the p-value nicer. 

# logACI by Elevation
aci_cont <- ggplot(leaf_traits_noB, aes(y = logACI, x = Elevation_m)) +
  geom_jitter(aes(color = Species), alpha = 0.6,
             size = 2.5)  +
  geom_smooth(aes(color = Species), method = lm, se = FALSE, na.rm = F) +
  geom_smooth(method = lm, se = TRUE, level = 0.95, na.rm = F, color = "black") +
  stat_poly_eq(formula = y ~ x, 
               aes(color = Species,
                   label = paste(#..eq.label.., 
                                 ..adj.rr.label.., 
                                 #"~italic(p) ==",
                                 ..p.value.label.., 
                                 #signif(..p.value.., digits = 4),
                                 sep = "*`,`~~")), 
               parse = TRUE,
               label.x.npc = "left",
               vstep = 0.08) +
  scale_color_manual(values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  theme_lfspp_continuous +
  theme(#plot.caption = element_text(size = 10),
        axis.text.y = element_text(size = 10)) +
  labs(size = 10,
    y = paste("log[Anthocyanin Content Index (ACI)]"),
    x = "Elevation (m.a.s.l.)"
    # caption = (expression(
    #   paste("ACI = % ", italic("A") [lambda], " at 931 nm / % ", italic("A") [lambda], " at 530 nm")
    # ))
  )+
  scale_y_continuous(expand = c(0.07, 0.07)) 

#aci_cont
#ggsave(filename = file.path(path, "figures/aci_cont.png"), plot = aci_cont, dpi=300, units=c("mm"), width=300, height=150)

# Joining plots
aci_elevation <- ggarrange("", aci_cont, "", aci, 
                           ncol = 1, nrow = 4, 
                           heights = c(0.1, 1.1, 0.1, 1.1),
                           common.legend = TRUE, 
                           legend = "bottom",
                           labels = c("A","", "B")
                           )

aci_elevation
ggsave(filename = file.path(path, "figures/supp_fig6.jpg"), plot = aci_elevation, dpi=300, units=c("mm"), width=250, height=250)
```

## Figure S7

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM27
#| fig-cap: "Figure S7. Change in log-transformed leaf lobe index (LBI) by species and elevation category. A) Compares logLBI means by species and elevation category: LOW (<1220 m), MID (1221 - 1828 m), HIGH (> 1829 m). The black trend line represents the fitted model for all data points, see main text for *R^2^~adj~* and *p* values. B) Change in logLPS per species as elevation increases. Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

#ggpubr package used for comparisons tests
lbi_compare_means <- compare_means(logLBI ~ Species,
         data = leaf_traits_noB,
         group.by = "Elevation_cat",
         method = "wilcox.test",
         p.adjust.method = "BH",
  paired = FALSE,
  alternative = "two.sided", 
  mu = 0, 
  var.equal = FALSE,
  conf.level = 0.95)

#lbi_compare_means


#Roughly the same outcome as pairwise_t_test.

#Report with `compare_means` because it is easier to plot with `ggpubr` package.
errbar_lims <- group_by(leaf_traits_noB, Species) |> 
              summarize(mean = mean(logLBI), se = sd(logLBI)/sqrt(n()), 
                        upper = mean + (2*se), lower = mean - (2*se))

#Comparisons
species_comparisons <- list(c("M. nasutus", "M. laciniatus"), c("M. nasutus", "M. guttatus"), c("M. laciniatus", "M. guttatus"))

lbi <- ggplot(leaf_traits_noB, aes(y = logLBI, x = Species)) +
 geom_point(
    aes(color = Species),
    alpha = 0.2,
    size = 1.7,
    position = position_jitter(w = 0.2)
  ) +
  geom_boxplot(aes(color = Species), alpha = 0.5, outlier.alpha = 0.0) +
  scale_color_manual(labels =c('M. laciniatus', 'M. nasutus', 'M. guttatus'),
                    values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  theme_lfspp_discrete +
  labs(y = "log[Leaf Lobe Index (LBI)", x = "") +
  #scale_x_discrete(labels = c('M. laciniatus', 'M. nasutus', 'M. guttatus')) +
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test", 
                     label = "p.signif") +
   facet_wrap(~Elevation_cat)

#lbi
#ggsave(filename = file.path(path, "figures/lbi_violin.png"), plot = lbi, dpi=300, units=c("mm"), width=300, height=150)

# By Elevation continuous comparison
# Linear model summary
#summary(lm(logLBI~ Elevation_m, data = leaf_traits_noB, na.action = na.omit))
elev_lbi_p <- format.p(cor.test(leaf_traits_noB$logLBI, leaf_traits_noB$Elevation_m, na.action = na.omit)$p.value) # This is just to make the p-value nicer. 

# logLBI by Elevation
lbi_cont <- ggplot(leaf_traits_noB, aes(y = logLBI, x = Elevation_m)) +
  geom_jitter(aes(color = Species), alpha = 0.6,
             size = 2.5)  +
  geom_smooth(aes(color = Species), method = lm, se = FALSE, na.rm = F) +
  geom_smooth(method = lm, se = TRUE, level = 0.95, na.rm = F, color = "black") +
  stat_poly_eq(formula = y ~ x, 
               aes(color = Species,
                   label = paste(#..eq.label.., 
                                 ..adj.rr.label.., 
                                 #"~italic(p) ==",
                                 ..p.value.label.., 
                                 #signif(..p.value.., digits = 4),
                                 sep = "*`,`~~")), 
               rr.digits = 2,
               p.digits = 3 ,
               parse = TRUE,
               label.x.npc = "left",
               vstep = 0.08) +
  scale_color_manual(labels =c('*M. laciniatus*', '*M. nasutus*', '*M. guttatus*'),
                    values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  theme_classic(base_size = 16) +
  theme(legend.title = element_text(face = "bold"),
        panel.border = element_rect(linetype = "blank", fill = NA),
        legend.text = element_markdown(size = 14), #Make legend text italic
        axis.text.x = element_text(size = 14),
        legend.position = "bottom") +
  labs(y = "log[Leaf Lobe Index (LBI)]", x = "Elevation (m.a.s.l.)") +
  scale_y_continuous(expand = c(0.07, 0.07)) 

#lbi_cont
#ggsave(filename = file.path(path, "figures/lbi_elev.png"), plot = lbi_cont, dpi=300, units=c("mm"), width=300, height=150)

# Joining plots
lbi_elevation <- ggarrange(lbi_cont, lbi, ncol = 1, nrow = 2, 
                           common.legend = TRUE, 
                           legend = "bottom",
                           labels = c("A", "B"),
                           heights = c(1, 1))

lbi_elevation

ggsave(filename = file.path(path, "figures/supp_fig7.jpg"), plot = lbi_elevation, dpi=300, units=c("mm"), width=250, height=250)
```

## Figure S8

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM28
#| fig-cap: "Figure S8. Change in log-transformed leaf-punch strength (LPS), a measure of leaf toughness, by species and elevation category. A) Compares logLPS means by species and elevation category: LOW (<1220 m), MID (1221 - 1828 m), HIGH (> 1829 m). The black trend line represents the fitted model for all data points, see main text for *R^2^~adj~* and *p* values. B) Change in logLPS per species as elevation increases. Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

#ggpubr package used for comparisons tests
lps_compare_means <- compare_means(logLPS ~ Species,
         data = leaf_traits_noB,
         group.by = "Elevation_cat",
         method = "wilcox.test",
         p.adjust.method = "BH",
  paired = FALSE,
  alternative = "two.sided", 
  mu = 0, 
  var.equal = FALSE,
  conf.level = 0.95)

#lps_compare_means

#Report with `compare_means` because it is easier to plot with `ggpubr` package.
errbar_lims <- group_by(leaf_traits_noB, Species) |> 
              summarize(mean = mean(logLPS), se = sd(logLPS)/sqrt(n()), 
                        upper = mean + (2*se), lower = mean - (2*se))


lps <- ggplot(leaf_traits_noB, aes(y = logLPS, x = Species)) +
 geom_point(
    aes(color = Species),
    alpha = 0.2,
    size = 1.7,
    position = position_jitter(w = 0.2)
  ) +
  geom_boxplot(aes(color = Species), alpha = 0.5, outlier.alpha = 0.0) +
  scale_color_manual(values = met.brewer(
    "Isfahan2",
    n = 4,
    type = "discrete",
    direction = c(-1)
  )) +
  theme_lfspp_discrete +
  labs(y = expression(paste("Leaf punch strength (LPS)(N mm", "^-1)")), x = "") +
  #scale_x_discrete(labels = c('M. laciniatus', 'M. nasutus', 'M. guttatus')) +
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test",
                     #label.y = c(0.4, 0.6, 1),
                     label = "p.signif") +
  facet_wrap( ~ Elevation_cat)
  
#lps
#ggsave(filename = file.path(path, "figures/lps_violin.png"), plot = lps, dpi=300, units=c("mm"), width=300, height=150)

# By Elevation continuous comparison
# Linear model summary
#summary(lm(logLPS ~ Elevation_m, data = leaf_traits_noB, na.action = na.omit))
elev_lps_p <- format.p(cor.test(leaf_traits_noB$logLPS, leaf_traits_noB$Elevation_m, na.action = na.omit)$p.value) # This is just to make the p-value nicer. 

# logLPS by Elevation
lps_cont <-
  ggplot(leaf_traits_noB, aes(y = logLPS, x = Elevation_m)) +
  geom_jitter(aes(color = Species), alpha = 0.6,
              size = 2.5)  +
  geom_smooth(aes(color = Species),
              method = lm,
              se = FALSE,
              na.rm = F) +
  geom_smooth(
    method = lm,
    se = TRUE,
    level = 0.95,
    na.rm = F,
    color = "black"
  )+
  stat_poly_eq(
    formula = y ~ x,
    aes(
      color = Species,
      label = paste(#..eq.label..,
        ..adj.rr.label..,
        #"~italic(p) ==",
        ..p.value.label..,
        #signif(..p.value.., digits = 4),
        sep = "*`,`~~")
    ),
    rr.digits = 2,
    p.digits = 3 ,
    parse = TRUE,
    label.x.npc = "left",
    vstep = 0.08
  ) +
  scale_color_manual(values = met.brewer(
    "Isfahan2",
    n = 4,
    type = "discrete",
    direction = c(-1)
  )) +
  theme_lfspp_continuous +
  labs(y = expression(paste("Leaf punch strength (LPS) (N mm", "^-1)")), x = "Elevation (m.a.s.l.)") +
  scale_y_continuous(expand = c(0.07, 0.07)) 

#lps_cont
#ggsave(filename = file.path(path, "figures/lps_elev.png"), plot = lps_cont, dpi=300, units=c("mm"), width=300, height=150)

# Joining plots
lps_elevation <- ggarrange("", lps_cont, "", lps, ncol = 1, nrow = 4, 
                           common.legend = TRUE, 
                           legend = "bottom",
                           labels = c( "A","", "B"),
                           heights = c(0.1, 1, 0.1, 1))

lps_elevation

ggsave(filename = file.path(path, "figures/supp_fig8.jpg"), plot = lps_elevation, dpi=300, units=c("mm"), width=250, height=250)
```

## Figure S9

```{r}
#| echo: false
#| eval: true
#| tidy: true
#| fig-width: 6
#| fig-height: 6
#| label: suppfigMIM29
#| fig-cap: "Figure S9. Change in log-transformed leaf thickness (LT) (𝜇m) by species and elevation category. A) Compares logLT means by species and elevation category: LOW (<1220 m), MID (1221 - 1828 m), HIGH (> 1829 m). The black trend line represents the fitted model for all data points, see main text for *R^2^~adj~* and *p* values. B) Change in logLT per species as elevation increases. Significance levels are represented by *ns* (not significant) and asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

#ggpubr package used for comparisons tests
lt_compare_means <- compare_means(logLT ~ Species,
         data = leaf_traits_noB,
         group.by = "Elevation_cat",
         method = "wilcox.test",
         p.adjust.method = "BH",
  paired = FALSE,
  alternative = "two.sided", 
  mu = 0, 
  var.equal = FALSE,
  conf.level = 0.95)

#lt_compare_means
#Roughly the same outcome as pairwise_t_test.

#Report with `compare_means` because it is easier to plot with `ggpubr` package.
errbar_lims <- group_by(leaf_traits_noB, Species) |> 
              summarize(mean = mean(logLT), se = sd(logLT)/sqrt(n()), 
                        upper = mean + (2*se), lower = mean - (2*se))

#Comparisons
species_comparisons <- list(c("M. nasutus", "M. laciniatus"), c("M. nasutus", "M. guttatus"), c("M. laciniatus", "M. guttatus"))

lt <- ggplot(leaf_traits_noB, aes(y= logLT, x = Species)) +
 geom_point(
    aes(color = Species),
    alpha = 0.2,
    size = 1.7,
    position = position_jitter(w = 0.2)
  ) +
  geom_boxplot(aes(color = Species), alpha = 0.5, outlier.alpha = 0.0) +
  scale_color_manual(values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  theme_lfspp_discrete +
  labs(y = expression(paste("log[Leaf thickness (LT) (",mu,"m)]")), x = "") +
  #scale_x_discrete(labels = c('M. laciniatus', 'M. nasutus', 'M. guttatus')) +
  stat_compare_means(comparisons = species_comparisons,
                     method = "wilcox.test",
                     #label.y = c(0.4, 0.6, 1),
                     label = "p.signif", hide.ns = TRUE) +
   facet_wrap(~Elevation_cat)

#lt
#ggsave(filename = file.path(path, "figures/lt_violin.png"), plot = lt, dpi=300, units=c("mm"), width=300, height=150)

# By Elevation continuous comparison
# Linear model summary
#summary(lm(logLT ~ Elevation_m, data = leaf_traits_noB, na.action = na.omit))
#elev_lt_p <- format.p(cor.test(leaf_traits_noB$logLT, leaf_traits_noB$Elevation_m, na.action = na.omit)$p.value) # This is just to make the p-value nicer. 

# logLT by Elevation
lt_cont <- ggplot(leaf_traits_noB, aes(y = logLT, x = Elevation_m)) +
  geom_jitter(aes(color = Species), alpha = 0.6,
             size = 2.5)  +
  geom_smooth(aes(color = Species), method = lm, se = FALSE, na.rm = F) +
  geom_smooth(method = lm, se = TRUE, level = 0.95, na.rm = F, color = "black") +
  stat_poly_eq(formula = y ~ x, 
               aes(color = Species,
                   label = paste(#..eq.label.., 
                                 ..adj.rr.label.., 
                                 #"~italic(p) ==",
                                 ..p.value.label.., 
                                 #signif(..p.value.., digits = 4),
                                 sep = "*`,`~~")), 
               rr.digits = 2,
               p.digits = 3 ,
               
               parse = TRUE,
               label.x.npc = "left",
               vstep = 0.08) +
  scale_color_manual(values = met.brewer("Isfahan2", n = 4, type = "discrete", direction = c(-1))) +
  theme_lfspp_continuous +
  labs(y = expression(paste("log[Leaf thickness (LT) (",mu,"m)]")), x = "Elevation (m.a.s.l.)") +
  scale_y_continuous(expand = c(0.07, 0.07)) 

#lt_cont
#ggsave(filename = file.path(path, "figures/lps_elev.png"), plot = lps_cont, dpi=300, units=c("mm"), width=300, height=150)

# Joining plots
lt_elevation <- ggarrange(lt_cont, lt, ncol = 1, nrow = 2, 
                           common.legend = TRUE, 
                           legend = "bottom",
                           labels = c("A", "B"),
                           heights = c(1, 1))
lt_elevation
ggsave(filename = file.path(path, "figures/supp_fig9.jpg"), plot = lt_elevation, dpi=300, units=c("mm"), width=250, height=250)
```


## Table S1
```{r, supp_table1}
# Sites and Species sampled
leaf_traits_noB |>
  group_by(Site,Species, Longitude, Latitude) |>
  summarise(n = as.integer(n()/3))|> # Dividing by 3 due to replicate measurements of each sample
  pivot_wider(names_from = Species, values_from = n)|>
  
  # huxtable::huxtable(number_format(fmt_pretty(.))) 
  # #huxtable::set_bold(1) |>
  # huxtable::set_italic(1, 4:6, TRUE) |>
  # set_all_padding(10) |>
  # set_outer_padding(0) |>
  # set_width(1.1) |>
  # theme_article() |>
  # huxtable::set_caption("Table 1: Number of individuals collected per population/site and species.") |>
  # huxtable::as_flextable(column_to_header = TRUE)
  knitr::kable(caption = "Table S1. Number of individuals collected per population/site and species.")
```


## Table S2
```{r, supp_table2}
phyloseq_ntaxa_by_tax(
  ps_clean_3,
  TaxRank ="Phylum",
  relative = F,
  add_meta_data = T
) |>
  as.data.frame() |>
  mutate(sum = sum(N.OTU),
          Phylum = str_remove(Phylum, "^.*?(_|-).")) |>
  group_by(Species, Phylum) |>
  summarise(Occurrence = n()) |>
  knitr::kable(caption = "Table S2. Prevalence of phyla in samples per species.")
```

## Table S3
```{r, supp_table3}
#| echo: false
#| eval: true
#| tidy: true
#| tbl-cap: "Table S3. Permutational ANOVA for the best-fit dbRDA model with hybrid cross-types. Significance levels are represented by asterisks [*p* < .05 (\\*), *p* < .01 (\\*\\*), *p* < .001 (\\*\\*\\*), and *p* < .0001 (\\*\\*\\*\\*)]."

#dbRDA table  (results from terms are the same as adonis2 but without R2)
load("clean_data/statistics/dbRDA_models.rda")

anova_table3 <- tidy(anova_m5hell_term) |>
  as.data.frame() |>
  mutate("p-value" = str_replace_all(p.value, "0.001", "***")) |>
  rename(
    Term = term,
    "Sum. of Sqs." = SumOfSqs,
    "F-value" = statistic,
    "Df" = df
  ) |> 
  select(Term, Df, "Sum. of Sqs.", "F-value", "p-value")

# Table

table3 <- flextable(anova_table3) |>
  flextable::set_table_properties(width = .5) |>
  flextable::set_table_properties(layout = "autofit") |>
  flextable::style(i = c(1:4),
                   j = 1,
                   pr_t = fp_text(bold = FALSE)) |>
  flextable::bold(
    i = c(1),
    j = c(1:5),
    bold= TRUE,
    part = "header"
  ) |>
  flextable::italic(
    i = c(4),
    j = 1,
    italic = TRUE,
    part = "body"
  ) |>
  flextable::font(
    i = NULL,
    j = NULL,
    fontname = "Tex Gyre Termes",
    part = "all",
    cs.family = fontname,
    hansi.family = fontname,
    eastasia.family = fontname
  ) |>
  flextable::set_caption("Table S3. Permutational ANOVA for the best-fit dbRDA model with hybrid cross-types. Significance levels are represented by asterisks")
  
table3
```

## Supplementary methods

### Illumina TruSeq adapters

The modified primers for the first PCR (adapter ligation and ITS1 amplification) were as follows: 5’-CAC TCT TTC CCT ACA CGA CGC TCT TCC GAT CTC TTG GTC ATT TAG AGG AAG TAA-3’ (forward) and 5’-GTG ACT GGA GTT CAG ACG TGT GCT CTT CCG ATC TGC TGC GTT CTT CAT CGA TGC-3’ (reverse).
See [@gardes1993] and [@white1990_inbook] for more details.

Supplementary Literature Cited