The gradient calculation would also yield the same results as seen here:
output = torch.randn(100, 100, requires_grad=True)
labels = torch.randn(100, 100)

torch.mean(torch.square(output - labels)).backward()
g0 = output.grad.clone()
print(output.grad)
output.grad = None

torch.mean((output - labels)**2).backward()
g1 = output.grad.clone()
print(output.grad)
output.grad = None

((output - labels)**2).mean().backward()
g2 = output.grad.clone()
print(output.grad)
output.grad = None

criterion = nn.MSELoss()
loss = criterion(output, labels)
loss.backward()
g3 = output.grad.clone()
print(output.grad)
output.grad = None

loss = F.mse_loss(output, labels)
loss.backward()
g4 = output.grad.clone()
print(output.grad)
output.grad = None

print(torch.allclose(g0, g1))
> True
print(torch.allclose(g0, g2))
> True
print(torch.allclose(g0, g3))
> True
print(torch.allclose(g0, g4))
> True

How reproducible is the difference, how many runs did you use for each approach, and what is the mean of the final loss as well as its stddev?