We present variational Monte Carlo calculations of the neutron matter equation of state using chiral nuclear forces. The ground-state wavefunction of neutron matter, containing non-perturbative many-body correlations, is obtained from auxiliary-field quantum Monte Carlo simulations of up to about 340 neutrons interacting on a 10^3 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) are then treated perturbatively. Our results for the equation of state are compared to previous quantum Monte Carlo simulations which 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 N3LO.