You can verify that the gradients and Hessian agree by direct calculation. Observe that both $\nabla f[\vec{\mathbf x}_0]$ and ${\mathbf H}_f[\vec{\mathbf x}_0]$ are constants, so the second term of $Q$ is a linear function of $\vec{\mathbf x}-\vec{\mathbf x}_0$ and the third term a quadratic form in $\vec{\mathbf x}-\vec{\mathbf x}_0$. (I use square brackets to indicate evaluation at a point in an attempt at disambiguation in the following.)
Analogous to the single-variable differentiation rule $(cx)'=c$ we have $$\nabla(\vec{\mathbf v}^T\vec{\mathbf x}) = \nabla\left(\sum_i v_ix_i\right) = \vec{\mathbf v},$$ therefore $$\nabla(\nabla f[\vec{\mathbf x}_0]^T(\vec{\mathbf x}-\vec{\mathbf x}_0)) = \nabla(\nabla f[\vec{\mathbf x}_0]^T\vec{\mathbf x})-\nabla(\nabla f[\vec{\mathbf x}_0]^T\vec{\mathbf x}_0) = \nabla f[\vec{\mathbf x}_0].$$ More generally, the differential of a linear function $A\mathbf x$ of $\mathbf x$ is $A$—the best linear approximation to a linear map is the map itself.
For the quadratic term, we can use the product rule $\nabla(u^Tv)=u^T\nabla v+(\nabla u)^Tv$: $$\begin{align} \nabla\left(\frac12(\vec{\mathbf x}-\vec{\mathbf x}_0)^T{\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0)\right) &= \frac12 (\vec{\mathbf x}-\vec{\mathbf x}_0)^T \nabla({\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0)) + \frac12(\nabla(\vec{\mathbf x}-\vec{\mathbf x}_0))^T{\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0) \\
&= \frac12 (\vec{\mathbf x}-\vec{\mathbf x}_0)^T {\mathbf H}_f[\vec{\mathbf x}_0] + \frac12 {\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0) \\
&= {\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0) \end{align}$$ since ${\mathbf H}_f$ is symmetric. This vanishes at $\vec{\mathbf x}_0$, therefore $\nabla Q[\vec{\mathbf x}_0] = \nabla f[\vec{\mathbf x}_0]$. (We could’ve gotten to this without actually differentiating by noting that the result must be a linear combination of $x_i-x_{0i}$ terms, all of which vanish.)
For ${\mathbf H}_Q[\vec{\mathbf x}_0]$ we can proceed in a similar fashion by differentiating $\nabla Q = \nabla f[\vec{\mathbf x}_0]+{\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0)$ directly, producing ${\mathbf H}_f[\vec{\mathbf x}_0]$ per the earlier observation regarding the differential of a linear function, or note that the Hessian of a constant or linear function vanishes, while the Hessian of a quadratic form is 2 times the (constant) matrix of the form (analogously to the single-variable $(ax^2)''=2a$), so that the only term that survives when computing the Hessian directly from $Q$ is ${\mathbf H}_f[\vec{\mathbf x}_0]$.