A density functional theory for bulk and surface elastic constants of biaxial nematic liquid crystals is developed. It is based on a functional Taylor expansion of the free energy of a distorted biaxial nematic with respect to the one-particle distribution function. Detailed microscopic expressions for the biaxial elastic constants of bulk and surface deformations are derived by expanding further the distribution functions into symmetry-adapted Wigner matrices. The final expressions depend on generalized orientational order parameters characterizing the biaxial nematic and on expansion coefficients of the direct pair correlation function. The case where the expansions are truncated at the lowest nontrivial order with respect to the momentum index of the Wigner matrices is analyzed in detail. It gives only six-distinct, nonzero bulk elastic constants. The mixed elastic constants, which measure distortions of more than one director, vanish within this approximation. As in the uniaxial case, a splay-bend degeneracy for all directors is apparent. The theory is next applied to the biaxial nematic phase recently studied by Biscarini et al. [Phys. Rev. Lett. 75, 1803 (1995)] providing numerical estimates of biaxial elastic constants for the case of thermodynamically stable biaxial ordering. It is shown that the values of the elastic constants connected with secondary directors are much lower than those associated with the primary one.