Let \(\nu = (\nu _1, \ldots , \nu _n) \in (-1/2, \infty )^n\) , with \(n \ge 1\) , and let \(\Delta _\nu \) be the multivariate Bessel operator defined by \( \Delta _{\nu } = -\sum _{j=1}^n\left( \frac{\partial ^2}{\partial x_j^2} - \frac{\nu _j^2 - 1/4}{x_j^2} \right) . \) In this paper, we develop the theory of Hardy spaces and BMO-type spaces associated with the Bessel operator \(\Delta _\nu \) . We then study the higher-order Riesz transforms associated with \(\Delta _\nu \) . First, we show that these transforms are Calderón-Zygmund operators. We further prove that they are bounded on the Hardy spaces and BMO-type spaces associated with \(\Delta _\nu \) .