We measured the production cross-sections of (86,87m,87(m+g),88)Y radioisotopes for proton-induced reactions on natural zirconium by using a stacked-foil activation technique in the energy range of 4-40 MeV at the MC50 cyclotron of the Korean Institute of Radiological and Medical Sciences. The activation method and the stacked-foil technique using the high-resolution HPGe gamma spectrometry were applied to determine the production cross-sections. The beam current was determined by the monitor reactions (nat)Cu(p,x)(62)Zn and (27)Al(p,x)(22,24)Na. This work reports new data for the above reactions. The results were also compared with ALICE-IPPE model calculations compiled in MENDL-2P file. This experiment is in partial agreement with model calculations. The Optimum energy range for the production of (86)Y is 40 - 30 MeV. The highest yield of (86)Y amounts to 33.37 MBq/mu Ah over this energy range and the levels of the (87)Y, (87m)Y and (88)Y impurities to 2.3, 26.55 and 0.33 %, respectively at EOB. The (nat)Zr(p,x) process may be a suitable route for the production of the medically important radioisotope (86)Y.