Homoskedasticity in the dependent variable does not imply homoskedasticity in the error terms of a linear regression Suppose you have $Var(Y_i|X_i)=\sigma^2, \forall i$, so that the conditional variance of $Y$ given $X$ is constant.
However, as long as the conditional expectation function is not linear, this does not imply that the residuals of the linear regression are homoskedastic.


*

*To see this is sufficient to write (in the case of one regressor plus the constant): 
$$E[(Y_i-\beta_0-\beta_1X_i)^2|X_i]=E\{[(Y_i-E[Y_i|X])+(E[Y_i|X]-\beta_0-\beta_1X_i)]^2|X_i\}=
Var[Y_i|X_i]+(E[Y_i|X_i]-\beta_0-\beta_1X_i)^2$$


which is in general different from $Var[Y_i|X_i]$ as long as $E[Y_i|X_i]$ is not linear.


*

*However my concern is the following: whatever is the functional form of the CEF, one can always write $Y_i=\beta_0+\beta_1X_i+u_i$, therefore:
$$Var[Y_i|X_i]=Var[(\beta_0+\beta_1X_i+u_i)|X_i]=Var[(\beta_1X_i+u_i)|X_i]\\=
Var[\beta_1X_i|X_i]+Var[u_i|X_i]+2\beta_1Cov[X_i,u_i|X_i]$$
since we have (is it true?) $Var[\beta_1X_i|X_i]=0$ and $2\beta_1Cov[X_i,u_i|X_i]=0$, we conclude that $\sigma^2=Var[Y_i|X_i]=Var[u_i|X_i]$ so that errors must be homoskedastic too.

 A: When you calculate the mean squared error of the residuals from a biased model, you must remember that MSE = variance + bias^2.
You are right you can express a working "mean model" (possibly misspecified) for $Y$ as $E[Y|X] = \beta_0 + \beta_1 X$. I use the expectation notation here, even when the mean model is possibly wrong, because the parameters $\beta_0$ and $\beta_1$ are well defined, and extend from the usual mathematical expression of expectation. 
You can make a somewhat stronger expression than you've stated. Using the working linear model, you can say the following: $Y = \beta_0 + \beta_1 X + e$ with $\mathbb{E}[e] = 0$ and $\mathbb{E}(eX) = 0$. The problem is that these conditions depend on the "design", meaning that we must take expectation over all the values of $X$, which is why I use the empirical DF notation here. For any particular observation, $i$, however it is not necessarily true that $E(e_i) =0$ nor is it true that $\text{Cov}(e_i, X_i)=0$. Suppose the truth is $Y = f(X) + e$, $e \sim_{iid} F$ with $\int_{-\infty}^\infty x F_x = 0$
Note I have used the squared expectation rather than the variance expression because while $X$ is not a random variable, its design impacts the variance. 
Then 
\begin{eqnarray}
E\left|Y- \beta_0 - \beta_1 X\right|^2 &=& E\left|Y- f(X) + f(X) - \beta_0 + \beta_1 X\right|^2 \\
&=&  E\left|Y- f(X)\right|^2 + E\left|f(X) - \beta_0 - \beta_1 X\right|^2\\
\end{eqnarray}
The reason why we must use the squared expectation here to correctly account for the addition to the residual error is the following: at each instance of $X_i$, the observed value of $Y_i$ will vary about a constant but that constant is not given by $\beta_0 + \beta_1 X_i$. So each computed residual will be offset by $f(X) - \beta_0 - \beta_1$ and this squared value will be a contribution to the residual. Put in other ways: the mean squared error will be the squared bias $E\left|f(X) - \beta_0 - \beta_1X\right|^2$ plus the variance $E\left|Y - f(X)\right|^2$.
