We investigate air entrainment and bubble statistics in three-dimensional breaking waves through novel direct numerical simulations of the two-phase air-water flow, resolving the length scales relevant for the bubble formation problem, the capillary length and the Hinze scale. The dissipation due to breaking is found to be in good agreement with previous experimental observations and inertial scaling arguments. The air entrainment properties and bubble size statistics are investigated for various initial characteristic wave slopes. For radii larger than the Hinze scale, the bubble size distribution, can be described by N(r,t)=B(V-0/2(pi))(epsilon(t- Delta tau)/Wg)r(-10/3) r(m)(-2/3) during the active breaking stages, where epsilon(t-Delta tau) is the time-dependent turbulent dissipation rate, with Delta tau the collapse time of the initial air pocket entrained by the breaking wave, W a weighted vertical velocity of the bubble plume, r(m) the maximum bubble radius, g gravity, V-0 the initial volume of air entrained, r the bubble radius and B a dimensionless constant. The active breaking time-averaged bubble size distribution is described by (N) over bar (r)=B(1/2 pi)(epsilon L-l(c)/Wg rho)r(-10/3)r(m)(-2/3), where epsilon(l) is the wave dissipation rate per unit length of breaking crest, rho the water density and L-c the length of breaking crest. Finally, the averaged total volume of entrained air, (V) over bar, per breaking event can be simply related to epsilon(l) by (V) over bar = B(epsilon L-l(c)/Wg rho), which leads to a relationship for a characteristic slope, S, of (V) over bar proportional to S-5/2. We propose a phenomenological turbulent bubble break-up model based on earlier models and the balance between mechanical dissipation and work done against buoyancy forces. The model is consistent with the numerical results and existing experimental results.