Should we always minimize squared deviations if we want to find the dependency of mean on features? It is known that for any distribution a minimization of squared deviations from a value $c$ gives mean of the distribution. In other words, if we generate many values $y_i$ using our distribution, and then try to find $c$ that minimizes the following expression:
$
\lim_{N \to \infty}  \frac{1}{N}\sum_{i}^N (y_i - c)^2
$
then we find that the minimum is achieved when $c = \mu$, where $\mu$ is the mean of the given distribution. And for any $N$ the above minimization gives the sample mean.

Now, let's assume tat we have a distribution which depends on a single real-valued feature $x$, and we have a data-set generated by this distribution: ($x_1$, $y_1$), ($x_2$, $y_2$), $\dots$, ($x_n$, $y_n$). Then we try to find a function that gives dependency of the mean of this distribution on the feature:
$
\mu = \mu(x, p_1, p_2, \dots, p_k)
$
where $p_i$ are model parameters that we want to find.
Obviously, that we can get the desired function (values of the parameters) if we choose them by minimization of squared deviations:
$
\frac{1}{N} \sum_{i}^{N} [y_i - \mu(x_i, p_1, p_2, \dots, p_k)]^2
$
In other words, we can find the correct answer (correct values of the parameters $p_i$) under condition that $N$ goes to infinity and that the real function, that we are searching, is covered by the space of functions that we consider.
To summarize, the minimization of squared deviations can give not only mean of a given constant distribution, it can also gives correct dependency of the mean on a feature.

In reality $N$ is never infinity and can be rather small. As a result minimization of squared deviation will not give exactly correct values of the model parameters $p_i$. For another sample of the same size we will give other estimations of the parameters $p_i$. In general case, we can say that our estimations of the parameters are somehow distributed.
If we, instead of squared deviations minimize some other measure, we will get other distributions of our model parameters $p_i$. So, different measure of accuracy give different distributions in the parameters space for a given fixed sample size.
These distributions are characterized by their "width" and possible a "shift" relative to the correct point (true values of the model parameters).
Can it be the case that a use of something different from squared deviation gives a better estimate of model parameters (for example more accurate (smaller width) and with smaller or no systematic from the correct model parameters?

Yes, it is possible for estimators obtained by minimizing some different than squared deviation to give a better estimator of model parameters.
The question of whether a given estimator can be beaten by others is studied in statistical decision theory. I'll lay out the basics below then give two examples.
A framework to compare estimators
Suppose we have a data vector $Z \in \mathbb{R}^{p}$, which gives the covariates and response of a single case. Suppose we are estimating a parameter $\mathbf{p}$. Using i.i.d. data $Z_1, \dots, Z_n$, we can estimate the parameter using $\delta(Z_1, \dots, Z_n)$ for some functional $\delta$. We can evaluate the closeness of an estimate to the parameter using the squared-error loss function $\|\delta(Z_1, \dots, Z_n) - \mathbf{p}\|^2$.
Notice the loss function depends on the data as well as the parameter. We can form the risk $R(\delta, \mathbf{p}) = \mathbb{E} \|\delta(Z_1, \dots, Z_n) - \mathbf{p}\|^2$ by averaging over the data. For a given estimation rule $\delta$, the risk $R$ tells us how close we can expect the estimator to be to the truth. Lower is better. Using a simple computation, we can prove that $$R(\delta, \mathbf{p}) = \| \mathrm{Bias}\, \delta(Z_1, \dots, Z_n) \|^2 + \mathrm{trace} \, \mathrm{Var} \, \delta(Z_1, \dots, Z_n),$$ i.e. that the risk trades off the bias (shift) and the variance (width) of the estimator. This decomposition shows that minimizing the variance among unbiased estimators does not lead to the estimator which is expected to be closest to the parameter - instead this tradeoff must be minimized.
We can compare the quality of different estimators $\delta_1$ and $\delta_2$ by comparing the risk curves $\mathbf{p} \mapsto R(\delta_1, \mathbf{p})$ and $\mathbf{p} \mapsto R(\delta_2, \mathbf{p})$. For example, if $R(\delta_1, \mathbf{p}) \leq R(\delta_2, \mathbf{p})$ for all parameters $\mathbf{p}$, this means that the estimation rule $\delta_1$ is will be closer on average to the parameter $\mathbf{p}$ than $\delta_2$ for all possible parameter values. This means that $\delta_1$ dominates $\delta_2$.
A first example
Let's consider the simple example given as the first in OP's question. Here the data $Z=Y \in \mathbb{R}^1$ so that there is only one variable. Let us further assume that $Y = \mu + \epsilon$ for normally distributed $\epsilon$.
The estimator formed by minimizing the empirical squared error is the sample mean $\bar{Y} = \frac{1}{n} \sum_{i=1}^n Y_i$. It is a classical result that the estimator $\delta(Y_1, \dots, Y_n) = \bar{Y}$ is admissible. This means that there does not exist any other estimator $\delta_2(Y_1, \dots, Y_n)$ which dominates the sample mean.
A second example
Now let's consider a linear regression example. Let the data be given by $Z=(y, x)$, where the outcome $y$ is scalar and the covariates $x \in \mathbb{R}^{p-1}$. Assume that $y = x^T \beta + \epsilon$, where $\epsilon$ is normally distributed. Let $\beta$ be the target of inference.
The estimate formed by minimizing the empirical squared error is the OLS estimator $\hat\beta$. When $p > 3$, it turns out this estimator is not admissible: that is, there are other estimators which are always closer on average to the true parameter value $\beta$, regardless of its (unknown) value. A classical example is the James-Stein estimator, which equals $s(Z_1, \dots, Z_n) \hat\beta$ for a suitably chosen data-dependent shrinkage term $s \in (0,1)$.
Conclusion
Basing an estimating equation on the loss function does not necessarily lead to finite sample optimality of the estimator. OP is right to question the basis of the procedure.
