Obtaining the FFT of a data set collected through an accelerometer
From you code it appears you do a low pass butterworth filter:

[b1,a1]= butter(order_of_butterworth_filter,wn,'low');
accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);


Then another low pass moving average filter which you minus from the previous.

accelerationd2_smoothned = smooth(accelerationd2_filtered,36.5);
accelerationd2_DC_shifted = accelerationd2_filtered - accelerationd2_smoothned;


Yes, this second part will remove DC, but I think this is the first thing you should do. That is, use a HIGH pass butterworth filter (with a few more orders perhaps)

cut_off_frequency = 5;
wn = cut_off_frequency/sampling_frequency;
order_of_butterworth_filter = 4;
[b1,a1]= butter(order_of_butterworth_filter,wn,'high');
accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);


Then apply the FFT. If there is two much high frequency noise, you can apply a low pass butterworth filter to the above output. 

As for you analysis, the FFT uses the entire signal to compute. So any frequency magnitudes represent the average over the whole time series. That is, there is no time information in there. Since you have a large sample of animal data, I imagine you would want to see how things change over time, so I suggest instead using the STFT or some other wavelet transform. 

In MATLAB you might get good results out of the box using spectrogram, no pre-filtering required.

window = 1000;
overlap = round(0.9*window);
s = spectrogram(accelerationd2,window,noverlap)