Rearranging a 'ggplot2' Legend to Mix and Match Different Levels

Maybe this what you are looking for. To group your phyla by domains you could

1. plot each domain separately (by splitting or filtering your dataset) and add a new fill scale for each domain.
2. As this is a lot of manual work (copy and pasting the code 5 times) I use a function to add the layers for each domain and use imap to loop over the domains.
3. To get the correct colors for the domains and phylas I put all colors inside a named vector.
4. The breaks for the fill scale are set such the domains come first which ensures that the domain is at the head of each legend
5. To set the order of the legends I use the order argument of guide_legend inside scale_fill_manual.

Finally, I was not sure what you meant by "... keep the legend titles 'Domain' and 'Phylum' in the legend". So I simply put the name of the domain in the legend title.

Domain_Colors <- c("orange", "blue", "grey", "green", "purple")
Phyla_Colors <- c("#FEEDDE", "#FDBE85", "#FD8D3C", "#E6550D", "#A63603", "#440154FF", 
                  "#443A83FF", "#31688EFF", "#21908CFF", "#35B779FF", "#8FD744FF", 
                  "#FDE725FF", "#4D4D4D", "#969696", "#C3C3C3", "#E6E6E6", "green1", 
                  "purple1")

Domain_Colors <- setNames(Domain_Colors, unique(df1$Domain))
Phyla_Colors <- setNames(Phyla_Colors, unique(df1$Phylum))

colors <- c(Domain_Colors, Phyla_Colors)
breaks_colors <- names(colors)

library(ggplot2)
library(ggnewscale)
library(purrr)
library(dplyr)

layers_help <- function(d, .domain, .order) {
  d <- filter(d, Domain == .domain)
  list(
    geom_rect(data = d, aes(fill = Domain, ymax = ymax, ymin = ymin, xmax = 2, xmin = 0)),
    geom_rect(data = d, aes(fill = Phylum, ymax = ymax, ymin = ymin, xmax = 4, xmin = 2)),
    scale_fill_manual(aesthetics = "fill", breaks = breaks_colors, values = colors, 
                      guide = guide_legend(title = .domain, order = as.integer(.order))),
    new_scale_fill()
  )
}

ggplot(df1) +
  imap(setNames(unique(df1$Domain), seq_along(unique(df1$Domain))), ~ layers_help(df1, .x, .y)) +
  coord_polar(theta = "y") +
  theme_void() +
  theme(legend.text=element_text(size = 4), legend.title=element_text(size = 6))

Another try. As the answer follows a different approach than my first I post it as a second answer. Moreover, even this second solution is not perfect. The major drawback is that it works only for monospaced fonts.

The general idea is to stick with your plotting code but fill up the data so that each domain contains an equal number of phyla or categories. Doing so makes it possible to arrange the phyla columnwise by domain. For the empty or auxilliary categories I set the fill color to transparent and use a labeller function which sets the labels to an empty string so that they do not appear in the legend.

However, the tricky part is to align the legend keys for the domain legend with the keys of the phyla legend. The issue is that the width of the legend is determined by the length of the labels and the chosen font. To make all labels the same length we can add additional spaces. However, the width of the labels will still vary depending on the chosen font (by default: "sans"). Unfortunately I'm not aware of an easy way to fix this, except for the case of a monospace font. In that case length is equivalent to width and simply adding additional spaces to the labels is enough to make sure that the items of the domain legend align with the items in the phyla legend.

library(ggplot2)
library(ggnewscale)
library(tidyr)
library(dplyr)
library(stringr)

d <- df1 %>% 
  group_by(Domain) %>% 
  mutate(rank = row_number(Phylum),
         max_len = max(nchar(Phylum))
         ) %>% 
  ungroup() %>% 
  complete(Domain, rank, fill = list(ymin = 0, ymax = 0, xmin = 0, xmax = 0, Phylum = "empty")) %>% 
  fill(max_len) %>% 
  mutate(Phylum = if_else(Phylum == "empty", paste0(Domain, rank), Phylum),
         Phylum = if_else(str_detect(Phylum, "\\d$"), str_pad(Phylum, max_len, "left"), str_pad(Phylum, max_len, "right")),
         Domain = str_pad(Domain, max_len, "right"))

is_empty <- str_detect(d$Phylum, "\\d$")
empty_cat <- d$Phylum[is_empty]
non_empty_cat <- d$Phylum[!is_empty]

Domain_Colors <- c("orange", "blue", "grey", "green", "purple")
Phyla_Colors <- c("#FEEDDE", "#FDBE85", "#FD8D3C", "#E6550D", "#A63603", "#440154FF", 
                  "#443A83FF", "#31688EFF", "#21908CFF", "#35B779FF", "#8FD744FF", 
                  "#FDE725FF", "#4D4D4D", "#969696", "#C3C3C3", "#E6E6E6", "green1", 
                  "purple1")

Domain_Colors <- setNames(Domain_Colors, unique(d$Domain))
Phyla_Colors <- setNames(Phyla_Colors, non_empty_cat)

colors <- c(Domain_Colors, Phyla_Colors)
colors <- setNames(colors, names(colors))
colors <- c(colors, setNames(rep("transparent", length(empty_cat)), empty_cat))

breaks_phylum <- d %>% 
  arrange(Domain, rank) %>%
  pull(Phylum)
breaks <- c(unique(d$Domain), breaks_phylum)

mylabels <- function(breaks) {
  if_else(str_detect(breaks, "\\d$"), "", breaks)
}

ggplot(d) +
  geom_rect(aes(fill = Domain, ymax = ymax, ymin = ymin, xmax = 2, xmin = 0)) +
  scale_fill_manual(aesthetics = "fill", breaks = breaks, values = colors, labels = mylabels, 
                    guide = guide_legend(title = "Domain", order = 1, nrow = 1)) +
  new_scale_fill() +
  geom_rect(aes(fill = Phylum, ymax = ymax, ymin = ymin, xmax = 2, xmin = 0)) +
  scale_fill_manual(aesthetics = "fill", breaks = breaks, values = colors, labels = mylabels,
                    guide = guide_legend(title = "Phylum", order = 0, nrow = 7)) +
  coord_polar(theta = "y") +
  theme_void(base_family = "mono") +
  theme(legend.text=element_text(size = 4), 
        legend.title=element_text(size = 6, vjust = 1), 
        legend.position = "bottom", 
        legend.box = "vertical",
        legend.box.just = "left")