Evidence is presented for the electroweak (EW) production of two jets (jj) in association with two Z bosons and constraints on anomalous quartic gauge couplings are set. The analysis is based on a data sample of proton-proton collisions at s=13TeV collected with the CMS detector in 2016–2018, and corresponding to an integrated luminosity of 137fb−1. The search is performed in the fully leptonic final state ZZ→ℓℓℓ′ℓ′, where ℓ,ℓ′=e,μ. The EW production of two jets in association with two Z bosons is measured with an observed (expected) significance of 4.0 (3.5) standard deviations. The cross sections for the EW production are measured in three fiducial volumes and the result is σEW(pp→ZZjj→ℓℓℓ′ℓ′jj)=0.33−0.10+0.11(stat)−0.03+0.04(syst)fb in the most inclusive volume, in agreement with the standard model prediction of 0.275±0.021fb. Measurements of total cross sections for jj production in association with two Z bosons are also reported. Limits on anomalous quartic gauge couplings are derived in terms of the effective field theory operators T0, T1, T2, T8, and T9.