ms.Rmd

---
title: Growth Dynamics of the Independent Gametophytes of *Pleurosoriopsis makinoi*
  (Polypodiaceae)
output:
  pdf_document:
    latex_engine: xelatex
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, cache = TRUE)

# Set order of microclimate variables for printing.
order_vars <- c("par_total", "rh_min", "temp_mean")

# Load objects from drake workflow into memory
tar_load(c(
  selected_vars,
  combined_morph,
  daily_microclimate,
  t_test_results,
  paired_microclimate,
  nested_paired_microclimate,
  combined_monthly_morph,
  linear_model_results,
  gameto_size)
)
```

Atsushi Ebihara^1^\*, Joel H. Nitta^1^\*\*, Yurika Matsumoto^2^, Yuri Fukazawa^2^, Marie Kurihara^2^, Hitomi Yokote^2^, Kaoru Sakuma^2^, Otowa Azakami^2^, Yumiko Hirayama1 and Ryoko Imaichi^2^

^1^Department of Botany, National Museum of Nature and Science, 4–1–1 Amakubo, Tsukuba, Ibaraki 305–0005, Japan

^2^Department of Chemical and Biological Sciences, Japan Women’s University, 2-8-1 Mejirodai, Bunkyo-ku, Tokyo 112-8681, Japan

\* E-mail: ebihara@kahaku.go.jp

\*\* Current affiliation: Department of Botany, National Museum of Natural History, Smithsonian Institute, Washington, D.C., 20013, U.S.A.

(Received 12 February 2019; accepted 27 March 2019)

**Abstract** Despite their importance to the fern life cycle, little is known about the ecology of fern gametophytes due to their cryptic morphology and small size. This study documents growth dynamics over five years in a population of *Pleurosoriopsis makinoi* gametophytes growing independently of sporophytes in Okutama-machi, Tokyo, Japan. Cover was measured in four 10 × 10 cm quadrats, and the number and size of gemmae born on gametophytes measured ca. every two months from March 2009 to March 2014. The population survived for the duration of the study, and no seasonality was observed in the cover measurements. Gemmae gradually grow throughout the spring and summer, and mostly detach in winter. Based on our preliminary comparison between the microenvironment of the independent gametophyte site and a site supporting sporophytes of the same species, a possible cause of suppression of sporophyte formation is low light conditions in winter due to planting of evergreen conifers.

**Key words**: fern, gemma, microenvironment, sporophyte.

Ferns possess an unusual life cycle consisting of two morphologically distinct, free-living stages: the diploid sporophyte and the haploid gametophyte. Unlike seed plants, in which the gametophyte is nutritionally dependent upon the sporophyte, in ferns gametophytes may grow separately from sporophytes. The phenomenon of ‘independent gametophytes’ includes fern gametophyte populations in various degrees of independence from sporophytes (Pinson *et al.*, 2017). Independent gametophyte species, lacking any known counterpart sporophyte, are rather extreme cases (Farrar, 1967, 1978, 1992). Conversely, facultatively independent species, which may or may not grow separately from sporophytes depending on environmental conditions, have been observed in many more species (Rumsey and Sheffield, 1996). Almost without exception, species capable of growing as independent gametophytes have noncordiform morphology with indeterminate growth (Ebihara *et al.*, 2013). Facultatively independent fern gametophytes represent an understudied but useful group for understanding the ecological correlates of the transition from haploid to diploid phase in land plants.

We recently conducted a survey of the fern gametophyte flora of Japan (Ebihara *et al.*, 2013). During the survey, we located a population of *Pleurosoriopsis makinoi* (Maxim. ex Makino) Fomin (Polypodiaceae) gametophytes in Unazawa Valley, Okutama-machi, Tokyo, Japan growing without any sporophytes nearby as far as we searched --- its sporophyte population has evidently been observed in the valley, but is very small and rare (S. Fujimoto, personal communication). We therefore considered it an ‘independent gametophyte’ in the broad sense. *Pleurosoriopsis makinoi* is a small-leaved fern species often growing with long-creeping rhizomes among bryophyte mats on rocks distributed in East Asia (China, Japan, Korea and Russian Far East) (Zhang *et al.*, 2013). It has ribbon-like gametophytes, and gemmae and sexual organs have been observed both in the wild and in cultivation from sown spores (Masuyama, 1975).

We conducted a continuous observational study on the wild independent gametophytes of *P. makinoi* in Unazawa Valley to answer the following questions:
1.	How stable is the size of independent gametophyte populations over long periods of time?
2.	Do independent gametophytes produce sexual organs (archegonia and antheridia) or possess the ability to produce sporophytes?
3.	Do the gemmae of independent gametophytes contribute to sustainment/expansion of populations?
4.	How do environmental factors affect gametophyte growth?

## Materials and Methods

### Location and habitat

The independent gametophyte site is located in Unazawa Valley, Okutama-machi, Tokyo, Japan, at ca. 315 m elevation. It grows on a small north-facing bank along a river (Fig. 1). The population is ca. 1 × 5 m in size.

### Measurement of cover

Four 10 cm × 10 cm quadrats (q1, q2, q3 and q4) were designated on the rockface (Fig. 1) where the independent gametophyte population grows to measure cover (growth area). As the ribbon-like gametophytes of *Hymenophyllum barbatum* (Bosch) Baker (Hymenophyllaceae) also grow on the rockface, we carefully situated the quadrats to include only *P. makinoi*. We were able to distinguish between gametophytes of the two species based on the results of Ebihara *et al.* (2013) and morphological observation --- unlike *P. makinoi*, gametophytes of *H. barbatum* do not produce gemmae. Each quadrat was photographed using a digital camera. The image files were printed out, and a 1 cm^2^ grid was overlayed on the photo with tracing paper. Cover scores for each 1 cm^2^ cell in the grid were estimated in four discrete levels ranging from zero to 1, i.e. 1.0 cm^2^ (100% ≥ cover > 50%), 0.5 cm^2^ (50% ≥ cover > 25%), 0.25 cm^2^ (25% ≥ cover > 0%), 0 (completely absent). Cover was calculated as the sum of the scores for each cell. The measurement was performed ca. every two months by placing the quadrats in exactly the same positions, and was repeated 37 times from March 2009 to March 2014.

### Gametophyte collection, observation and measurements

Living gametophytes were collected from the peripheries of the quadrats, and were brought back to the laboratory in small zipped plastic bags with wet papers. Collections were made 38 times from March 2009 to March 2014. The gametophytes were washed using brushes under a stereomicroscope (Leica M420). The number of gemmae was counted and length of each gemma (Figs. 2a, b) was measured under a light microscope (Olympus BX51) in ten individuals per collection, and mean of the number and the length of gemmae calculated.

We conducted a survey for reproductive structures in 2015. We collected 90 individuals, measured their length and width, and checked for the presence of antheridia or archegonia. During the 2015 survey, no reproductive structures were observed, so we also inspected images of all vouchers of *P. makinoi* collected at this site for archegonia or antheridia.

In this study, an individual gametophyte is defined as a single distinct thallus.

### Measurement of microenvironment

HOBO Micro Station Data Loggers H21-002 (Onset) with sensors S-TMB-M002 (for temperature), S-THB-M002 (for temperature and humidity) and S-LIA-M003 for photosynthetic light (PAR) were used to measure microenvironment. In Unazawa Valley, the logger was installed from September 22, 2012 to July 22, 2014. To compare the environment of the independent gametophyte habitat with that of conspecific sporophytes, another set of the equipment was installed at a site where sporophytes were growing. The sporophyte site was at Uratakaomachi, Hachioji-shi, Tokyo, at 220 m elevation (from August 20, 2013 to November 22, 2014). Temperature and humidity were recorded every 30 minutes. PAR was sampled every 4 minutes and the average value was recorded every 30 minutes.

The data recorded by the loggers were collected by connecting the logger to a notebook computer (Windows 7) and downloading using BoxCar Pro 4.3 and HOBOWare Pro every two to four months. The collected data were summarized by site as the following values: 1) daily mean temperature, 2) daily minimum humidity, 3) daily light integral (DLI). Daily light integral was calculated by summing the total PAR measurments over each day.

### Statistical analysis

Comparisons of microclimate between the sporophyte site and gametophyte site were made using a paired *t*-test. To test for the influence of microclimate on total cover (sum of all four quadrats), mean gemmae length, and mean gemmae number, linear models were fit with each of these as the dependent variable and the microclimatic variables as the independent variable. All statistical analysis was done using R 3.5.1 (R Core Team 2018). Plots were made with the “ggplot2” package (Wickham 2016). Raw data and code to replicate all analyses and figures are available at https://github.com/joelnitta/pleurosoriopsis.

## Results

### Measurement of cover

The cover of Pleurosoriopsis gametophytes measured in four quadrats over five years is shown in Fig. 3a. Cover irregularly fluctuated lower than 40% of original size, but never reached zero. Cover in three of the four quadrats (q1, q2, and q4) decreased during the last two years of the survey (2012–2014).

### Measurement of number and size of gemmae

The mean number and length of gemmae per individual are summarized in Figs. 3b and 3c. The highest number of gemmae tend to be observed in late summer to autumn and the lowest number tend to be observed in spring (Fig. 3b). Longer gemmae are frequently observed during November to January (Fig. 3c).

### Measurement of microenvironment

The measured temperature, humidity and DLI are summarized in Fig 4. At the gametophyte site, data collection was interrupted by technical malfunction during the following periods: March 8–22, 2014, May 4–6, 2014, May 11–July 13, 2014.

The sporophyte site and gametophyte site differed significantly in nearly all measured microclimatic variables (Fig. 5, Table 1). The gametophyte site had lower light levels, especially in fall (0.21 ± 0.06 mmol per day, vs. 0.33 ± 0.12 at the sporophyte site; *p* < 0.001, paired *t*-test) and winter (0.36 ± 0.27 mmol per day, vs. 1.51 ± 0.74; *p* < 0.001), and tended to be drier, particularly in spring (50 ± 22% min. rel. humidity, vs. 69 ± 21%; *p* < 0.001) and winter (78 ± 18%, vs. 91 ± 10%; *p* < 0.001). Differences in temperature were significant but small (all estimated differences < 1°C; Table 1).

### Relationship between growth and climate

No significant relationship was observed between gametophyte growth and microclimate except for length of gemmae, which was negatively related to temperature (linear model; *r*^2^ = 0.51, *p* = 0.013) (Fig. 6, Table 2).

### Reproductive structures

None of the gametophytes surveyed in 2015 possesed any sexual organs (*n* = 90). When we searched all other images of specimens collected from the Okutama site, archegonia (Figs. 2c, d) were observed in 32 individuals (Fig. 7). Antheridia (Figs. 2e, f) were much less common than archegonia --- we only found them in four individuals. Females tended to be larger than asexuals (it was not possible to conduct statistical analysis of size differences between sexes due to the low number of individuals with antheridia and inconsistency in sampling season and effort). No bisexual individuals were observed. We did not observe any juvenile sporophytes on the gametophytes.

## Discussion

The gametophyte patches in each of the four quadrats survived through all five years of our study. Even though we observed change in cover over time in each quadrat, there seems to be no seasonality in this fluctuation (Fig. 3a). In contrast, clear seasonal patterns were observed in gemmae formation and growth (Figs. 3b, c). Gemmae gradually grow and mostly detach in winter, the coldest and driest season. Considering the fact that no obvious new patches of gametophytes appeared during our study, the function of gemmae as asexual propagules was not demonstrated by our results.  

Although our measurement data of microenvironment are too short-termed to discuss the factors causing the independent gametophyte phenomenon, an interesting trend was observed in the DLI values --- they are higher in the sporophyte site than the gametophyte site, particularly in the fall and winter (Fig. 5, Table 1). The difference is undoubtedly caused by the different tree species covering the *P. makinoi* habitats --- summer-green broad-leaved trees at the sporophyte site (Figs. 8a, b) versus evergreen conifer [Chamaecyparis obtusa (Siebold et Zucc.) Endl.] plantation in Unazawa Valley, Okutama, the independent gametophyte site (Figs. 8c, d). It is highly likely that sporophytes of *P. makinoi*, which have winter-green life history, are probably adapted to the summer-green forest floor environment and are not able to grow and survive on the evergreen forest floor due to a shortage of light. The stable presence of the gametophyte patches in Okutama suggests that the environment is suitable for gametophytes, but not sporophytes. Since we only observed a very small number of individuals with gametangia and no juvenile sporophytes, it seems likely that suppression of sporophyte production occurs prior to fertilization, or shortly thereafter. Archegoniate (female) individuals tended to be larger than asexual individuals (Fig. 7), suggesting that gametophytes may need to reach a size threshold before they are able to sexually reproduce. The low winter light conditions at Okutama may prevent gametophytes from reaching this size consistently and thereby inhibit fertilization.

“Inter-generational ecological niche separation”, or different ecological preferences between sporophytes and gametophytes of the same species, has been reported in several ferns (Sato and Sakai, 1981, Rumsey and Sheffield 1996, Pinson *et al.*, 2017). In most cases, the gametophyte tends to tolerate a wider range of environments than the conspecific sporophyte (Watkins *et al.*, 2007). This seems also to be the case for *P. makinoi* with regards to light condition. A possible scenario for formation of the independent gametophyte population in Okutama is as follows: 1) natural vegetation of the area was summer-green forest preferred by both sporophytes and gametophytes of *P. makinoi*; 2) after evergreen conifer plantations were established in the mid-20th century, sporophyte growth in the plantation areas was hampered by light shortage in winter and eventually ceased; 3) gametophytes still survive by indeterminate growth and vegetative reproduction by gemmae.

We could not obtain any evidence for past presence of sporophytes at the site, but we found a similar case in Kyoto City, Kyoto Prefecture, Japan. At that site as well, evergreen conifer plantations dominated the forests starting from almost the same era. At least six herbarium specimens of *Pleurosoriopsis* sporophytes collected in the mountains surrounding Kyoto from 1928 to 1958 deposited in KYO, TNS and TOFO imply that it was not a very rare species at that time. Nowadays it is quite difficult to find *Pleurosoriopsis* sporophytes in the same area (Y. Kazumi, personal communication), but independent gametophytes are occasionally found (A. Ebihara, unpublished data). Such gametophytes whose sporophyte formation has been suppressed by recent environmental change during the past 100 years may resume their sexual life cycles when or if the environment returns to its prior condition. Further studies are necessary to determine the conditions needed for sporophyte formation, in particular frequencies of sexual organ formation.

#### Acknowledgements

This study was supported by a Grant-in-Aid for Scientific Research (KAKENHI 16370046 to Ryoko Imaichi) from the Japan Society of the Promotion of Science.

## References

Ebihara, A., Yamaoka, A., Mizukami, N., Sakoda, A., Nitta, J. H. and Imaichi, R. 2013. A survey of the fern gametophyte flora of Japan: frequent independent occurrences of noncordiform gametophytes. Amer. J. Bot. 100: 735-743.

Farrar, D. R. 1967. Gametophytes of four tropical fern genera reproducing independently of their sporophytes in the southern Appalachians. Science 155: 1266–1267.

Farrar, D. R. 1978. Problems in the identity and origin of the Appalachian *Vittaria* gametophyte, a sporophyteless fern of the eastern United States. American Journal of Botany 65: 1–12.

Farrar, D. R. 1992. *Trichomanes intricatum*: the independent *Trichomanes* gametophyte in the eastern United States. American Fern Journal 82: 68–74.

Masuyama, S. 1975. Gametophyte of *Pleurosoriopsis makinoi* (Maxim.) Fomin. Journal of Japanese Botany 50: 105-114.

Pinson, J. B., Chambers, S. M., Nitta, J. H., Kuo, L.-Y. and Sessa, E. B. 2017. The separation of generations: biology and biogeography of long-lived sporophyteless fern gametophytes. International Journal of Plant Sciences 178: 1-18.

R Core Team. 2018. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Rumsey, F. J. and Sheffield, E. 1996. Inter-generational ecological niche separation and the 'independent gametophyte' phenomenon. In: Camus, J. M., Gibby, M. and Johns, R. J. (eds.), Pteridology in Perspective, pp. 563-570. Royal Botanic Gardens, Kew, London.

Sato, T. and Sakai, A. 1981. Cold tolerance of gametophytes and sporophytes of some cool temperate ferns native to Hokkaido. Canadian Journal of Botany 59: 604-608.

Watkins Jr, J. E., Mack, M. K., and Mulkey, S. S. 2007. Gametophyte ecology and demography of epiphytic and terrestrial tropical ferns. American Journal of Botany 94: 701-708.

Wickham, H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org/.

Zhang, X. C., Lu, S. G., Lin, Y. X., Qi, X. P., Moore, S., Xing, F. W, Wang, F. G., Hovenkamp, P. H., Gilbert, M. G., Nooteboom, H. P., Parris, B. S., Haufler, C., Kato, M. and Smith, A. R. 2013. Polypodiaceae. In: Wu, Z. Y., Raven, P. H. and Hong, D. Y. (eds.), Flora of China, Vol. 2-3 (Pteridophytes), pp. 758-850. Science Press, Beijing, Missouri Botanical Garden Press, St. Louis.

```{r, results='asis', echo = FALSE, eval=(knitr::opts_knit$get('rmarkdown.pandoc.to') == 'latex')}
# Add pagebreak if rendering to pdf
cat('\\pagebreak')
```

## Figures

```{r fig-1, fig.width=7, fig.height=7, fig.cap = "The quadrats for measurement of gametophyte cover of *Pleurosoriopsis makinoi* in Unazawa Valley, Okutama-machi, Tokyo, Japan. a. A full view of the site with four quadrat positions (q1, q2, q3 and q4). b. A 10 cm × 10 cm quadrat frame."}
knitr::include_graphics(here("images/Fig1_quadrats.pdf"))

