The shape evolution and shape coexistence phenomena in neutron-rich nuclei at N ≈ 60, including Kr, Sr, Zr, and Mo isotopes, are studied in the covariant density functional theory (DFT) with the new parameter set PC-PK1. Pairing correlations are treated using the BCS approximation with a separable pairing force. Sharp rising in the charge radii of Sr and Zr isotopes at N = 60 is observed and shown to be related to the rapid changing in nuclear shapes. The shape evolution is moderate in neighboring Kr and Mo isotopes. Similar as the results of previous Hartree-Fock-Bogogliubov (HFB) calculations with the Gogny force, triaxiality is observed in Mo isotopes and shown to be essential to reproduce quantitatively the corresponding charge radii. In addition, the coexistence of prolate and oblate shapes is found in both 98 Sr and 100 Zr. The observed oblate and prolate minima are related to the low single-particle energy level density around the Fermi surfaces of neutron and proton respectively. Furthermore, the 5-dimensional (5D) col-