I think the problem is just inaccuracy of the numerical computation of the freuqency response at very small frequencies. Keep in mind that tf_tf, as it's defined, would evaluate to 0/0 at z = exp(j*w) = 1 and that tf_tf as a pole/zero cancellation at z = 1. I don't think (don't know for sure) that bode() tries to deal with pole/zero cancellation. However, i think dcgain does try to deal with that, so it gives the expected answer. The desired, or expected, Bode plot can be obtained by using minreal() on tf_tf to remove the common pole and zero at z = 1
M = 8; %Decimation factor
N = 1; %Order of the filter
fs = 25e6; %Sampling frequency
T = 1/fs; %Sampling period
z = tf('z',T);
gain = 1/(M);
num_tf = 1-z^(-M);
den_tf = 1-z^(-1);
tf_tf = (gain*num_tf/den_tf)^N
mrtf_tf = minreal(tf_tf)
bode(tf_tf,mrtf_tf)