```

```{r fig-2, fig.width=7, fig.height=7, fig.cap = "Gametophytes of *Pleurosoriopsis makinoi* collected in Unazawa Valley, Okutama-machi, Tokyo, Japan. a. A typical ribbon-like gametophyte with gemmae on the apical part (arrowhead), collected in August. b. Gemmae (arrowheads) on a gametophyte collected in November. c. A gametophyte with archegonia (arrowheads), collected in October. d. Archegonia. e. A gametophyte with antheridia (arrowhead) collected in December. f. An antheridium. Scale bars: 2 mm for a, 200 μm for b and d, 1 mm for c and e, 100 μm for f."}
knitr::include_graphics(here("images/Fig2_gemma_sexualorgans.pdf"))
```

```{r fig-3, fig.width=8, fig.height=7, fig.cap="Change in growth of *Pleurosoriopsis makinoi* gametophytes observed from May 2009 to May 2014. a. Cover (cm^2^) of gametophytes in four 10 × 10 cm quadrats. b. Mean number of gemmae per individual. c. Mean gemma length (μm). For b and c, *n* = 10 individuals per observation; grey bars show SD. Shaded bars show alternating years."}
#' Set x-limits manually since the automatically set limits will vary
#' when we subset the data.
start_date <- combined_morph %>%
  pull(date) %>%
  min

