In this paper we consider a coupled bulk-surface PDE in two space dimensions. The model consists of a PDE in the bulk that is coupled to another PDE on the surface through general nonlinear boundary conditions. For such a system we propose a novel method, based on coupling a virtual element method (Beirão Da Veiga et al. in Math Models Methods Appl Sci 23(01):199–214, 2013. https://doi.org/10.1051/m2an/2013138) in the bulk domain to a surface finite element method (Dziuk and Elliott in Acta Numer 22:289–396, 2013. https://doi.org/10.1017/s0962492913000056) on the surface. The proposed method, which we coin the bulk-surface virtual element method includes, as a special case, the bulk-surface finite element method (BSFEM) on triangular meshes (Madzvamuse and Chung in Finite Elem Anal Des 108:9–21, 2016. https://doi.org/10.1016/j.finel.2015.09.002). The method exhibits second-order convergence in space, provided the exact solution is $$H^{2+1/4}$$
H
2
+
1
/
4
in the bulk and $$H^2$$
H
2
on the surface, where the additional $$\frac{1}{4}$$
1
4
is required only in the simultaneous presence of surface curvature and non-triangular elements. Two novel techniques introduced in our analysis are (i) an $$L^2$$
L
2
-preserving inverse trace operator for the analysis of boundary conditions and (ii) the Sobolev extension as a replacement of the lifting operator (Elliott and Ranner in IMA J Num Anal 33(2):377–402, 2013. https://doi.org/10.1093/imanum/drs022) for sufficiently smooth exact solutions. The generality of the polygonal mesh can be exploited to optimize the computational time of matrix assembly. The method takes an optimised matrix-vector form that also simplifies the known special case of BSFEM on triangular meshes (Madzvamuse and Chung 2016). Three numerical examples illustrate our findings.