We obtained an analytical solution for the effective thermal conductivity of composites composed of orthotropic matrices and spherical inhomogeneities with interfacial thermal resistance using a micromechanics-based homogenization. We derived the closed form of a modified Eshelby tensor as a function of the interfacial thermal resistance. We then predicted the heat flux of a single inhomogeneity in the infinite media based on the modified Eshelby tensor, which was validated against the numerical results obtained from the finite element analysis. Based on the modified Eshelby tensor and the localization tensor accounting for the interfacial resistance, we derived an analytical expression for the effective thermal conductivity tensor for the composites by a mean-field approach called the Mori-Tanaka method. Our analytical prediction matched very well with the effective thermal conductivity obtained from finite element analysis with up to 10% inhomogeneity volume fraction.