With the motivation of exploring new high-strength ceramics, ab initio evolutionary simulations are performed to search for all the stable compounds in the Zr-O system. We have found that not only the traditional compound ZrO2, but also the ordered suboxides R3¯-Zr6O, R3¯c-Zr3O, P3¯1m-Zr2O and P6¯2m-ZrO are stable at zero pressure. The crystal structure of semimetallic P6¯2m-ZrO consists of Zr-graphene layers and can be described as an intercalated version of the ω-Zr structure. An interesting massive Dirac cone is found in the three-dimensional (3D) band structure of P6¯2m-ZrO at the F-point. The elastic properties, the hardness and the correlation between the mechanical properties of Zr-O compounds and the oxygen content have been systematically investigated. Surprisingly, the hardest zirconium oxide is not ZrO2, but ZrO. Both P6¯2m-ZrO and P3¯1m-Zr2O exhibit relatively high hardness values of 14 GPa and 10 GPa, respectively. The anisotropic Young's modulus E, torsion shear modulus Gt and linear compressibility β have been derived for P6¯2m-ZrO and P3¯1m-Zr2O. Further analyses of the density of states, the band structure and the crystal orbital Hamilton population indicate that the electronic structure of Zr-O compounds is directly related to their mechanical properties. The simultaneous occurrence of the 3D-framework of Zr-O and the strong Zr-Zr bonds in P6¯2m-ZrO explains its high hardness.