Inconsistencies between different results of region extractions

The rules of including/excluding boundaries of a spectral region during region extraction might be confusing, consider this example:

```
import astropy.units as u
import numpy as np
from specutils import Spectrum1D, SpectralRegion
from specutils.manipulation import extract_region

flux = np.random.randn(500) * u.Jy

spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.AA, flux=flux)
extract_region(spec1d, SpectralRegion(5000 * u.AA, 5100 * u.AA)) # [ 5000.0 Angstrom, ..., 5100.0 Angstrom ]

spec1d = Spectrum1D(spectral_axis=np.arange(0, 500) * u.pix, flux=flux)
extract_region(spec1d, SpectralRegion(0 * u.pix, 100 * u.pix)) # [ 0.0 pix, ..., 99.0 pix ]
```

One would expect that the spectral regions are treated as closed intervals in both cases. Consider another example where the upper bound is placed between the spectral axis values:

```
spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.AA, flux=flux)
extract_region(spec1d, SpectralRegion(5000 * u.AA, 5100.5 * u.AA)) # [ 5000.0 Angstrom, ..., 5100.0 Angstrom ]

spec1d = Spectrum1D(spectral_axis=np.arange(0, 500) * u.pix, flux=flux)
extract_region(spec1d, SpectralRegion(0 * u.pix, 100.5 * u.pix)) # [ 0.0 pix, ..., 100.0 pix ]
```

The resulting spectral axis is changed only in the "pixel" case. Final example:

```
spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.AA, flux=flux)
extract_region(spec1d, SpectralRegion(0 * u.pix, 5499 * u.pix)) # [ 5000.0 Angstrom, ..., 5499.0 Angstrom ]

spec1d = Spectrum1D(spectral_axis=np.arange(5000, 5500) * u.pix, flux=flux)
extract_region(spec1d, SpectralRegion(5000 * u.pix, 5499 * u.pix)) # [ 5000.0 pix, ..., 5498.0 pix ]
```
In the first case, `5499.0 Angstrom` is included in the resulting spectral axis. In the second case, `5499.0 pix` is excluded.

One possible solution would be to treat _all_ spectral regions as closed intervals during region extraction.

Start at specutils.manipulation.extract_region, using the issue examples to compare wavelength and pixel spectral axes. Add regression coverage for boundary inclusion and verify that extraction treats the stated intervals consistently; done means the examples produce the expected spectral axes.

Written by the indexing model from the issue text.