Creating a mask for natural earth R giving unexpected results

I am wanting to create a mask to color all parts of the image that are NOT within the specified land areas. I can get this to work successfully for two countries (e.g Burkina Faso and Nigeria). However, the same code does not work when I try a longer list of 9 countries in the same area. The error I am receiving is : Error in [[<-.data.frame(tmp, attr(x, "sf_column"), value = list( : replacement has 2 rows, data has 1

Being relatively new to GIS, I am not quite sure what the issue is or how to resolve it. I would really appreciate some help if people could offer some insight. The method I used was based on the accepted answer from this question

Relevant code is below. Thankyou in advance for anyone who is able to help.

This code works perfectly fine

library(sf)
library(dplyr)
library(tidyverse)
library(rnaturalearth)
library(rnaturalearthdata)

# https://stackoverflow.com/questions/49266736/clip-spatial-polygon-by-world-map-in-r
x_coord = c(-20,-20,27,27)
y_coord = c(-1,28,28,-1)

# Create rectalgular polygon for a mask (slightly bigger than the image)
polygon <-  cbind(x_coord, y_coord) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_sfc(crs = 4326, check_ring_dir = TRUE) %>%
  st_sf() %>%
  st_wrap_dateline()

land = ne_countries(scale = "medium", 
             returnclass = "sf") %>% 
     filter(admin %in% c("Nigeria", "Burkina Faso"))
    #filter(admin %in% c("Niger", "Burkina Faso", "Nigeria", "Mali", "Chad", "Mauritania",
    #              "Gambia", "Senegal", "Guinea", "Cameroon"))


# Get the difference between the bounding box and the area I want to plot
# This works fine when using only two countries 
land = st_union(land)
polygon_land_diff = st_difference(polygon, land)

# Plot the data 
plot(st_geometry(land))
plot(st_geometry(polygon), add=TRUE)
plot(st_geometry(polygon_land_diff),add=TRUE, col="green")

This code produces an error and I am not sure how to resolve it. The only difference is the number of countries selected

x_coord = c(-20,-20,27,27)
y_coord = c(-1,28,28,-1)

# Create rectalgular polygon for a mask (slightly bigger than the image)
polygon <-  cbind(x_coord, y_coord) %>%
  st_linestring() %>%
  st_cast("POLYGON") %>%
  st_sfc(crs = 4326, check_ring_dir = TRUE) %>%
  st_sf() %>%
  st_wrap_dateline()

land = ne_countries(scale = "medium", 
             returnclass = "sf") %>% 
    filter(admin %in% c("Niger", "Burkina Faso", "Nigeria", "Mali", "Chad", "Mauritania",
                  "Gambia", "Senegal", "Guinea", "Cameroon"))


# get the difference between the bounding box and the area I want to plot
# not quite working. The rror occurs when running  `polygon_land_diff = st_difference(polygon, land)`

land = st_union(land)
polygon_land_diff = st_difference(polygon, land)

# Plot the data 
plot(st_geometry(land))
plot(st_geometry(polygon), add=TRUE)
plot(st_geometry(polygon_land_diff),add=TRUE, col="green")