**Edit**: Now I have found something that appears to answer my own question. It is section 2 in the paper "On Submanifolds and Quotients of Poisson and Jacobi Manifolds" by Ch.-M. Marle. (There, he transfer all on the cotangent bundle via the Legendre map and uses a different, but really equivalent, construction of the nonholonomic bracket.) Thank you.

In "On The Hamiltonian Formulation of Nonholonomic Mechanical Systems" of Van der Schaft and Maschke (this is a link), a mechanical system subject to linear constraints on the velocities is provided with a bracket $\lbrace\cdot,\cdot\rbrace_{nh}$, which is almost Poisson because it lacks only the Jacobi identity.

Notwithstanding this loss, it equally describes the time evolution of the system according with d'Alembert's principle, that is $\frac{d}{dt}f=\lbrace f,E\rbrace_{nh}$, where $E$ is the total energy of the system.

My interest was attracted above all by the following result:

$\lbrace\cdot,\cdot\rbrace_{nh}$ measures the holonomic character of the linear constraint, that is the bracket satisfies the Jacobi identity iff the linear constraint is an integrable distribution.

I am writing this question for the following reason:

Being their presentation deeply depending on coordinates, I was in search for a more geometric, invariant demonstration of this last result.

In order to convey more information, below I briefly sketched the geometric context as far as I have understood until now. Thank you.

Let us impose on mechanical system having configuration space $Q$ and Lagrangian $L$ a constraint linear on the velocities represented by $D\subset TQ$, a tangent distribution on $Q$.

Let us introduce $J:TTQ\to TTQ$ the endomorphism of $\tau_{TQ}:TTQ\to TQ$ defined locally by $$J(\frac{\partial}{\partial x_i})=\frac{\partial}{\partial v_i},\ J(\frac{\partial}{\partial v_i})=0$$ where (x,v) are the standard coordinates on $TQ$ associated to local coordinates $x$ on $Q$. Below $J^\ast$ will denote the transpose of $J$.

Let $\omega$ be the symplectic form on $TQ$ given by the pull-back through $\mathbb{F}L$ of the canonical symplectic form on $T^\ast Q$. Here $\mathbb{F}L:TQ\to T^|ast Q$ is the Legendre trasformation associated to $L$ which locally acts as $(x,v)\mapsto (x,\frac{\partial L}{\partial v})$. Below $\flat$ and $\sharp$ will denote the musical isomorphisms corresponding to $\omega$.

The Chetaev bundle $F$ of the reaction forces is the tangent distribution on $TQ$ such that $\flat(F)=J^\ast((TD)^0)$. Here $(TD)^0$ denotes the annihilator of $TD$.

If the Lagrangian is natural then $H:=F^\perp\cap TD$ is a symplectic distribution on $(TQ,\omega)$ , and we can define an endomorphism $P:T_DTQ\to T_DTQ$ by projecting fiberwise $T_DTQ=H^\perp\oplus H$ on $H$ along $H^\perp$.

For any pair $f,g$ of smooth functions on $D$, their bracket $\lbrace f,g\rbrace_{nh}$ is defined as $\omega(P\circ\sharp(d\tilde{f}),P\circ\sharp(d\tilde{g}))$ where $\tilde{f}$ and $\tilde{g}$ are arbitrary smooth extensions on $TQ$ of $f$ and $g$ respectively.

Such a bracket is skew-simmetric, bilinear and satisfies the Leibniz rule in each one of its argument, but couldn't be a Poisson bracket on $D$ because of the lack of Jacobi identity.