A novel two-level Newton-Krylov-Schwarz (NKS) solver is constructed and analyzed for implicit time discretizations of the bidomain reaction-diffusion system in three dimensions. This multiscale system describes the bioelectrical activity of the heart by coupling two degenerate parabolic equations with several ordinary differential equations at each point in space. Together with a finite element discretization in space, the proposed NKS Bidomain solver employs an outer inexact Newton iteration to solve the nonlinear finite element system originating at each time step of the decoupled implicit discretization. The Jacobian update during the Newton iteration is solved by a Krylov method employing a two-level overlapping Schwarz preconditioner. A convergence rate estimate is proved for the resulting preconditioned operator, showing that its condition number is independent of the number of subdomains (scalability) and bounded by the ratio of the subdomains characteristic size and the overlap size. This theoretical result is confirmed by several parallel simulations employing up to more than 2000 processors for scaled and standard speedup tests in three dimensions. The results show the scalability of the proposed NKS Bidomain solver in terms of both nonlinear and linear iterations, in both Cartesian slabs and ellipsoidal cardiac domains.