Why in Box-Cox method we try to make x and y normally distributed, but that's not an assumption for linear regression? In Sheather's book, it states that 

The Box-Cox procedure aims to find a transformation that makes the transformed variable close to normally distributed.

To be specific:


Also, when x and y are normally distributed, the maximum likelihood estimates of $\beta_0$ and $\beta_1$ are the same as the least squares estimates.
But in simple linear regression, actually we don't assume this to be necessarily true. Why is that? 
(Since based on the picture above, it seems only when x and y are normally distributed will Y on X to be close to linear, which is just the linear regression model)
Also, Box-Cox method's goal is to make X and Y more normally distributed, yet usually when people use this method for data transformation, they actually want to make the errors(or std.residuals) normally distributed. How does these two relate to each other? 
 A: In reality the box-cox transformation finds a transformation that homogenize variance.  And constant variance is really an important assumption!  The comment of @whuber: The Box-Cox transform is a data transformation (usually for positive data) defined by $Y^{(\lambda)}= \frac{y^\lambda - 1}{\lambda}$ (when $\lambda\not=0$ and its limit $\log y$ when $\lambda=0$). This transform can be used in different ways, and the Box-Cox method usually refers to likelihood estimation of the transform parameter $\lambda$. $\lambda$ could potentially be chosen in other ways, but this post (and the question) is about this likelihood method of choosing $\lambda$. 
What happens is that boxcox transform maximizes a likelihood function constructed from a constant variance normal model.  And the main contribution in maximizing that likelihood comes from homogenizing the variance! ( * ) You could construct some similar likelihood function from some other location-scale family (maybe, for example, constructed from $t_{10}$, say) and the constant variance assumption, and it would give similar results.  Or you could construct a boxcox-like criterion function from robust regression, again with constant variance. It would give similar results. (eventually, I want to come back here showing this with some code).
( * ) This shouldn't really be surprising. By drawing a few figures you can convince yourself that changing the scale of a density is a much larger change, influencing density values (that is, likelihood values) much more than just changing the basic form a little, but keeping the scale. 
I once built (with Xlispstat) a slider demonstration showing this convincingly, but what you should do is simply to make some simple examples and you will see this result for yourself. 
What happens is simply that the contribution to the likelihood function from constant variance assumption greatly overshadows changes to the likelihood by small changes to the form of the basic density $f_0$ used to generate the location-scale family.  
A: I'm assuming you're referring to Box-Cox normality plots by "method" in your question. It is true that normality assumption in OLS is not required for the method to be useful. For instance, regardless of the error distribution it will produce the coefficients that are unbiased under certain other conditions.
With that said though normality assumption is not useless. For instance, in small samples without normality assumption you can't say much about the probability distribution of the ceofficients beyond the variance and covariance. With normality assumption you can estimate this probability distribution. On large samples under certain conditions you can do this without normality assumption using central limit theorem. Normality assumption makes maximum likelihood estimation (MLE) to produce the same coefficients as OLS, and share many properties of the estimators in (again) small samples.
finally, many people use Box-Cox transformation not to normalize the data, but to stabilize the variance. Sometimes variance increases for larger levels of the dependent variable. In this case Box-Cox transformation can help make variance uniform across the sample. This is related to the assumption of homoscedasticity in OLS
A: Sorry that my question is a bit unorganized, but one of my question(and most confused part) is why do we want our predictors and responses variable to be symmetric or normally distributed. And after dwelling on this for two days, I think now I quite got the answer. 
Here is what I found is most useful: https://stats.stackexchange.com/a/123252/161581
The core idea is: 

the log-or-power-transformed, more normally distributed variables
  are more likely to fulfill linear regression's assumptions, particularly linearity, homoscedasticity, and normally distributed
  residual.

As for the reason, the quote picture in my question can answer for the linearity part. Or as @Penguin_Knight said, skewed independent variable would have some data points with very high leverage, potentially able to bias the regression slope.
For others, in the link above, there are two pictures(which I copied below) that show how transformation can help to  make the variance of errors more like a constant and residual plot a better look (i.e. no discernible pattern). 

