should I apply low-pass filter when calculating central derivative? Let us say that we have a discrete signal $I_n$, $n=0, 1, 2, ...$. According to Nyquist theorem the maximum frequency for such discretization is $f_{max} = 0.5$.

Now imagine that I want to calculate the derivative $D_n$ for this signal. The simplest approximation is the right hand side derivative

$$D_n = I_{n+1} - I_n$$

But what if I want to use a central derivative instead?

$$D_n = \frac{1}{2}(I_{n+1} - I_{n-1})$$

Just looking at this relation, it seems like I am increasing the sampling time from $1$ to $2$. Does this mean that before calculating central derivative I should apply low-pass filter to kill all the frequencies $f &gt; 0.25$? Or am I being wrong and the simple fact that I am using the spacing = 2 for calculating derivative does not mean that I increase the sampling distance, because I can still calculate this derivative for each value of $n$?

No matter what the correct answer is (yes or no), can you please explain in more detail why? If the low-pass filter should be applied, can you explain if there is a method, how given discrete time relation I can infer to what degree I should smoothen the signal before using such relation?
No, you are not increasing the sampling time.  The sample rate stays the same.  The Central difference derivative filter is simply longer than the First difference derivative filter.  A longer filter just means more delay before you get your result.

No, you should not low pass filter without a clear reason for why it is needed to meet your signal processing objectives.  Your reasoning is flawed in this case, so don't apply a LPF.

Note that the Ideal derivative filter multiplies a signal's spectrum by $j\omega \cdot e^{j\omega \mu}$ in the frequency domain.  In the time domain, this ideal filter is an infinitely long filter with taps:

$$h[n] = \dfrac{1}{n+\mu}\left[\cos(\pi \left[n+\mu\right])-\mathrm{sinc}(n+\mu)\right]$$

Where

$n$ is the filter tap index in $(-\infty, \infty).$

$\mu$ corresponds to an inter-sample shift of the filter taps in the range $[0.0, 1.0]$ samples.

Given the above expression for $h[n]$, note that:

The Central difference is $h[n]$ truncated to 3 taps with $\mu = 0$.

The First difference is a scaled version of $h[n]$ truncated to 2 taps with $\mu = 0.5$.

Update to add filter magnitude responses plot