  In the present paper the axial quasi-normal modes of neutron stars in $f(R)$
gravity are examined using a large set of equations of state. The numerical
calculations are made using two different approaches -- performing time
evolution of the perturbation equations and solving the time-independent
representation of the equations as a boundary value problem. According to the
results the mode frequencies and the damping times decrease with the increase
of the free parameter of the theory in comparison to the pure general
relativistic case. While the frequencies deviate significantly from Einstein's
theory for all realistic neutron star masses (say above $1M_\odot$), the
damping times reach non-negligible differences only for the more massive
models. We have constructed as well universal (equation of state independent)
gravitational wave asteroseismology relations involving the frequencies and the
damping times. It turns out that the equation of state independence is
preserved using the same normalization as in pure general relativity and the
qualitative differences of the phenomenological relations with respect to
Einstein's theory of gravity can be large for large values of the free
parameter in $f(R)$ gravity.
