While (as I said in a comment above) no examples of (continuously) flexible smooth (by *smooth* I mean of class $C^2$ or above) surfaces are known, there are known examples of closed smooth surfaces that are *infinitesimally flexible*, (i.e. they admit infinitesimal isometric deformations, or infinitesimal bendings, defined in the next paragraph). Thus one can ask whether the volumes of these surfaces are preserved under infinitesimal bendings. I found some literature on this question for surfaces of revolution.

[Let $\Sigma$ be some closed surface and $f:\Sigma\rightarrow\mathbb{R}^3$ an embedding into $\mathbb{R}^3$. An infinitesimal bending of $(\Sigma,f)$ is a vector field $u$ so that the Riemannian metric of the embedding $f+\epsilon u$ induced by that of $\mathbb{R}^3$ differs from that of $f$ by $o(\epsilon)$.]

In Sabitov's survey article "Local theory of bendings" in the book "Geometry III: Theory of Surfaces" (1992), he writes about the bellows conjecture in section 9.8 ("Conjecture on the Invariance of the Volume of a Bendable Polyhedron"). After describing the fact that the bellows conjecture (for polyhedra) was then still unproven (note that he was the first to give a proof, in 1996), he writes:

On the contrary, for smooth surfaces there is an encouraging result: for all non-rigid surfaces of revolution their volume remains stationary under an infinitesimal bending: see V.A. Aleksandrov (1989)[, researchgate link]. Although up to now there are no examples of regular closed bendable surfaces, nevertheless information on the behaviour of their volume under a priori assumed bendings can be useful.

The paper of Alexandrov is titled "Sabitov's conjecture that volume is stationary under infinitesimal bending of a surface" and the relevant theorem is:

THEOREM 2. Let the meridian of a connected closed surface of rotation be a
$C^1$-smooth curve and not contain a segment perpendicular to the axis of rotation. Then the flux through the given surface of any $C^1$-smooth bending field is equal to zero.

The first condition rules out surfaces of revolution which contain planar pieces, as these do admit infinitesimal bendings which change the volume to first order due to the large number of infinitesimal bendings of a planar domain. Note also that the flux condition is equivalent to the stationarity of the volume to first order.

I also found a 1999 paper "On variation of the volume under infinitesimal bending of a closed rotational surface" by Lj. S. Velimirović which treats the piecewise smooth case and gives an example calculation of the volume variation of "Belov's rotational toroid".

In a comment, the OP asks about the $C^1$ embeddings of surfaces. The answer in this case is that bellows **do exist**: given a closed Riemannian surface $\Sigma$, there exist (many) continuous families of $C^1$ isometric embeddings of $\Sigma$ that do not preserve volume. This is a consequence of "the parametric $h$-principle for isometric $C^1$ embeddings", due to Nash-Kuiper and Gromov; see chapter 21 of the book "Introduction to the h-principle" by Eliashberg and Mishachev for a proof of the non-parametric version and an assurance that it can be adapted to families of embeddings.

Here's the idea. Nash and Kuiper originally proved that any embedding of a Riemannian surface into $\mathbb{R}^3$ which is short can be $C^0$-approximated by an isometric $C^1$ embedding (actually they proved it in general dimensions). Gromov reformulated this in terms of the h-principle and showed as a consequence that this result extended to continuous families of embeddings. Suppose we have a 2-dimensional submanifold $\Sigma\subset\mathbb{R}^3$, then we can construct a 1-parameter family of short embeddings of $\Sigma$ by uniformly shrinking the initial embedding. From the parametric h-principle, we get a 1-parameter family of isometric $C^1$ embeddings that approximates this shrinking family, and it's then not hard to see that the volume of these approximating isometric embeddings is decreasing and not constant.

You may enjoy the work of the Hévéa project, a collaboration that has been working to make constructions of $C^1$ isometric embeddings explicit enough that one can compute images. See this recent paper of Bartzos, Borrelli, Denis, Lazarus, Rohmer and Thibert which gives full details on how to isometrically embed the unit sphere into arbitrarily small balls:

I suspect that their techniques will allow for the explicit construction of families of isometric embeddings, but it may take quite a bit of work to do this.