We perform a numerical simulation of marine controlled-source electromagnetic （MCSEM） responses in a layered anisotropic seabed using integral equation method. Firstly, we use the operator theory to deduce the contraction integral equation, which is always convergent regardless of any parameters, for it satisfies the contraction mapping condition. Then we introduce the subdomain multigrid and quasi-linear approximation to improve the computation efficiency. We also validate the accuracy and efficiency of the algorithm by synthetic numerical modelings. Finally, we simulate the three-dimensional MCSEM responses in the layered anisotropic seabed, and investigate the effects of anisotropy.