We present the first lattice QCD calculation of the charm quark contribution to the nucleon electromagnetic form factors $G^c_{E,M}(Q^2)$ in the momentum transfer range $0leq Q^2 leq 1.4$ $rm GeV^2$. The quark mass dependence, finite lattice spacing and volume corrections are taken into account simultaneously based on the calculation on three gauge ensembles including one at the physical pion mass. The nonzero value of the charm magnetic moment $mu^c_M=-0.00127(38)_{rm stat}(5)_{rm sys}$, as well as the Pauli form factor, reflects a nontrivial role of the charm sea in the nucleon spin structure. The nonzero $G^c_{E}(Q^2)$ indicates the existence of a nonvanishing asymmetric charm-anticharm sea in the nucleon. Performing a nonperturbative analysis based on holographic QCD and the generalized Veneziano model, we study the constraints on the $[c(x)-bar{c}(x)]$ distribution from the lattice QCD results presented here. Our results provide complementary information and motivation for more detailed studies of physical observables that are sensitive to intrinsic charm and for future global analyses of parton distributions including asymmetric charm-anticharm distribution.