We present the results of the first Dalitz plot analysis of the decay D0→K-π+η. The analysis is performed on a data set corresponding to an integrated luminosity of 953 fb-1 collected by the Belle detector at the asymmetric-energy e+e- KEKB collider. The Dalitz plot is well described by a combination of the six resonant decay channels K̄∗(892)0η, K-a0(980)+, K-a2(1320)+, K̄∗(1410)0η, K∗(1680)-π+ and K2∗(1980)-π+, together with Kπ and Kη S-wave components. The decays K∗(1680)-→K-η and K2∗(1980)-→K-η are observed for the first time. We measure ratio of the branching fractions, B(D0→K-π+η)B(D0→K-π+)=0.500±0.002(stat)±0.020(syst)±0.003(BPDG). Using the Dalitz fit result, the ratio B(K∗(1680)→Kη)B(K∗(1680)→Kπ) is measured to be 0.11±0.02(stat)-0.04+0.06(syst)±0.04(BPDG); this is much lower than the theoretical expectations (≈1) made under the assumption that K∗(1680) is a pure 13D1 state. The product branching fraction B(D0→[K2∗(1980)-→K-η]π+)=(2.2-1.9+1.7)×10-4 is determined. In addition, the πη′ contribution to the a0(980)± resonance shape is confirmed with 10.1σ statistical significance using the three-channel Flatté model. We also measure B(D0→K̄∗(892)0η)=(1.41-0.12+0.13)%. This is consistent with, and more precise than, the current world average (1.02±0.30)%, deviates with a significance of more than 3σ from the theoretical predictions of (0.51-0.92)%.