Hello Bob,
thank you for your answer! I'm going to empirically test the difference between the two cases for sure.
In the meantime, if I run 3dDeconvolve with all my regressors as "signal" and adding "polort", i.e.:

3dDeconvolve -input data.nii.gz  -num_stimts 7 -polort 5 \
-stim_file 1 A.1D -stim_label 1 CO2 \
-stim_file 2 M.par[1] -stim_label 2 Motion1 \
-stim_file 3 M.par[2] -stim_label 3 Motion2 \
-stim_file 4 M.par[3] -stim_label 4 Motion3 \
-stim_file 5 M.par[4] -stim_label 5 Motion4 \
-stim_file 6 M.par[5] -stim_label 6 Motion5 \
-stim_file 7 M.par[6] -stim_label 7 Motion6 \
[...]

I get an output similar to the one in the attached picture.
I have two follow up questions to keep the transatlantic enquiring:
1. The first brick in the output is the "Full R^2" of the model. Does that R^2 consider the whole model (signal+baseline), or only the part of the model attributed to signal?
2. In case the Full R^2 considers only the signal part of the model, would there be a way to include the polynomials in the signal part of the model as well?

Cheers,
Stefano