  The linear collisionless plasma response to a zonal density perturbation in
quasisymmetric stellarators is studied, including the initial
geodesic-acoustic-mode oscillations and the final Rosenbluth--Hinton residual
flow. While the geodesic-acoustic-mode oscillations in quasiaxisymmetric
configurations are similar to tokamaks, they become almost nonexistent in
quasi-helically symmetric configurations when the effective safety factor is
small. Compared with concentric circular tokamaks, the Rosenbluth--Hinton
residual is also found to be multiplied by a geometric factor $\mathcal{C}$
that arises from the flux-surface averaged classical polarization. Using the
near-axis-expansion framework, we derive an analytic expression for
$\mathcal{C}$, which varies significantly among different configurations. These
analytic results are compared with numerical simulation results from the global
gyrokinetic particle-in-cell code GTC, and good agreement in the
Rosenbluth--Hinton residual level is achieved only when the quasisymmetry error
is small enough. Since zonal flows can be important for regulating turbulent
transport, these results suggest possible relation between the transport level
and the stellarator geometric parameters via nonlinear interactions with zonal
flows.
