diff --git a/R/functions.R b/R/functions.R
index da84b37..aab3c5f 100644
--- a/R/functions.R
+++ b/R/functions.R
@@ -2252,6 +2252,79 @@ relabel_bioregions_by_lat <- function(bioregions, shape_ferns, cutoff = 1) {
 
 }
 
+
+#' Relabel bioregions in order of latidudinal position
+#' 
+#' Keep the old labels in the dataframe.
+#' 
+#' Should eventually use this as a replacement for relabel_bioregions_by_lat()
+#'
+#' @param bioregions Bioregions of Japanese ferns; output of combine_bioregions()
+#' @param shape_ferns Simple features data frame; grid of fern richness, abundance, and redundancy
+#' @param cutoff Cutoff value to use for filtering 'major' bioregions. Only bioregions
+#' with more grid-cells than the cutoff value will be relabeled.
+#'
+#' @return Dataframe; bioregions relabeled in descending order by mean latitude
+#' 
+relabel_bioregions_by_lat_keep_old <- function(bioregions, shape_ferns, cutoff = 1) {
+
+  # Summarize bioregions data: cells per cluster, mean latitude of each bioregion
+  bioregions_summary_all <-
+    shape_ferns %>%
+    left_join(bioregions, by = "grids") %>%
+    # Set CRS to JDG2000 = EPSG code 4612 so that lat and long are in degrees
+    st_transform(4612) %>%
+    sf_add_centroids() %>%
+    st_drop_geometry() %>%
+    select(contains("cluster"), lat, grids) %>%
+    as_tibble() %>%
+    pivot_longer(
+      names_to = "cluster_type",
+      values_to = "cluster",
+      -c(grids, lat)) %>%
+    group_by(cluster_type, cluster) %>%
+    summarize(mean_lat = mean(lat), n = n(), .groups = "drop")
+
+  # Map original cluster labels to new cluster labels
+  cluster_mapping <-
+    bioregions_summary_all %>%
+    # New labels are in descending order by latitude starting with major regions
+    # (n grid-cells above cutoff)
+    mutate(major_region = n > cutoff) %>%
+    arrange(cluster_type, desc(major_region), desc(mean_lat)) %>%
+    group_by(cluster_type) %>%
+    mutate(new_cluster = sort(cluster)) %>%
+    ungroup() %>%
+    mutate(new_cluster = case_when(
+      major_region == FALSE ~ "Other",
+      TRUE ~ new_cluster
+    ))
+
+  # Convert the old bioregion labels to new ones
+  bioregions %>%
+    left_join(
+      filter(cluster_mapping, cluster_type == "phylo_cluster") %>% 
+        select(cluster, new_phylo_cluster = new_cluster),
+      by = c("phylo_cluster" = "cluster")) %>%
+    left_join(
+      filter(cluster_mapping, cluster_type == "taxonomic_cluster") %>% 
+        select(cluster, new_taxonomic_cluster = new_cluster),
+      by = c("taxonomic_cluster" = "cluster")) %>%
+    rename(
+      old_taxonomic_cluster = taxonomic_cluster,
+      old_phylo_cluster = phylo_cluster,
+    ) %>%
+    rename(
+      taxonomic_cluster = new_taxonomic_cluster, 
+      phylo_cluster = new_phylo_cluster
+    ) %>%
+    # Make sure the number of grid cells per cluster is the same 
+    # (only the labels changed)
+    assert(not_na, taxonomic_cluster, phylo_cluster) %>%
+    verify(nrow(.) == nrow(bioregions))
+
+}
+
 # Conservation ----
 
 #' Crop biodiversity data to protected areas in Japan