end_date <- combined_morph %>%
  pull(date) %>%
  max

#' Make tibble of years for shading. This dataframe will be passed to `geom_rect()` when
#' plotting.
shading_dates <-
  tibble(date = seq.Date(from = start_date, to = end_date, by = "day")) %>%
  mutate(year = year(date)) %>%
  # shade by alternate years
  filter(year %in% c(2009, 2011, 2013)) %>%
  group_by(year) %>%
  summarize(
    year_start = min(date),
    year_end = max(date)
  )

#' Make list of subplots. Build each one at a time, then add common
#' formatting.
subplots <- list()

subplots[["cover"]] <-
  combined_morph %>%
  ungroup %>%
  select(starts_with("q"), date) %>%
  filter(complete.cases(.)) %>%
  gather(plot, area, -date) %>%
  ggplot(aes(x = date, y = area, linetype = plot)) %>%
  shade_years(shading_dates) +
  geom_line() +
  scale_linetype_manual(
    name="Quadrat",
    breaks=c("q_1", "q_2", "q_3", "q_4"),
    labels=c("1", "2", "3", "4"),
    values=c("solid", "dashed", "dotted", "dotdash")
  ) +
  labs(
    x = "",
    y = expression("Cover ("~cm^2~")"),
    subtitle = "a"
  ) +
  theme(legend.position = "bottom")

