In this article, we propose and analyze a fully coupled, nonlinear, and energy-stable virtual element method (VEM) for solving the coupled Poisson-Nernst-Planck (PNP) and Navier-Stokes (NS) equations modeling microfluidic and electrochemical systems (diffuse transport of charged species within incompressible fluids coupled through electrostatic forces). A mixed VEM is employed to discretize the NS equations whereas classical VEM in primal form is used to discretize the PNP equations. The stability, existence and uniqueness of solution of the associated VEM are proved by fixed point theory. Global mass conservation and electric energy decay of the scheme are also proved. Also, we obtain unconditionally optimal error estimates for both the electrostatic potential and ionic concentrations of PNP equations in the H 1 -norm, as well as for the velocity and pressure of NS equations in the H 1 -and L 2 -norms, respectively. Finally, several numerical experiments are presented to support the theoretical analysis of convergence and to illustrate the satisfactory performance of the method in simulating the onset of electrokinetic instabilities in ionic fluids, and studying how they are influenced by different values of ion concentration and applied voltage. These tests relate to applications in the desalination of water.Keywords Coupled Poisson-Nernst-Planck/Navier-Stokes equations • mixed virtual element method • optimal convergence • charged species transport • electrokinetic instability • water desalination.