diff --git a/_targets.R b/_targets.R
index ccde252..520f733 100644
--- a/_targets.R
+++ b/_targets.R
@@ -334,6 +334,9 @@ tar_plan(
   bioregion_cutoff = 2, # later, lump bioregions into "other" if they have fewer than this number of grid-cells # nolint
   bioregions = relabel_bioregions_by_lat(
     bioregions_raw, shape_ferns, cutoff = bioregion_cutoff),
+  # - once more, retaining original labels for plotting
+  bioregions_all_labels = relabel_bioregions_by_lat_keep_old(
+    bioregions_raw, shape_ferns, cutoff = bioregion_cutoff),
 
   # Traits analysis ----
 
diff --git a/ms/SI.Rmd b/ms/SI.Rmd
index 001415c..f9aa166 100644
--- a/ms/SI.Rmd
+++ b/ms/SI.Rmd
@@ -90,6 +90,7 @@ tar_load(c(
   "deer_range",
   "biodiv_ferns_spatial_lumped",
   "bioregion_cutoff",
+  "bioregions_all_labels",
   "ms_functions"
 ))
 ```
@@ -406,7 +407,7 @@ a + b + c +
   theme(plot.tag = element_text(face = "bold"))
 ```
 
-**`r s_figure("grain-size")`**. Effect of grain size on sampling redundancy. (A), Percent of grid cells with redundancy > 0.5 by grain size (side length of square grid cells set to 10 km, 20 km, 30 km, or 40 km). (B), Percent of grid cells with redundancy = 0 by grain size. (C), Distribution of redundancy values by grain size.
+**`r s_figure("grain-size")`**. Effect of grain size on sampling redundancy. (A) Percent of grid cells with redundancy > 0.5 by grain size (side length of square grid cells set to 10 km, 20 km, 30 km, or 40 km). (B) Percent of grid cells with redundancy = 0 by grain size. (C) Distribution of redundancy values by grain size.
 
 \clearpage
 
@@ -437,18 +438,18 @@ plot_env <- function(part, label, var_select, fill_lab = "", labels = ggplot2::w
     )
 }
 
-a <- plot_env("a", "Area", "lat_area", "km^2")
-b <- plot_env("b", "% apomictic taxa", "percent_apo", labels = function(x) scales::percent(x, accuracy = 1))
-c <- plot_env("c", "Precipitation", "precip", "mm")
-d <- plot_env("d", "Precip. season.", "precip_season" , "")
-e <- plot_env("e", "Temperature", "temp", "°C", function(x) x/10)
-f <- plot_env("f", "Temp. season.", "temp_season") +
+a <- plot_env("A", "Area", "lat_area", "km^2")
+b <- plot_env("B", "% apomictic taxa", "percent_apo", labels = function(x) scales::percent(x, accuracy = 1))
+c <- plot_env("C", "Precipitation", "precip", "mm")
+d <- plot_env("D", "Precip. season.", "precip_season" , "")
+e <- plot_env("E", "Temperature", "temp", "°C", function(x) x/10)
+f <- plot_env("F", "Temp. season.", "temp_season") +
   annotation_scale(location = "tl", width_hint = 0.2, pad_y = unit(0.8, "cm"))
 
 a + b + c + d + e + f + plot_layout(ncol = 2, nrow = 3)
 ```
 
-**`r s_figure("envir")`**. Reproductive mode and environmental data in 20 km × 20 km grid cells across Japan. (A), rolling mean area (km^2^) in 20 km latitudinal bands.  (B), % apomictic taxa. (C), annual precipitation. (D), percipitation seasonality (coefficient of variation). (E), mean annual temperature. (F), temperature seasonality (standard deviation × 100). Data excluded from models not shown in this and subsequent plots unless otherwise mentioned.
+**`r s_figure("envir")`**. Reproductive mode and environmental data in 20 km × 20 km grid cells across Japan. (A) rolling mean area (km^2^) in 20 km latitudinal bands. (B) % apomictic taxa. (C) annual precipitation. (D) percipitation seasonality (coefficient of variation). (E) mean annual temperature. (F) temperature seasonality (standard deviation × 100). Data excluded from models not shown in this and subsequent plots unless otherwise mentioned.
 
 \clearpage
 
@@ -513,7 +514,7 @@ a + b +
   plot_layout(nrow = 1)
 ```
 
-**`r s_figure("abundance")`**. Observed number of specimens (A) and sampling redundancy (B) per 20 km grid cell in the ferns of Japan.
+**`r s_figure("abundance")`**. Observed number of specimens (A) and sampling redundancy (B) per 20 km × 20 km grid cell in the ferns of Japan.
 
 \clearpage
 
