Three-dimensional grain growth behaviors under anisotropic (misorientation-dependent) grain boundary energy and mobility are investigated via phase-field simulations. Based on a multi-phase-field model and parallel graphics-processing unit computing on a supercomputer, very large-scale simulations with more than three million grains are achieved, enabling reliable statistical evaluation of anisotropic grain growth. The anisotropic boundary properties are introduced by the classical Read-Shockley and sigmoidal models; the threshold misorientation angle, Delta theta(h), included in these models is used as a quantity to determine the anisotropy strength of the system. Systematic simulations are performed for different Delta theta(h) values, through which the correlations between the anisotropy strength and grain growth characteristics such as grain size and misorientation distributions are examined. The obtained results show that anisotropic grain growth reaches the steady-state regime irrespective of the Delta theta(h) value. However, the kinetics and microstructural morphology during the steady-state growth are largely dependent on Delta theta(h). Furthermore, by comparison with the simulated results, the applicability of analytical grain growth theories to anisotropic systems are tested. The tests reveal that the steady-state microstructure in anisotropic growth cannot be well captured by the existing theories, which is likely because the basic assumptions of the theories do not hold for anisotropic systems.