The Virtual Element Method (VEM) is a novel family of numerical methods for approximating partial differential equations on very general polygonal or polyhedral computational grids. This work aims to propose a Balancing Domain Decomposition by Constraints (BDDC) preconditioner that allows using the conjugate gradient method to compute the solution of the saddle-point linear systems arising from the VEM discretization of the three-dimensional Stokes equations. We prove the scalability and quasi-optimality of the algorithm and confirm the theoretical findings with parallel computations. Numerical results with adaptively generated coarse spaces confirm the method's robustness in the presence of large jumps in the viscosity and with high-order VEM discretizations.