@@ -565,7 +566,7 @@ a + b + c +
   theme(plot.tag = element_text(face = "bold"))
 ```
 
-**`r s_figure("raw-lat")`**. (A) Raw taxonomic diversity (B), phylogenetic diversity, and (C) functional diversity of the ferns of Japan plotted by latitude. For (B) and (C), raw branchlengths were transformed to relative values by dividing each branchlength by the sum of all branchlengths. Trendlines fit with general additive model using `r func("geom_smooth")` function in the `r pack("ggplot2")` R package [@Wickham2016].
+**`r s_figure("raw-lat")`**. (A) Raw taxonomic diversity, (B) phylogenetic diversity, and (C) functional diversity of the ferns of Japan plotted by latitude. For (B) and (C), raw branchlengths were transformed to relative values by dividing each branchlength by the sum of all branchlengths. Trendlines fit with general additive model using `r func("geom_smooth")` function in the `r pack("ggplot2")` R package [@Wickham2016].
 
 \clearpage
 
@@ -596,7 +597,7 @@ a + b + c + d +
   theme(plot.tag = element_text(face = "bold"))
 ```
 
-**`r s_figure("raw-biplot")`**. Relationships between observed functional and phylogenetic diversity and taxonomic richness in the ferns of Japan. (A), Phylogenetic diversity (PD). (B), Relative phylogenetic diversity (RPD). (C), Functional diversity (FD). (D), Relative functional diversity (RFD). Trendlines in (a, c) fit with general additive model using `r func("geom_smooth")` function in the `r pack("ggplot2")` R package [@Wickham2016].
+**`r s_figure("raw-biplot")`**. Relationships between observed phylogenetic and functional diversity and taxonomic richness in the ferns of Japan. (A) Phylogenetic diversity (PD). (B) Relative phylogenetic diversity (RPD). (C) Functional diversity (FD). (D) Relative functional diversity (RFD). Trendlines in (A, C) fit with general additive model using `r func("geom_smooth")` function in the `r pack("ggplot2")` R package [@Wickham2016].
 
 \clearpage
 
@@ -613,7 +614,7 @@ ggplot() +
   theme(legend.title = element_blank())
 ```
 
-**`r s_figure("endemism-restricted")`**. Phylogenetic endemism of the ferns of Japan measured using CANAPE (categorical analysis of neo- and paleo-endemism), restricted dataset including only taxa endemic to Japan. 'paleo' indicates paleoendemic areas (those with an overabundance of long branches); 'neo' indicates neodendemic areas (those with an overabundance of short branches); 'mixed' indicates significantly endemic sites that have neither an overabundance of short nor long branches; 'super' indicates highly significantly endemic sites that have neither an overabundance of short nor long branches. For details, see Methods.
+**`r s_figure("endemism-restricted")`**. Phylogenetic endemism of the ferns of Japan measured using CANAPE (categorical analysis of neo- and paleo-endemism), restricted dataset including only taxa endemic to Japan. 'neo' indicates neodendemic areas (those with an overabundance of short branches); 'paleo' indicates paleoendemic areas (those with an overabundance of long branches); 'mixed' indicates significantly endemic sites that have neither an overabundance of short nor long branches; 'super' indicates highly significantly endemic sites that have neither an overabundance of short nor long branches. For details, see Materials and Methods.
 
 \clearpage
 
@@ -670,7 +671,7 @@ a + b +
   theme(plot.tag = element_text(face = "bold"))
 ```
 
-**`r s_figure("k-plot")`**. Selection of *k* for bioregions. (A), taxonomic bioregions. (B), phylogenetic bioregions. Color indicates optimal *k* selected selected by the `r func("optimal_phyloregion")` function in the R package `r pack("phyloregion")` (black, non-optimal; red, optimal).
+**`r s_figure("k-plot")`**. Selection of *k* for bioregions. (A) taxonomic bioregions. (B) phylogenetic bioregions. Color indicates optimal *k* selected selected by the `r func("optimal_phyloregion")` function in the R package `r pack("phyloregion")` (black, non-optimal; red, optimal).
 
 \clearpage
 
@@ -678,18 +679,29 @@ a + b +
 # Extract dendrogram from phyloregions analysis, convert to phylogeny
 dendro_tax <- attributes(k_taxonomy$df)$meta$hclust.obj %>% as.phylo() %>% ladderize()
 
-# Map taxonomic clusters to grid cells
+# Map taxonomic clusters to grid cells. Include original ("old") clusters
+# so that grouping of "other" works correctly even if not monophyletic
 clusters_tax <-
