TL;DR - Your current results are incorrect and unsafe to use. Carry out the analysis again, including all load in the P-delta analysis, not just the axial load.
If you apply axial load "P", alone, to a perfectly straight beam, then there is no "delta" for the P-delta effect to occur.
There are various methods for creating an initial "delta"; I have never used SAP 2000 so have no idea whether there is a built-in method, or whether you have to do it manually to some extent. In software I have used, I would carry out an eigenvalue analysis, pick the first mode-shape, and use it with a code-specified magnitude for imperfections as the initial shape for my P-delta analysis.
The issue with this is selecting an appropriate mode shape. If you pick a hogging mode shape, then P-delta will tend to increase the hogging; if you pick a sagging mode shape then it will increase sagging. Numerically, for a perfectly straight beam, the software could easily choose either hogging or sagging, as they are both as numerically valid as each other.
We now come to your case. An important thing to learn about non-linear analysis is that you cannot carry out linear superposition - it isn't valid. You need to have all loads in the non-linear case.
Going back through the previous methodology, if you carry out an eigenvalue analysis with axial load and with transverse load causing sagging, then the mode shape will be a sagging one (caused by the transverse load). In the P-delta analysis, again including all load, the transverse load will increase the delta still further, giving an even greater P-delta effect.
Your current results showing the axial load reducing the bending moment are incorrect and unsafe to use.