legend <- get_legend(subplots[["cover"]])

subplots[["cover"]] <-
  subplots[["cover"]] +
  theme(legend.position = "none")

subplots[["gemmae_count"]] <-
  combined_morph %>%
  select(count_mean, count_sd, date) %>%
  drop_na %>%
  ggplot(aes(x = date, y = count_mean)) %>%
  shade_years(shading_dates) +
  geom_errorbar(
    aes(ymin=count_mean-count_sd,
        ymax=count_mean+count_sd),
    width=.1,
    color = "dark grey") +
  geom_line() +
  geom_point(color = "black") +
  labs(
    x = "",
    y = expression("Gemmae count"),
    subtitle = "b"
  )

subplots[["gemmae_length"]] <-
  combined_morph %>%
  select(length_mean, length_sd, date) %>%
  drop_na %>%
  ggplot(aes(x = date, y = length_mean)) %>%
  shade_years(shading_dates) +
  geom_errorbar(
    aes(ymin=length_mean-length_sd,
        ymax=length_mean+length_sd),
    width=.1,
    color = "dark grey") +
  geom_line() +
  geom_point(color = "black") +
  labs(
    x = "",
    y = expression(paste("Gemmae len. (", mu, "m)")),
    subtitle = "c"
  )

subplots <- subplots %>%
  map(~ . +
        theme(
          axis.text.x = element_text(
            angle = 30,
            hjust = 1,
            vjust = 0.5,
            margin=margin(-10,0,0,0)
          ),
          plot.margin = margin(0,10,0,10)
        ) +
        scale_x_date(
          date_labels = "%b %Y",
          date_breaks = "6 months",
          limits = c(start_date, end_date)
        )
  )