-  biodiv_ferns_spatial %>%
-  as_tibble() %>%
-  select(grids, cluster = taxonomic_cluster)
+bioregions_all_labels %>%
+  # Filter to grid-cells used in modeling
+  inner_join(
+    select(biodiv_ferns_spatial, grids),
+    by = "grids"
+  ) %>%
+  add_count(old_taxonomic_cluster) %>%
+  mutate(
+    taxonomic_cluster = case_when(
+      n <= bioregion_cutoff ~ paste(taxonomic_cluster, old_taxonomic_cluster, sep = "_"),
+      TRUE ~ taxonomic_cluster
+      )
+  ) %>%
+  select(-n) %>%
+  select(grids, cluster = taxonomic_cluster) %>%
+  arrange(grids)
 
 cluster_nodes_tax_tbl <- map_df(unique(clusters_tax$cluster), ~map_cluster_to_nodes(dendro_tax, clusters_tax, .)) %>%
   # Lump regions with fewer grid-cells than cutoff
-  mutate(cluster = as.factor(cluster) %>% fct_lump_min(bioregion_cutoff))
-
-# Make sure cutoff matches caption
-assert_that(bioregion_cutoff == 2, msg = "Bioregions cutoff for lumping not 2!")
+  mutate(cluster = as.factor(cluster) %>% fct_lump_min(bioregion_cutoff + 1))
 
 # Use the ggtree::`%<+%` operator to map the bootstrap values onto the ggtree object
 ggtree(dendro_tax, size = 0.15) %<+% 
@@ -702,22 +714,36 @@ ggtree(dendro_tax, size = 0.15) %<+%
   theme(plot.subtitle = element_markdown())
 ```
 
-**`r s_figure("dendrogram")`**. Bioregion dendrograms. Tips are sites (names not shown). Dendrograms generated by clustering taxonomic (Sørensen) or phylogenetic (PhyloSor) distances between sites (see Materials and Methods). Color shows bioregion; bioregions not consisting of more than `r english2(bioregion_cutoff)` grid-cells each are lumped into the "Other" category. (A), taxonomic bioregions. (B), phylogenetic bioregions. Dendrograms visualized with the `r pack("ggtree")` R package [@Yu2017].
+**`r s_figure("dendrogram")`**. Bioregion dendrograms. Tips are sites (names not shown). Dendrograms generated by clustering taxonomic (Sørensen) or phylogenetic (PhyloSor) distances between sites (see Materials and Methods). Color shows bioregion; bioregions not consisting of more than `r english2(bioregion_cutoff)` grid-cells each are lumped into the "Other" category. (A) taxonomic bioregions. (B) phylogenetic bioregions. Dendrograms visualized with the `r pack("ggtree")` R package [@Yu2017].
 \clearpage
 
 ```{r dendrogram-b, fig.width = 7, fig.height = 9}
 # Extract dendrogram from phyloregions analysis, convert to phylogeny
 dendro_phy <- attributes(k_phylogeny$df)$meta$hclust.obj %>% as.phylo() %>% ladderize()
 
-# Map phylogenetic clusters to grid cells
+# Map phylogenetic clusters to grid cells. Include original ("old") clusters
+# so that grouping of "other" works correctly even if not monophyletic
 clusters_phy <-
-  biodiv_ferns_spatial %>%
-  as_tibble() %>%
-  select(grids, cluster = phylo_cluster)
+bioregions_all_labels %>%
+  # Filter to grid-cells used in modeling
+  inner_join(
+    select(biodiv_ferns_spatial, grids),
+    by = "grids"
+  ) %>%
+  add_count(old_phylo_cluster) %>%
+  mutate(
+    phylo_cluster = case_when(
+      n <= bioregion_cutoff ~ paste(phylo_cluster, old_phylo_cluster, sep = "_"),
+      TRUE ~ phylo_cluster
+      )
+  ) %>%
+  select(-n) %>%
+  select(grids, cluster = phylo_cluster) %>%
+  arrange(grids)
 
 cluster_nodes_phy_tbl <- map_df(unique(clusters_phy$cluster), ~map_cluster_to_nodes(dendro_phy, clusters_phy, .)) %>%
   # Lump regions with fewer grid-cells than cutoff
-  mutate(cluster = as.factor(cluster) %>% fct_lump_min(bioregion_cutoff))
+  mutate(cluster = as.factor(cluster) %>% fct_lump_min(bioregion_cutoff + 1))
 
 # Use the ggtree::`%<+%` operator to map the bootstrap values onto the ggtree object
 ggtree(dendro_phy, size = 0.15) %<+% 
