I came up with a rather unsophisticated solution, where I have padded the initial ds with Nan values for latitude -80 to -90, and then clip the ds for the min/max in ESPG:3857 before the reprojection. When I generate an image from this, and put it on top of a 3857 map, it seems fairly correct.

However, I would be very happy if someone would comment on my approach.

ds = ds.sel(latitude=slice(85.05112878, -85.05112878))

ds.rio.bounds()
(-180.0416666672555, -85.04181929195629, 179.95833842018519, 85.04166670406566)
ds_mercator.rio.bounds()
(-20037508.232311152, -20026561.83042367, 20031128.729270592, 20025310.01109047)

Reprojected after new padding and clipping:

most likely because of lon=180 crossing over, instead of spanning [-180, 180] your image is shifted by half a pixel left and hence spans from just under lon=-180 to just under lon=180. X coordinate of the data is supposed to contain location of the center of the pixel, not of it's left edge, same for Y. Adjust your coords when constructing ds, such that bounds are reported in the correct range for espg4326