You can consider getting the NUTS regions via {giscoR} package - it serves formally valid NUTS regions.

nuts <- giscoR::gisco_get_nuts(year = '2016')
all(sf::st_is_valid(nuts))
[1] TRUE

You can consider getting the NUTS regions via {giscoR} package - it serves formally valid NUTS regions.

nuts <- giscoR::gisco_get_nuts(year = '2016')
all(sf::st_is_valid(nuts))
[1] TRUE

The example is now:

library(sf)
library(giscoR)
NUTS_WGS84 = gisco_get_nuts(year = '2016')
y_coord <- c(46.976, 46.948)
x_coord <- c(7.483, 7.45)
coordinates_as_df <- data.frame(x_coord, y_coord)
coordinates_WGS84 <- st_as_sf(x = coordinates_as_df, 
                                  coords = c("x_coord", "y_coord"), 
                                  crs = 4326)
overlay <- st_join(coordinates_WGS84, NUTS_WGS84, join=st_intersects)
print(overlay[overlay$LEVL_CODE == 3, ])

And the output:

Simple feature collection with 2 features and 9 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 7.45 ymin: 46.948 xmax: 7.483 ymax: 46.976
Geodetic CRS:  WGS 84
  NAME_LATN MOUNT_TYPE CNTR_CODE URBN_TYPE COAST_TYPE   FID NUTS_NAME NUTS_ID LEVL_CODE             geometry
4      Bern          3        CH         2          0 CH021      Bern   CH021         3 POINT (7.483 46.976)
8      Bern          3        CH         2          0 CH021      Bern   CH021         3  POINT (7.45 46.948)