We establish that the Riesz transforms of all orders corresponding to the Grusin operator \(H_N=-\nabla _{x}^2-|x|^{2N}\,\nabla _{y}^2\), and the first-order operators \((\nabla _{x},x^\nu \,\nabla _{y})\) where \(x\in \mathbf{R}^n\), \(y\in \mathbf{R}^m\), \(N\in \mathbf{N}_+\), and \(\nu \in \{1,\ldots ,n\}^N\), are bounded on \(L_p(\mathbf{R}^{n+m})\) for all \(p\in \langle 1,\infty \rangle \) and are also weak-type (1, 1). Moreover, the transforms of order less than or equal to \(N+1\) corresponding to \(H_N\) and the operators \((\nabla _{x}, |x|^N\nabla _{y})\) are bounded on \(L_p(\mathbf{R}^{n+m})\) for all \(p\in \langle 1,\infty \rangle \). But if N is odd all transforms of order \(N+2\) are bounded if and only if \(p\in \langle 1,n\rangle \). The proofs are based on the observation that the \((\nabla _{x},x^\nu \,\nabla _{y})\) generate a finite-dimensional nilpotent Lie algebra, the corresponding connected, simply connected, nilpotent Lie group is isometrically represented on the spaces \(L_p(\mathbf{R}^{n+m})\) and \(H_N\) is the corresponding sublaplacian.