A mixed virtual element method (mixed-VEM) for a pseudostress-velocity formulation of the two-dimensional Navier–Stokes equations with Dirichlet boundary conditions is proposed and analyzed in this work. More precisely, we employ a dual-mixed approach based on the introduction of a nonlinear pseudostress linking the usual linear one for the Stokes equations and the convective term. In this way, the aforementioned new tensor together with the velocity constitute the only unknowns of the problem, whereas the pressure is computed via a postprocessing formula. In addition, the resulting continuous scheme is augmented with Galerkin type terms arising from the constitutive and equilibrium equations, and the Dirichlet boundary condition, all them multiplied by suitable stabilization parameters, so that the Banach fixed-point and Lax–Milgram theorems are applied to conclude the well-posedness of the continuous and discrete formulations. Next, we describe the main VEM ingredients that are required for our discrete analysis, which, besides projectors commonly utilized for related models, include, as the main novelty, the simultaneous use of virtual element subspaces for [Formula: see text] and [Formula: see text] in order to approximate the velocity and the pseudostress, respectively. Then, the discrete bilinear and trilinear forms involved, their main properties and the associated mixed virtual scheme are defined, and the corresponding solvability analysis is performed using again appropriate fixed-point arguments. Moreover, Strang-type estimates are applied to derive the a priori error estimates for the two components of the virtual element solution as well as for the fully computable projections of them and the postprocessed pressure. As a consequence, the corresponding rates of convergence are also established. Finally, we follow the same approach employed in previous works by some of the authors and introduce an element-by-element postprocessing formula for the fully computable pseudostress, thus yielding an optimally convergent approximation of this unknown with respect to the broken [Formula: see text]-norm.