Bulk-surface partial differential equations (BS-PDEs) are prevalent in many applications such as cellular, developmental and plant biology as well as in engineering and material sciences. Novel numerical methods for BS-PDEs in three space dimensions (3D) are sparse. In this work, we present a bulk-surface virtual element method (BS-VEM) for bulk-surface reaction-diffusion systems, a form of semilinear parabolic BS-PDEs in 3D. Unlike previous studies in two space dimensions (2D), the 3D bulk is approximated with general polyhedra, whose outer faces constitute a flat polygonal approximation of the surface. For this reason, the method is restricted to the lowest order case where the geometric error is not dominant. The BS-VEM guarantees all the advantages of polyhedral methods such as easy mesh generation and fast matrix assembly on general geometries. Such advantages are much more relevant than in 2D. Despite allowing for general polyhedra, general nonlinear reaction kinetics and general surface curvature, the method only relies on nodal values without needing additional evaluations usually associated with the quadrature of general reaction kinetics. This latter is particularly costly in 3D. The BS-VEM as implemented in this study retains optimal convergence of second order in space.