Summary: Let \(X_n = \{x^j\}_{j=1}^n\) be a set of \(n\) points in the \(d\)-cube \({\mathbb {I}}^d:=[0,1]^d\), and \(\Phi _n = \{\varphi _j\}_{j =1}^n\) a family of \(n\) functions on \({\mathbb {I}}^d\). We consider the approximate recovery of functions \(f\) on \({{\mathbb {I}}}^d\) from the sampled values \(f(x^1),\dotsc, f(x^n)\), by the linear sampling algorithm \( L_n(X_n,\Phi _n,f) := \sum _{j=1}^n f(x^j)\varphi _j. \) The error of sampling recovery is measured in the norm of the space \(L_q({\mathbb {I}}^d)\)-norm or the energy quasi-norm of the isotropic Sobolev space \(W^\gamma _q({\mathbb {I}}^d)\) for \(1 < q < \infty \) and \(\gamma > 0\). Functions \(f\) to be recovered are from the unit ball in Besov-type spaces of an anisotropic smoothness, in particular, spaces \(B^{\alpha ,\beta }_{p,\theta }\) of a “hybrid” of mixed smoothness \(\alpha > 0\) and isotropic smoothness \(\beta \in {\mathbb {R}}\), and spaces \(B^a_{p,\theta }\) of a nonuniform mixed smoothness \(a \in {\mathbb {R}}^d_+\). We constructed asymptotically optimal linear sampling algorithms \(L_n(X_n^\ast,\Phi _n^\ast,\cdot )\) on special sparse grids \(X_n^\ast\) and a family \(\Phi _n^\ast\) of linear combinations of integer or half integer translated dilations of tensor products of B-splines. We computed the asymptotic order of the error of the optimal recovery. This construction is based on B-spline quasi-interpolation representations of functions in \(B^{\alpha ,\beta }_{p,\theta }\) and \(B^a_{p,\theta }\). As consequences, we obtained the asymptotic order of optimal cubature formulas for numerical integration of functions from the unit ball of these Besov-type spaces.