In this paper, we focus on model reduction of large-scale bilinear systems. The main contributions are threefold. First, we introduce a new framework for interpolatory model reduction of bilinear systems. In contrast to the existing methods where interpolation is forced on some of the leading subsystem transfer functions, the new framework shows how to enforce multipoint interpolation of the underlying Volterra series. Then, we show that the first-order conditions for optimal H2 model reduction of bilinear systems require multivariate Hermite interpolation in terms of the new Volterra series interpolation framework; and thus we extend the interpolation-based first-order necessary conditions for H2 optimality of LTI systems to the bilinear case. Finally, we show that multipoint interpolation on the truncated Volterra series representation of a bilinear system leads to an asymptotically optimal approach to H2 optimal model reduction, leading to an efficient model reduction algorithm. Several numerical examples illustrate the effectiveness of the proposed approach.