subplots[1:2] <- subplots[1:2] %>%
  map(~ . + theme(axis.text.x = element_blank())
  )

subplots[2:3] <- subplots[2:3] %>%
  map(~ . + theme(plot.margin = margin(-10,10,0,10))
  )

#' Combine subplots and write out
subplots[[1]] + subplots[[2]] + subplots[[3]] +
  legend +
  plot_layout(ncol = 1, heights = c(1,1,1,0.2))

```

```{r fig-4, fig.width=8, fig.height=7, fig.cap = "Change in microclimate at site with *Pleurosoriopsis makinoi* independent gametophytes (Okutama) and sporophytes (Uratakao) from 2013 to 2014. a--b. Daily light integral (DLI) (mol m^-2 day^-1). c--d. Daily minimum relative humidity (%). e--f. Daily mean temperature (°C). Shaded bars show alternating years. Data collected at Okutama from September 22, 2012 to July 22, 2014 and at Uratakao from August 20, 2013 to November 22, 2014. Some data missing due to equipment malfunction (see Results)."}

#' Extract common start and end dates so x-axis are same across subplots.
start_date <- daily_microclimate %>% pull(date) %>% min
end_date <- daily_microclimate %>% pull(date) %>% max
all_days <- tibble(
  date = seq(start_date, end_date, by = "day")
)

#' Make tibble of years for shading.
shading_dates <-
  tibble(date = seq.Date(from = start_date, to = end_date, by = "day")) %>%
  mutate(year = year(date)) %>%
  # shade by odd years
  filter(year %in% c(2009, 2011, 2013)) %>%
  group_by(year) %>%
  summarize(
    year_start = min(date),
    year_end = max(date)
  )

#' PAR differs enough between sites that automatically set y-axis limits
#' would be different. Use the common min and max between them.
min_par <- daily_microclimate %>% pull(par_total) %>% min
max_par <- daily_microclimate %>% pull(par_total) %>% max

#' Insert NAs into microclimate for each site so that missing days
#' don't get connected by lines.
daily_microclimate_okutama <- filter(daily_microclimate, site == "okutama") %>%
  right_join(all_days)

daily_microclimate_takao <- filter(daily_microclimate, site == "uratakao") %>%
  right_join(all_days)

daily_microclimate_with_na <- bind_rows(
  daily_microclimate_okutama, daily_microclimate_takao)

#' Make list of subplots.
subplots <- list()

subplots[[1]] <-
  daily_microclimate_okutama %>%
  ggplot(aes(x = date, y = par_total)) %>%
  shade_years(shading_dates) +
  geom_line() +
  scale_y_continuous(limits = c(min_par, max_par)) +
  labs(y = expression(paste("DLI (", mol~m^-2~day^-1, ")", sep = "") ),
       x = "",
       title = "Okutama",
       subtitle = "a")

subplots[[2]] <-
  daily_microclimate_takao %>%
  ggplot(aes(x = date, y = par_total)) %>%
  shade_years(shading_dates) +
  geom_line() +
  scale_y_continuous(limits = c(min_par, max_par)) +
  labs(y = "",
       x = "",
       title = "Uratakao",
       subtitle = "b")

subplots[[3]] <-
  daily_microclimate_okutama %>%
  ggplot(aes(x = date, y = rh_min)) %>%
  shade_years(shading_dates) +
  geom_line() +
  scale_y_continuous(limits = c(0, 100)) +
  labs(y = "Min. Rel. Humidity (%)",
       x = "",
       subtitle = "c")

subplots[[4]] <-
  daily_microclimate_takao %>%
  ggplot(aes(x = date, y = rh_min)) %>%
  shade_years(shading_dates) +
  geom_line() +
  scale_y_continuous(limits = c(0, 100)) +
  labs(y = "",
       x = "",
       subtitle = "d")

subplots[[5]] <-
  daily_microclimate_okutama %>%
  ggplot(aes(x = date, y = temp_mean)) %>%
  shade_years(shading_dates) +
  geom_line() +
  labs(y = "Mean Temp. (°C)",
       x = "",
       subtitle = "e")

subplots[[6]] <-
  daily_microclimate_takao %>%
  ggplot(aes(x = date, y = temp_mean)) %>%
  shade_years(shading_dates) +
  geom_line() +
  labs(y = "",
       x = "",
       subtitle = "f")

#' Apply common formatting.
subplots <- subplots %>%
  map(~ . +
        theme(
          axis.text.x = element_text(
            angle = 30,
            hjust = 1,
            vjust = 0.5,
            margin=margin(-10,0,0,0)
          )
        ) +
        scale_x_date(
          date_labels = "%b %Y",
          date_breaks = "6 months",
          limits = c(start_date, end_date)
        )
  )

#' Remove x-axis lables from upper plots.
subplots[1:4] <- subplots[1:4] %>%
  map(~ . + theme(axis.text.x = element_blank())
  )

#' Close gaps on top margins for lower plots.
subplots[c(3,5)] <- subplots[c(3,5)] %>%
  map(~ . + theme(plot.margin = margin(-10,10,0,10))
  )

#' Close gaps on left side for RHS plots.
subplots[c(2,4,6)] <- subplots[c(2,4,6)] %>%
  map(~ . + theme(plot.margin = margin(-10,10,0,-40))
  )

