We propose and analyze zeroth-order algorithms for optimization over Riemannian manifolds, where we observe only potentially noisy evaluations of the objective function. Our approach is based on estimating the Riemannian gradient from the objective function evaluations. We consider three settings for the objective function: (i) deterministic and smooth, (ii) stochastic and smooth, and (iii) composition of smooth and non-smooth parts. For each of the setting, we characterize the oracle complexity of our algorithm to obtain appropriately defined notions of $\epsilon$-stationary points. Notably, our complexities are independent of the ambient dimension of the Euclidean space in which the manifold is embedded in, and only depend on the intrinsic dimension of the manifold. As a proof of concept, we demonstrate the applicability of our method to the problem of black-box attacks to deep neural networks, by providing simulation and real-world image data based experimental results.