@@ -784,7 +810,7 @@ f <- plot_genus("F", "Diplazium") +
 a + b + c + d + e + f + plot_layout(ncol = 2, nrow = 3)
 ```
 
-**`r s_figure("genus-map")`**. Map of species richness for selected genera of ferns of Japan. (A), *Equisetum* is an example of a possible refugial lineage with greater species richness in northern Japan. (B), *Deparia* is an example of a genus with possible recent radiation in southern Honshu. (C)--(F) show genera with high rates of apomictic taxa (`r s_figure("apo-genus")`), which tend to be greatest in southern Honshu.
+**`r s_figure("genus-map")`**. Map of species richness for selected genera of ferns of Japan. (A) *Equisetum* is an example of a possible refugial lineage with greater species richness in northern Japan. (B) *Deparia* is an example of a genus with possible recent radiation in southern Honshu. (C)--(F) show genera with high rates of apomictic taxa (`r s_figure("apo-genus")`), which tend to be greatest in southern Honshu.
 
 \clearpage
 
@@ -939,10 +965,10 @@ make_protected_plot <- function(part, variable, ...) {
     map_theme()
 }
 
-a <- make_protected_plot("a", "Richness", richness_obs_p_upper > 0.95)
-b <- make_protected_plot("b", "PD", pd_obs_p_upper > 0.95)
-c <- make_protected_plot("c", "FD", fd_obs_p_upper > 0.95)
-d <- make_protected_plot("d", "PE", pe_obs_p_upper > 0.95) +
+a <- make_protected_plot("A", "Richness", richness_obs_p_upper > 0.95)
+b <- make_protected_plot("B", "PD", pd_obs_p_upper > 0.95)
+c <- make_protected_plot("C", "FD", fd_obs_p_upper > 0.95)
+d <- make_protected_plot("D", "PE", pe_obs_p_upper > 0.95) +
   annotation_scale(location = "tl", width_hint = 0.2, pad_y = unit(0.8, "cm"))
 
 a + b + c + d + plot_layout(ncol = 2, nrow = 2) +
@@ -976,10 +1002,10 @@ make_deer_plot <- function(part, variable, ...) {
     map_theme()
 }
 
-a <- make_deer_plot("a", "Richness", richness_obs_p_upper > 0.95)
-b <- make_deer_plot("b", "PD", pd_obs_p_upper > 0.95)
-c <- make_deer_plot("c", "FD", fd_obs_p_upper > 0.95)
-d <- make_deer_plot("d", "PE", pe_obs_p_upper > 0.95) +
+a <- make_deer_plot("A", "Richness", richness_obs_p_upper > 0.95)
+b <- make_deer_plot("B", "PD", pd_obs_p_upper > 0.95)
+c <- make_deer_plot("C", "FD", fd_obs_p_upper > 0.95)
+d <- make_deer_plot("D", "PE", pe_obs_p_upper > 0.95) +
   annotation_scale(location = "tl", width_hint = 0.2, pad_y = unit(0.8, "cm"))
 
 a + b + c + d + plot_layout(ncol = 2, nrow = 2) +
diff --git a/ms/manuscript.Rmd b/ms/manuscript.Rmd
index a3fcd27..e427615 100644
--- a/ms/manuscript.Rmd
+++ b/ms/manuscript.Rmd
@@ -165,17 +165,17 @@ Here, we characterize multiple facets and drivers of biodiversity to understand
 
 **Methods:** We compiled a community dataset of `r nrow(comm_ferns) %>% scales::number(big.mark = ",")` 20 km × 20 km grid-cells including `r ncol(comm_ferns)` taxa based on > `r occ_sum$filtered$n_specimens$value %>% plyr::round_any(1000, f = floor) %>% number(big.mark = ",")` specimen records.
 We combined this with a phylogeny and functional traits to analyze taxonomic, phylogenetic, and functional diversity, and modeled biodiversity metrics in response to environmental factors and reproductive mode.
-Hierarchical clustering was used to analyze bioregions.
+Hierarchical clustering was used to delimit bioregions.
 Conservation status and threats were assessed by comparing the overlap of significantly diverse grid-cells with conservation zones and range maps of native Japanese deer.
 
 **Results:** Taxonomic richness was highest at mid-latitudes.
 Phylogenetic and functional diversity and phylogenetic endemism were highest in small southern islands.
-Relative phylogenetic and functional diversity were high at both high and low latitudes, and low at mid-latitudes.
-Grid-cells were grouped in to `r maj_cluster_count$phylo_cluster$n %>% english2` (phylogenetic) or `r maj_cluster_count$taxonomic_cluster$n %>% english2` (taxonomic) major bioregions.
-Temperature and apomixis were identified as drivers of taxonomic and phylogenetic diversity.
-Conservation status was generally high for grid-cells with significantly high biodiversity, but the threat due to herbivory was greater for taxonomic richness than other metrics.
+Relative phylogenetic and functional diversity were high at high and low latitudes, and low at mid-latitudes.
+Grid-cells were grouped into `r maj_cluster_count$phylo_cluster$n %>% english2` (phylogenetic) or `r maj_cluster_count$taxonomic_cluster$n %>% english2` (taxonomic) major bioregions.
+Temperature and apomixis were identified as drivers of biodiversity patterns.
+Conservation status was generally high for grid-cells with significantly high biodiversity, but the threat due to herbivory by deer was greater for taxonomic richness than other metrics.
 
-**Conclusions:** The integrative approach used here reveals previously undetected patterns and drivers of biodiversity in the ferns of Japan.
+**Conclusions:** Our integrative approach reveals previously undetected patterns and drivers of biodiversity in the ferns of Japan.
 Future conservation efforts should recognize that threats can vary by biodiversity metric and consider multiple multiple metrics when establishing conservation priorities.
 
 <!-- may include title words -->
@@ -244,8 +244,8 @@ Furthermore, incorporating these two frameworks---the categorization of areas in
 
 However, a comprehensive understanding of the relationships between richness and other metrics requires densely sampled data, and such datasets are rare on the regional (country) scale.
 The ferns of Japan are excellent model system because they have been the target of intense botanical interest for several decades and are densely sampled [reviewed in @Ebihara2019b].
-The ferns of Japan include `r total_fern_taxa` native, non-hybrid taxa (including species and varieties) and hundreds of hybrids [@Ebihara2019b].
-The availability of detailed distribution data [distribution maps at the ca. 10 km scale for nearly all species\; @Ebihara2016b; @Ebihara2017], trait data [multiple quantitative and qualitative traits compiled for identification of nearly all species\; @Ebihara2019b], and DNA sequences for > 97% of species [@Ebihara2010; @Ebihara2019b\; all coverage statistics exclude hybrids] make the ferns of Japan an ideal system for investigating the relationships between, and drivers of, multiple dimensions of biodiversity.
+The ferns of Japan include `r total_fern_taxa` native, non-hybrid taxa (including species and varieties) and hundreds of hybrids [nothotaxa\; @Ebihara2019b].
+The availability of detailed distribution data [distribution maps at the ca. 10 km scale for nearly all species\; @Ebihara2016b; @Ebihara2017], trait data [multiple quantitative and qualitative traits compiled for identification of nearly all species\; @Ebihara2019b], and DNA sequences for > 97% of species [@Ebihara2010; @Ebihara2019b\; all coverage statistics exclude nothotaxa] make the ferns of Japan an ideal system for investigating the relationships between, and drivers of, multiple dimensions of biodiversity.
 
 One particularly valuable characteristic of the Japanese fern flora is the availability of data on reproductive mode, which is known for `r repro_count$known$total` native fern taxa excluding hybrids [`r repro_count$known$percent`\; @Ebihara2019b].
 Reproductive mode is likely to affect population-level genetic diversity [@Bengtsson2003], and thereby higher-level biodiversity [@Krueger-Hadfield2019].
@@ -274,14 +274,14 @@ Map data for Japan were downloaded from the Geospatial Information Authority of
 ## Occurrence data
 
 We used a list of native, non-hybrid fern specimens deposited at herbaria in Japan to quantify occurrences [@Ebihara2016b; @Ebihara2017; @Ebihara2019b].
-We chose to use this over other sources (`r eg`, GBIF) because of its high quality, which obviates many cleaning steps otherwise needed when working with large, publicly available datasets.
+We chose to use this over other sources such as the Global Biodiversity Information Facility (GBIF) because of its high quality, which obviates many cleaning steps otherwise needed when working with large, publicly available datasets.
 Furthermore, the overall sampling completeness of this dataset is also high (`r s_figure("sampling-curve")` in Appendix S1; see Supplemental Data with this article).
 The original list comprises `r occ_sum$unfiltered$n_specimens$value %>% scales::number(big.mark = ",")` specimens representing `r occ_sum$unfiltered$n_taxa$value` terminal taxa at the rank of species, subspecies, form, or variety, excluding non-native taxa and nothotaxa.
-For the purposes of analysis, we treated infraspecific taxa as distinct species and hereafter use "taxon" and "species" interchangeably to refer to entities at this rank, unless otherwise indicated.
+We treated infraspecific taxa as distinct species during analysis and hereafter use "taxon" and "species" interchangeably to refer to entities at this rank, unless otherwise indicated.
 All names are standardized to a common taxonomy, the Japan Green List v. 1.01 (http://www.rdplants.org/gl/).
 
 Occurrence records were georeferenced and thoroughly vetted as follows [@Ebihara2016b; @Ebihara2017; @Ebihara2019b].
-For specimens that lacked latitude and longitude data, georeferencing was done by mapping collection site names to a set of standard ca. 10 × 10 km grid squares defined by the Statistics Bureau of Japan (the "second-degree mesh"; http://www.stat.go.jp/english/data/mesh/05.html), and the centroid of the second-degree mesh cell used as the specimen location.
+For specimens that lacked latitude and longitude data, georeferencing was done by mapping collection site names to a set of standard ca. 10 km × 10 km grid squares defined by the Statistics Bureau of Japan (the "second-degree mesh"; http://www.stat.go.jp/english/data/mesh/05.html), and the centroid of the second-degree mesh cell used as the specimen location.
 In the case that the collection site could not be mapped to a single second-degree mesh cell, it was excluded.
 Next, occurrence maps were generated for each taxon showing presence or absence in each second-degree mesh cell.
 In the case that a given taxon appeared insufficiently sampled (`r ie`, not present on the map where it would typically be expected to occur), AE and members of the Nippon Fernist Club (local botanists familiar with the ferns of Japan) searched for additional specimens, which were then added to the list.
@@ -422,6 +422,7 @@ Furthermore, use of community phylogenies generated by only by sampling the spec
 Because one goal of the study is to understand the distribution of paleo- vs. neo-endemic areas (defined in units of time), we require an ultrametric tree.
 Therefore, to obtain a robust, ultrametric tree, we combined the *rbcL* data of @Ebihara2019b (`r n_japan_fern_taxa` taxa) with a globally sampled plastid dataset including all fern sequences on GenBank for four widely-sequenced plastid genes (*atpA*, *atpB*, *rbcL*, and *rps4*; `r n_ftol_fern_taxa` taxa) and `r n_plastome_genes` other coding, single-copy plastid genes extracted from `r n_plastome_fern_taxa` complete fern plastomes (Nitta et al., in prep.) (`r ape::Ntip(plastome_tree) %>% scales::number(big.mark = ",")` OTUs total, including outgroup taxa).
 This gene sampling is comparable to a recent global fern phylogeny that resolved relationships across ca. 4,000 taxa using six plastid markers [@Testo2016a], and the addition of plastome data can be expected to increase support along the backbone.
+
 We conducted maximum-likelihood phylogenetic analysis with `r software("IQ-TREE")` v1.6.12 [@Nguyen2015] under the GTR+I+G model (all sites concatenated into a single data matrix), and assessed node support with 1,000 ultrafast bootstrap replicates [@Minh2013].
 Molecular dating was conducted with `r software("treePL")` v1.0 [@Smith2012] using `r n_fossil_points` fossil calibration points after @Testo2016a, with the exception of treating fossil *Kuylisporites* as belonging to crown *Cyathea*, not *Cyathea* + *Alsophila* [@Loiseau2020].
 The root of the tree (most recent common ancestor of all land plants) was fixed to 475 mya [@Testo2016a].
@@ -1075,7 +1076,7 @@ If high elevation habitats were available, the southern islands would be expecte
 For example, Taiwan, which is just east of the southernmost Japanese islands and harbors high mountains (> 3,000 m), shares many fern species with the main Japanese islands that are not found in southern Japanese islands.
 This supports that the sharp taxonomic turnover between southern and main islands is more likely due to ecological factors than dispersal limitation.
 
-One element of biodiversity that we did not consider here is hybrids (nothotaxa).
+One element of biodiversity that we did not consider here is hybrids.
 The Japanese fern flora includes `r n_hybrid_fern_taxa` nothotaxa [@Ebihara2019b].
 Although hybrids are often thought of as evolutionary "dead-ends", in some cases they are able to interbreed with diploids and thus contribute to diversification [@Barrington1989].
 Future studies should consider the distribution of hybrids in Japan and their possible effects on biodiversity.
@@ -1093,6 +1094,7 @@ The comprehensive molecular sampling of sporophytes in Japan makes this area ide
 The authors thank Members of the Iwasaki Lab at the University of Tokyo and two anonymous reviewers for providing comments that improved the manuscript.
 Shawn Laffan provided advice about CANAPE.
 This study was supported by JSPS KAKENHI Grant Number 16H06279.
+The authors are grateful to members of the Nippon Fernist Club for their cooperation and efforts towards quantifying and conserving the biodiversity of the ferns of Japan.
 
 ## AUTHOR CONTRIBUTIONS
 
@@ -1677,7 +1679,7 @@ ggsave(
 endemism_plot
 ```
 