#' Combine plots and write out.
subplots[[1]] + subplots[[2]] + subplots[[3]] +
  subplots[[4]] + subplots[[5]] + subplots[[6]] +
  plot_layout(ncol = 2)

```

```{r fig-5, evaluate = FALSE, fig.width=7, fig.height=8, fig.cap="Density plots showing difference between microclimate at *Pleurosoriopsis makinoi* independent gametophyte site (Okutama) and sporophyte site (Uratakao). Significant differences (paired t-test) shown by asterisks: \\*\\*\\*, *p* < 0.001; \\*\\*, *p* < 0.01; \\*, *p* < 0.05.  Y-axis is relative probability density; curves scaled so that total probability of each scales to 1. Vertical lines shown medians."}
#' Extract p-values for plotting, and put into data frame
#' with xvals specifying season and yvals specifying height
#' to plot asterisks.
p_vals <-
  t_test_results %>%
  mutate(pval = case_when(
    pval < 0.0001 ~ "***",
    pval < 0.001 ~ "**",
    pval < 0.05 ~ "*",
    TRUE ~ ""
  )) %>%
  select(season, var, pval) %>%
  spread(var, pval) %>%
  rename(
    par_total_pval = par_total,
    rh_min_pval = rh_min,
    temp_mean_pval = temp_mean
  )

# Simple function to calculate maximum y-position for plotting asterisks.
# Should be just above the mean + error bar height (1 sd).
maxy <- function(x) {mean(x, na.rm = TRUE) + sd(x, na.rm = TRUE)}

pval_plot_data <-
  paired_microclimate %>%
  select(-date, -month) %>%
  group_by(season, site) %>%
  summarize_all(
    maxy
  ) %>%
  group_by(season) %>%
  summarize_all(
    max
  ) %>%
  select(-site) %>%
  left_join(p_vals)

#' Format data for plotting: need columns to be mean and sd, by site and season.
paired_means <-
  paired_microclimate %>%
  select(-date) %>%
  group_by(site, season) %>%
  summarize_all(
    list(~mean(., na.rm = TRUE),
         ~sd(., na.rm = TRUE))
  )

#' Make plot including t-test on three microclimatic variables.
#'
#' The final plot includes multiple subplots. Each is a ridge plot showing
#' the distribution of one of the variables between the two sites, with the
#' results of a paired t-test showing if they are different or not.
#'
#' Need to manually set x ranges to keep x-axis scale the same across plots
x_ranges <- tibble(var = selected_vars) %>%
  mutate(
    min = map_dbl(var, ~ pull(paired_microclimate, .) %>% min ),
    max = map_dbl(var, ~ pull(paired_microclimate, .) %>% max )
  )

# light blue / light green colorblind-friendly from colorbrewer
# paired_colors <- c("#a6cee3", "#b2df8a")
# OR dark blue / dark green colorblind-friendly from colorbrewer
# paired_colors <- c("#1f78b4", "#33a02c")
# OR dark grey / light grey for bw
paired_colors <- c("grey55", "grey90")

subplots <- list()

#' Make ridge plots by mapping across nested datasets.
plot_results <-
  nested_paired_microclimate %>%
  right_join(x_ranges) %>%
  mutate(min = min - 0.05 * min,
         max = max + 0.05 * max,
         plot = pmap(list(var, data, x_min = min, x_max = max),
                     plot_ridges,
                     y_var = "site"),
         pval = map2(var, data, run_t_test),
         pval = case_when(
           pval < 0.0001 ~ "***",
           pval < 0.001 ~ "**",
           pval < 0.05 ~ "*",
           TRUE ~ ""
         ),
         plot = map2(
           plot, pval,
           ~ .x + annotate("text", Inf, Inf,
                           label = .y, hjust = 1, vjust = 1))

         # If needed, add annotation for season to make sure subplots
         # are positioned correctly

         # ,plot = map2(
         #   plot, season,
         #   ~ .x + annotate("text", -Inf, Inf,
         #                   label = .y, hjust = -1, vjust = 1)),
         # season = as_factor(season),
         # season = fct_relevel(season, "winter", "spring", "summer", "fall")

  ) %>%
  mutate(
    season = factor(season, levels = c("winter", "spring", "summer", "fall"))
  ) %>%
  arrange(season, var)

subplots <- plot_results$plot

#' Manually some subplots to remove redundant labels, etc.
subplots[1:9] <- map(subplots[1:9],
                     ~ . + theme(axis.title.x = element_blank(),
                                 axis.text.x = element_blank() ))


subplots[1:12] <- map(subplots[1:12],
                      ~ . + scale_fill_manual(
                        values = paired_colors,
                        name="Site",
                        breaks=c("okutama", "uratakao"),
                        labels=c("Okutama", "Uratakao"))
)

subplots[[10]] <- subplots[[10]] +
  labs(x = expression(paste("DLI (", mol~m^-2~day^-1, ")", sep = "")) )

subplots[[11]] <- subplots[[11]] +
  labs(x = expression("Min. Rel. Humidity (%)"))

subplots[[12]] <- subplots[[12]] +
  labs(x = expression("Mean Temp. (°C)"))

#' Add titles.
subplots <- map(subplots, ~ . + labs(title = "") )

subplots[[2]] <- subplots[[2]] +
  labs(title = "Winter\n")

subplots[[5]] <- subplots[[5]] +
  labs(title = "Spring\n")

subplots[[8]] <- subplots[[8]] +
  labs(title = "Summer\n")

subplots[[11]] <- subplots[[11]] +
  labs(title = "Fall\n")

