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?

NO
It is important to keep in mind that an estimator of a parameter can take on many forms. In fact, constants can be estimators! Consequently, we might find that calculating something other than the empirical mean might have desirable properties for estimating the mean.
In this post, I give an example where estimating the conditional median, by minimizing MAE, gives better estimates of the regression parameters than the OLS estimates, in the sense that the estimator is unbiased (as is the case with the OLS estimator) but has lower variance than the OLS estimator.
Another answer here mentions the Gauss-Markov theorem. As Richard Hardy explains in the answer to my linked question, the MAE minimizer is nonlinear. Thus, Gauss-Markov does not apply. It is fine for minimization of MAE to result in an unbiased estimator that has lower-variance than OLS.
EDIT
Another answer of mine shows when minimizing in-sample MAE results in lower out-of-sample MSE than minimizing in-sample MSE.
EDIT 2
Let's check out a simulation and visualization. In the code, I simulate a heavy-tailed $t_{1.1}$ error term. At each iteration, I calculate the OLS coefficients (with lm) and the MAE-minimizing coefficients (with rq). Then I plot those regression lines, along with the true regression line.
library(quantreg)
set.seed(2022)

N <- 50
B <- 100
beta0 <- 2
beta1 <- -3
x <- seq(0, 1, 1/(N - 1))
yhat <- beta0 + (beta1 * x)
Q0 <- Q1 <- L0 <- L1 <- rep(NA, B)
for (i in 1:B){
  
  y <- yhat + rt(N, 1.1)
  L <- lm(y ~ x)
  L0[i] <- summary(L)$coef[1, 1]
  L1[i] <- summary(L)$coef[2, 1]
  
  Q <- quantreg::rq(y ~ x, tau = 0.5)
  Q0[i] <- summary(Q)$coef[1, 1]
  Q1[i] <- summary(Q)$coef[2, 1]
  
}

par(mfrow = c(2, 1))
plot(
  x, 
  yhat,
  type = 'l', 
  lty = 2, 
  ylim = c(min(L0 + L1*x, Q0 + Q1*x), max(L0 + L1*x, Q0 + Q1*x)),
  main = "OLS"
  )
for (i in 1:B){
  lines(x, L0[i] + L1[i] * x, col = 'red')
}
lines(x, yhat, type = 'l', lty = 2)
#
plot(
  x, 
  yhat,
  type = 'l', 
  lty = 2, 
  ylim = c(min(L0 + L1*x, Q0 + Q1*x), max(L0 + L1*x, Q0 + Q1*x)),
  main = "MAE"
)
for (i in 1:B){
  lines(x, Q0[i] + Q1[i] * x, col = 'red')
}
lines(x, yhat, type = 'l', lty = 2)
par(mfrow = c(1, 1))


The red estimates of the conditional means are much more reasonable when absolute loss is minimized.
EDIT 4
Another example is in "classification" problems with discrete outcomes (say binary for now). The typical loss function minimized is log loss ("crossentropy" in some circles), which corresponds to maximum likelihood estimation in logistic regression. Our Frank Harrell has a strong opinion about minimizing this loss function as opposed to minimizing square loss.
$$
\text{Log Loss}\\
L(y, p) = -\dfrac{1}{N} \sum_{i = 1}^N \bigg[
y_i\log(p_i) + (1 - y_i)\log(1 - p_i)
\bigg]
$$
EDIT 5
Finally, there is the James-Stein estimator, which shows that the OLS solution to linear regression is inadmissible for any reasonable sample size for doing regression, despite the Gaussian conditional distribution. That is, even the maximum likelihood estimator is inadmissible due to being dominated by James-Stein.