-**`r figure("endemism")`**. Phylogenetic endemism of the ferns of Japan measured using CANAPE (categorical analysis of neo- and paleo-endemism). 'paleo' indicates areas with an overabundance of rare long branches; 'neo' indicates areas with an overabundance of rare short branches; 'mixed' indicates significantly endemic sites that have neither an overabundance of short nor long branches; 'super' indicates highly significantly endemic sites that have neither an overabundance of short nor long branches. For details, see Materials and Methods.
+**`r figure("endemism")`**. Phylogenetic endemism of the ferns of Japan measured using CANAPE (categorical analysis of neo- and paleo-endemism). 'neo' indicates areas with an overabundance of rare short branches; 'paleo' indicates areas with an overabundance of rare long branches; 'mixed' indicates significantly endemic sites that have neither an overabundance of short nor long branches; 'super' indicates highly significantly endemic sites that have neither an overabundance of short nor long branches. For details, see Materials and Methods.
 
 `r pagebreak_pdf()`
 
diff --git a/ms/references.yaml b/ms/references.yaml
index 2ba478e..fc2b63b 100644
--- a/ms/references.yaml
+++ b/ms/references.yaml
@@ -1833,7 +1833,7 @@ references:
   author:
   - literal: R Core Team
   issued:
-  - year: 2020
+  - year: 2021
   accessed:
     year: 2021
     month: 09