#' Extract legend
legend <- paired_microclimate %>%
  ggplot(aes(x = temp_mean, y = site, fill = site)) +
  geom_density_ridges() +
  theme(legend.position = "bottom") +
  scale_fill_manual(
    values = paired_colors,
    name="Site",
    breaks=c("okutama", "uratakao"),
    labels=c("Okutama", "Uratakao"))

legend <- get_legend(legend)

#' Combine and write out.
subplots[[1]] +  subplots[[2]] +  subplots[[3]] + # Winter
subplots[[4]] +  subplots[[5]] +  subplots[[6]] + # Spring
subplots[[7]] +  subplots[[8]] +  subplots[[9]] + # Summer
subplots[[10]] + subplots[[11]] + subplots[[12]] + # Fall
plot_spacer() + legend + plot_spacer() + plot_layout(ncol = 3, heights = c(1, 1, 1, 1, 0.2))

```

```{r fig-6, fig.width=7, fig.height=8, fig.cap="Relationships between growth of *Pleurosoriopsis makinoi* gametophytes (total cover, mean gemma number, and mean gemma length) and microclimate at the Okutama site. Microclimatic variables are monthly means of daily light integral (DLI) (mol m^-2^ day^-1^), daily minimum relative humidity (%), and daily mean temperature (°C). Grey bars show SD (except for total cover, which is a sum, not a mean). Significant relationships (linear model) shown with a line."}
#' Loop through all combination of dep and indep variables and their models,
#' and make plots for each.
var_grid <- make_var_grid (c("rh_min", "par_total", "temp_mean"))

subplots <-
  pmap(var_grid,
       quick_plot_with_stats,
       plot_data = combined_monthly_morph,
       model_data = linear_model_results)

#' Add axis titles
#'
#' First, remove all default titles, then add formatted titles.
subplots <- map(subplots, ~ . + labs(x = "", y = ""))
subplots[[1]] <- subplots[[1]] + labs(y = expression("Total cover ("~cm^2~")"))
subplots[[4]] <- subplots[[4]] + labs(y = "Gemmae count")
subplots[[7]] <- subplots[[7]] +
  labs(
    y = expression(paste("Gemmae len. (", mu, "m)", sep = "")),
    x = expression(paste("DLI (", mol~m^-2~day^-1, ")", sep = ""))
  )
subplots[[8]] <- subplots[[8]] + labs(x = expression("Min. Rel. Humidity (%)"))
subplots[[9]] <- subplots[[9]] + labs(x = expression("Mean Temp. (°C)"))

#' Set amount of gap to close on left and top
left_close <- -15
top_close <- -10

#' Close gaps on left only: 2,3
subplots[2:3] <- subplots[2:3] %>%
  map(~ . + theme(plot.margin = margin(0,0,0,left_close))
  )

#' Close gaps on left side and top: 5,6,8,9
subplots[c(5,6,8,9)] <- subplots[c(5,6,8,9)] %>%
  map(~ . + theme(plot.margin = margin(top_close,0,0,left_close))
  )

#' Close gaps on top only: 4,7
subplots[c(4,7)] <- subplots[c(4,7)] %>%
  map(~ . + theme(plot.margin = margin(top_close,0,0,0))
  )

#' Add legend for middle bottom plot, combine, and write out.
subplots[[8]] <- subplots[[8]] +
  theme(legend.position="bottom") +
  labs(color = "Month", shape = "Year")

subplots[[1]] + subplots[[2]] + subplots[[3]] +
  subplots[[4]] + subplots[[5]] + subplots[[6]] +
  subplots[[7]] + subplots[[8]] + subplots[[9]] +
  plot_layout(ncol = 3)

#' Black and white version if needed.
# subplots <-
#   map(subplots,
#       ~ . + scale_color_grey()
#   )
#
# subplots[[1]] + subplots[[2]] + subplots[[3]] +
#   subplots[[4]] + subplots[[5]] + subplots[[6]] +
#   subplots[[7]] + subplots[[8]] + subplots[[9]] +
#   plot_layout(ncol = 3)
```

```{r fig-7, fig.width=7, fig.height=4, fig.cap = "Sizes of gametophytes with antheridia (male), archegonia (female), or no reproductive structures (asexual). a. Length (mm). b. Width (mm). Boxplots show medians and extend to the first and third quartile; whiskers extend to 1.5× the interquartile range. A single archegoniate (female) outlier with length of 46.6 mm was removed prior to plotting. Points randomly scattered within each x-axis category to reduce overlap."}
#' Note there is one outlier with length 46 mm,
#' then the rest all less than 19 mm. We will remove this
#' before plotting and make a note in the figure caption.
gameto_size_no_outlier <-
  gameto_size %>%
  # Change sex to uppercase for nicer plotting
  mutate(sex = stringr::str_to_title(sex)) %>%
  # Remove one outlier
  filter(length < 40)

#' Make jitter plots overlayed with box plots. Make two
#' subplots: one for width and one for length, then combine.

# Use colorbrewer palette: 3 diverging, color-blind friendly colors
cols <- c("#1b9e77","#d95f02","#7570b3")

# Try to make sure dots stay in the same place, but doesn't work
set.seed(1000)

# Color version

subplots <- list()

subplots[["length"]] <-
  ggplot(gameto_size_no_outlier, aes(y = length, x = sex, color = sex)) +
  labs(
    subtitle = "a",
    y = "Length (mm)",
    x = "",
    color = "")

subplots[["width"]] <-
  ggplot(gameto_size_no_outlier, aes(y = width, x = sex, color = sex)) +
  labs(
    subtitle = "b",
    y = "Width (mm)",
    x = "")

