The means are equal
Suppose $f$ is analytic so that we can give it a Taylor series that works everywhere such that $f(x) = \sum a_n x^n$, and then let us let this be bounded too. To show that the expectations are the same consider taking the expectation of $U_t$
\begin{equation}
\mathbb{E}(U_t) = \mathbb{E}\left(\int_0^t f(W_s) \,\mathrm{d}s\right)\\
\end{equation}
as $f$ if bounded we can use Fubini's theorem to move the expectation inside the integral, giving
\begin{align}
\mathbb{E}(U_t) & = \int_0^t \mathbb{E}\left(f(W_s)\right) \,\mathrm{d}s\\
& = \int_0^t \sum_n a_n \mathbb{E}\left(W_s^n\right) \,\mathrm{d}s \\
& = \int_0^t \sum_n a_n s^{n/2} M_n \,\mathrm{d}s
\end{align}
where $M_n$ denotes the $n$-th moment of a standard Gaussian distribution. Identically we could do this for the expectation of $V_t$ and see that the two are equal, and hence that $U_t$ has the same expectation of $V_t$.
The variances are not equal
For the variance of $U_t$ we need to evaluate the expected value of
\begin{align}
U_t^2 = \left(\int_0^t f(W_s) \,\mathrm{d}s \right)^2
&= \int_0^t \int_0^t f(W_s) f(W_u) \,\mathrm{d}u \,\mathrm{d}s \\
&= 2 \int_0^t \int_0^s f(W_s) f(W_u) \,\mathrm{d}u \,\mathrm{d}s.
\end{align}
Taking the expectation and putting this inside the integral again gives
\begin{equation}
\mathbb{E}\left(U_t^2\right) =
2 \int_0^t \int_0^s \sum_{n, m, k} a_n a_m \binom{n}{k} \mathbb{E}(W_s^n W_u^m) \,\mathrm{d}u \,\mathrm{d}s.
\end{equation}
Within the integrand $u \leq s$ and so expand $W_s = W_u + (W_s - W_u)$ where as $W$ is a Weiner process the increment $W_s - W_u$ is independent to the value of $W_u$ and has variance $s - u$. This then simplifies to give
\begin{equation}
\mathbb{E}\left(U_t^2\right) =
2 \int_0^t \int_0^s \sum_{n, m, k} a_n a_m \binom{n}{k} u^{\frac{n + m}{2}}(s - u)^{\frac{n-k}{2}}M_{n+m}M_{n-k} \,\mathrm{d}u \,\mathrm{d}s.
\end{equation}
An identical process can be applied for $V_t$ but instead of expanding $W_s$ we just use $X_t = \sqrt{t}Z$ and obtain
\begin{equation}
\mathbb{E}\left(V_t^2\right) =
2 \int_0^t \int_0^s \sum_{n, m, k} a_n a_m \binom{n}{k} s^{\frac{n}{2}}u^{\frac{m}{2}}M_{n+m} \,\mathrm{d}u \,\mathrm{d}s.
\end{equation}
The second expression is more easily integrated out than the first, but both are doable and I think they evaluate to different quantities.
A numeric demonstration
Take $f \equiv \sin$ and you can simulate this, as the example python code does, the key result is then:
Variable = U, Mean = 0.00206326180593 +/- 0.00228420660865, Standard deviation = 0.228420660865
Variable = V, Mean = -0.0177668651186 +/- 0.0380135632954, Standard deviation = 3.80135632954
We see this confirms that the means are the same but variances differ.

Python code
import numpy as np
from scipy.stats import norm
f = np.sin
M = 10000
T = 2*np.pi
dt = T/M
N = 10000
mean_u, mean_v = 0, 0
var_u, var_v = 0, 0
for _ in range(N):
w = np.cumsum(np.concatenate([[0], norm.rvs(scale=dt, size=M)]))
z = norm.rvs()
t = np.cumsum([0] + [dt]*M)
x = np.sqrt(t) * z
u = np.sum(f(w))*dt
v = np.sum(f(x))*dt
mean_u += u
mean_v += v
var_u += u ** 2
var_v += v ** 2
mean_u /= N
mean_v /= N
var_u /= N
var_v /= N
var_u -= mean_u ** 2
var_v -= mean_v ** 2
std_u = np.sqrt(var_u)
std_v = np.sqrt(var_v)
std_error_u = np.sqrt(var_u / N)
std_error_v = np.sqrt(var_v / N)
for n, m, s, se in [["U", mean_u, std_u, std_error_u], ["V", mean_v, std_v, std_error_v]]:
print("Variable = {},\tMean = {} +/- {},\tStandard deviation = {} ".format(n, m, se, s))