@@ -2440,7 +2440,7 @@ references:
   page: 1451-1456
   source: DOI.org (Crossref)
   title: 'iNEXT: an R package for rarefaction and extrapolation of species diversity
-    ( <span style="font-variant:small-caps;" smallcaps="smallcaps">H</span> ill numbers)'
+    (Hill numbers)'
   title-short: iNEXT
   type: article-journal
   volume: '7'
@@ -2484,9 +2484,6 @@ references:
   container-title: Methods in Ecology and Evolution
   container-title-short: Methods Ecol Evol
   DOI: 10.1111/2041-210x.12628
-  editor:
-  - family: McInerny
-    given: Greg
   ISSN: 2041-210X, 2041-210X
   issue: '1'
   issued:
@@ -2495,8 +2492,8 @@ references:
   language: en
   page: 28-36
   source: DOI.org (Crossref)
-  title: '<span style="font-variant:small-caps;" smallcaps="smallcaps">ggtree</span>
-    : an <span style="font-variant:small-caps;" smallcaps="smallcaps">r</span> package
+  title: '<span style="font-variant:small-caps;" smallcaps="smallcaps">ggtree</span>:
+    an <span style="font-variant:small-caps;" smallcaps="smallcaps">R</span> package
     for visualization and annotation of phylogenetic trees with their covariates and
     other associated data'
   title-short: ''
diff --git a/ms/references_other.yaml b/ms/references_other.yaml
index b3b095a..a2aa36e 100644
--- a/ms/references_other.yaml
+++ b/ms/references_other.yaml
@@ -69,7 +69,7 @@ references:
   author:
   - literal: R Core Team
   issued:
-  - year: 2020
+  - year: 2021
   accessed:
     year: 2021
     month: 09