You can consider getting the NUTS regions via [{giscoR}](https://dieghernan.github.io/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)