We present variational Monte Carlo calculations of the neutron matter equation of state using chiral nuclear forces. The ground-state wave function of neutron matter, containing nonperturbative many-body correlations, is obtained from auxiliary-field quantum Monte Carlo simulations of up to about 340 neutrons interacting on a 103 discretized lattice. The evolution Hamiltonian is chosen to be attractive and spin independent in order to avoid the fermion sign problem and is constructed to best reproduce broad features of the chiral nuclear force. This is facilitated by choosing a lattice spacing of 1.5 fm, corresponding to a momentum-space cutoff of Lambda = 414 MeV/c, a resolution scale at which strongly repulsive features of nuclear two-body forces are suppressed. Differences between the evolution potential and the full chiral nuclear interaction (Entem and Machleidt Lambda = 414 MeV [L. Coraggio et al., Phys. Rev. C 87, 014322 (2013)]) are then treated perturbatively. Our results for the equation of state are compared to previous quantum Monte Carlo simulations that employed chiral two-body forces at next-to-next-to-leading order (N2LO). In addition, we include the effects of three-body forces at N2LO, which provide important repulsion at densities higher than 0.02 fm(-3) as well as two-body forces at next-to-next-to-next-to-leading order.