A search for the production of three massive vector bosons, VVZ(V=W,Z), in proton–proton collisions at s=13 TeV is performed using data with an integrated luminosity of 140 fb−1 recorded by the ATLAS detector at the Large Hadron Collider. Events produced in the leptonic final states WWZ→ℓνℓνℓℓ (ℓ=e,μ), WZZ→ℓνℓℓℓℓ, ZZZ→ℓℓℓℓℓℓ, and the semileptonic final states WWZ→qqℓνℓℓ and WZZ→ℓνqqℓℓ, are analysed. The measured cross section for the pp→VVZ process is 660−90+93(stat.)−81+88(syst.) fb, and the observed (expected) significance is 6.4 (4.7) standard deviations, representing the observation of VVZ production. In addition, the measured cross section for the pp→WWZ process is 442±94(stat.)−52+60(syst.) fb, and the observed (expected) significance is 4.4 (3.6) standard deviations, representing evidence of WWZ production. The measured cross sections are consistent with the Standard Model predictions. Constraints on physics beyond the Standard Model are also derived in the effective field theory framework by setting limits on Wilson coefficients for dimension-8 operators describing anomalous quartic gauge boson couplings.