subplots <-
  map(subplots,
      ~ . +
        geom_jitter(size = 3) +
        geom_boxplot(
          outlier.shape = NA,
          fill = NA,
          color = "dark grey",
          alpha = 0.7
        ) +
        scale_color_manual(values = cols) +
        theme(legend.position = "none")
  )

subplots[[1]] + subplots[[2]] + plot_layout(ncol = 2)

# Black and white version
# subplots <-
#   map(subplots,
#       ~ . +
#         scale_color_grey() +
#         theme(legend.position = "none")
#   )

```

```{r fig-8, fig.width=7, fig.height=7, fig.cap = "A comparison of vegetation of the independent gametophyte site with that of a sporophyte site of Pleurosoriopsis. a--b. Typical vegetation of sporophyte habitat covered with summer-green broad-leaved trees (Chichibu-shi, Saitama Pref., Japan). c--d. Vegetation of the independent gametophyte site in Unazawa Valley covered with evergreen cypress trees. a, c. in August, b, d. in December.  "}
knitr::include_graphics(here("images/Fig8_vegetation.pdf"))
```

```{r, results='asis', echo = FALSE, eval=(knitr::opts_knit$get('rmarkdown.pandoc.to') == 'latex')}
# Add pagebreak if rendering to pdf
cat('\\pagebreak')
```

## Tables

Table 1. Results of paired *t*-test for difference in means between microclimatic variables (daily light integral [DLI], daily minimum relative humidity, and daily mean temperature) at two sites in Japan hosting independent gametophytes or sporophytes of *Pleurosoriopsis makinoi* (Okutama and Uratakao, respectively) by season.

```{r table-1}
#' Format Table 1 and output to RTF using `gt` package.
t_test_results %>%
  arrange(factor(var, levels = order_vars)) %>%
  mutate(season = stringr::str_to_title(season)) %>%
  mutate(var = case_when(
    var == "rh_min" ~    "Min. Rel. Humidity",
    var == "par_total" ~ "DLI",
    var == "temp_mean" ~ "Mean Temp."
  )) %>%
  group_by(season) %>%
  gt() %>%
  fmt_number(
    columns = vars(estimate, statistic, conf.low, conf.high,
                   okutama_mean, okutama_sd,
                   uratakao_mean, uratakao_sd)
  ) %>%
  cols_move(
    columns = vars(okutama_mean, okutama_sd),
    after = vars(conf.high)
  ) %>%
  tab_spanner(
    label = "Okutama",
    columns = vars(okutama_mean, okutama_sd)
  ) %>%
  tab_spanner(
    label = "Uratakao",
    columns = vars(uratakao_mean, uratakao_sd)
  ) %>%
  cols_label(
    var = "",
    statistic = md("*t*"),
    pval = md("*p*"),
    parameter = "df",
    conf.low = "Lower",
    conf.high = "Upper",
    okutama_mean = "Mean",
    okutama_sd = "SD",
    uratakao_mean = "Mean",
    uratakao_sd = "SD"
  ) %>%
  tab_style(
    style = cells_styles(
      bkgd_color = "white"),
    locations = cells_data(
      columns = everything()
    )
  ) %>%
  tab_header(
    title = "Table 1",
    subtitle = md("&nbsp;") # work-around for https://github.com/rstudio/gt/issues/197
  )

```

```{r, results='asis', echo = FALSE, eval=(knitr::opts_knit$get('rmarkdown.pandoc.to') == 'latex')}
# Add pagebreak if rendering to pdf
cat('\\pagebreak')
```

Table 2. Results of linear models testing for relationship between growth traits (total cover, count of gemmae, and length of gemmae) as dependent variable on microclimate (monthly averages of daily light integral [DLI], daily minimum relative humidity, and daily mean temperature) as independent variable in gametophytes of *Pleurosoriopsis makinoi* at Unazawa Valley, Okutama-machi, Tokyo, Japan. Results presented by independent variable.

```{r table-2}
cols_numeric <- c("r.squared", "adj.r.squared", "sigma", "statistic", "logLik", "AIC", "BIC")

linear_model_results %>%
  select(-fits) %>%
  arrange(factor(x_var, levels = order_vars)) %>%
  mutate(
    x_var = case_when(
      x_var == "rh_min" ~    "Min. Rel. Humidity",
      x_var == "par_total" ~ "DLI",
      x_var == "temp_mean" ~ "Mean Temp."
    ),
    y_var = case_when(
      y_var == "total_cover" ~ "Total cover",
      y_var == "count_mean" ~ "Gemmae count",
      y_var == "length_mean" ~ "Gemmae length"
    )
  ) %>%
  group_by(x_var) %>%
  gt() %>%
  cols_move_to_start(vars(y_var)) %>%
  cols_move(
    columns = vars(df.residual),
    after = vars(df)
  ) %>%
  fmt_number(
    columns = colnames(.) %in% cols_numeric
  ) %>%
  fmt_number(
    columns = vars(p.value),
    decimals = 3
  ) %>%
  fmt_number(
    columns = vars(deviance),
    decimals = 0
  ) %>%
  cols_label(
    y_var = "",
    r.squared = md("r2"),
    adj.r.squared = md("adj. r2"),
    statistic = md("*t*"),
    p.value = md("*p*"),
    df.residual = "df resid."
  ) %>%
  tab_style(
    style = cells_styles(
      bkgd_color = "white"),
    locations = cells_data(
      columns = everything()
    )
  ) %>%
  tab_header(
    title = "Table 2",
    subtitle = md("&nbsp;") # work-around for https://github.com/rstudio/gt/issues/197
  )

```