I will assume that you already understand the basic principles of multiple regression. In particular, I assume you already understand that all the regression coefficients will change when you add a new predictor variable to a regression, unless the new variable is uncorrelated with the old variables. So I won't spent time on that and will concentrate specifically on the type of interaction model that you have fitted.
The basic explanation here is that interaction effects are non-linear in the original variables and it is always possible for non-linear effects to average out to zero on the scale of the original linear variables. For this reason, it is quite common for an interaction to be significant even when neither (or only one) of the interacting variables were significant without the interaction.
Any variation on this scenario is possible. Even if the interaction isn't quite strong enough to be statistically significant itself, it nevertheless can easily soak up enough variability to push a previously borderline main effect to become significant.
In your case the interaction term is simply a new covariate equal to the product of the original variables. In other words, you added a third predictor variable equal to income_thousands $\times$ edu_uni. This new covariate is obviously nonlinear in the original variables but nevertheless also correlated with the original variables. The new product covariate just fails to make 0.05 significance (p = 0.06), but is nevertheless strong enough to push edu_uni from borderline significant (p = 0.10) to significance (p = 0.018).
One possible explanation in your case is that the original linear regression might be underestimating slightly the obesity for cities that are simultaneously above average for both income and education. The original linear regression predicts a reduction in obesity with both income and education, but this trend might perhaps taper off for cities with both variables high, and that would lead to this sort of underestimate. It might be that either income or education need to be high to get the obesity reduction but that having both high doesn't give a further reduction. In this scenario, the interaction term would get a positive coefficient to pick up this effect at high levels of both variables, but then the coefficients of the original variables would have to become more negative to compensate for the collinear interaction term in the main body of the regression. The original variables would therefore become more significant. The results you show are consistent with this scenario.
Note that this is an essentially nonlinear phenomenon. If obesity depended on income and education in a purely linear fashion, then these effects are not possible.
It is not true that you need to include the interaction term in the model simply because income and education are themselves correlated. The multiple regression arithmetic already accounts for correlations. In fact the type of model you have fitted:
$$E(y) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3 x_1x_2$$
is one that generally should be avoided when $x_1$ and $x_2$ are continuous covariates. It makes sense to try nonlinear terms, but including the product term without either of the square terms very seldom makes sense.
To check whether $x_1$ and $x_2$ have a nonlinear relationship with $y$, it would be reasonable to try a quadratic response surface. But to fit that, you need the quadratic terms as well:
$$E(y) = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3 x_1x_2+\beta_4x_1^2 +\beta_5x_2^2$$
You could then use an anova approach to test whether the quadratic regression is a worthwhile improvement over the linear regression.
The 3-variable regression that you've fitted is a sort of half-way house between the linear and quadratic regressions without any clear interpretation. To see what is wrong with the 3-variable regression, suppose that you redefined your income and education variables by subtracting off the country-wide averages. The new income and education variables would be the same as the old, just shifted by a constant value. Doing so would make no difference to the linear regression -- only the intercept would change. It would also make no difference to the quadratic regression in that fitted values should be unchanged and the p-values for the quadratic terms would stay the same. However your 3-variable regression would change completely, giving different fitted values from before. The trouble is that the 3-variable regression is too special: it isn't invariant under sensible rescalings of the predictor variables and so doesn't have any general interpretation.
When you fit in interaction model
~A*B in R with factors A and B, the interaction model takes into account all possible combinations of the factors A and B. Unfortunately the
A*B formula doesn't do this when A and B are both continuous covariates. In that case
A*B just adds a single product term, which accounts for only a small part of the possible bivariate effects. Despite the terminology in R, this should not truly be considered to fully represent possible ``interactions'' between the covariates.